Added QMCChem module

devel/qmcchem/EZFIO.cfg

@@ -0,0 +1,6 @@
+[ci_threshold]
+type: Threshold
+doc: Threshold on the CI coefficients as computed in QMCChem
+interface: ezfio,provider,ocaml
+default: 1.e-8
+

devel/qmcchem/NEED

@@ -0,0 +1,4 @@
+determinants
+davidson_undressed
+fci
+zmq

devel/qmcchem/README.rst

@@ -0,0 +1,12 @@
+==============
+QmcChem Module
+==============
+
+For multi-state calculations, to extract state 2 use:
+
+``
+QP_STATE=2 qp_run save_for_qmcchem x.ezfio
+``
+(state 1 is the ground state).
+
+

devel/qmcchem/expand_spindets_qmcchem.irp.f

@@ -0,0 +1,8 @@
+program densify
+  implicit none
+  read_wf = .True.
+  touch read_wf
+  call generate_all_alpha_beta_det_products()
+  call diagonalize_ci
+  call save_wavefunction
+end

devel/qmcchem/pot_ao_pseudo_ints.irp.f

@@ -0,0 +1,131 @@
+
+BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size)  ]
+ implicit none
+ BEGIN_DOC
+! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
+!
+! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
+! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
+ END_DOC
+ ! l,m    : Y(l,m) parameters
+ ! c(3)   : pseudopotential center
+ ! a(3)   : Atomic Orbital center
+ ! n_a(3) : Powers of x,y,z in the Atomic Orbital
+ ! g_a    : Atomic Orbital exponent      
+ ! r      : Distance between the Atomic Orbital center and the considered point
+ double precision, external :: ylm_orb
+ integer :: n_a(3)
+ double precision :: a(3), c(3), g_a
+ integer :: i,j,k,l,m,n,p
+ double precision :: dr, Ulc
+ double precision :: y
+ double precision, allocatable :: r(:)
+
+ allocate (r(pseudo_grid_size))
+ dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+ r(1) = 0.d0
+ do j=2,pseudo_grid_size
+   r(j) = r(j-1) + dr
+ enddo
+
+ ao_pseudo_grid = 0.d0
+ do j=1,pseudo_grid_size
+   do k=1,nucl_num
+     c(1:3) = nucl_coord(k,1:3)
+     do l=0,pseudo_lmax
+       do i=1,ao_num
+         a(1:3) = nucl_coord(ao_nucl(i),1:3)
+         n_a(1:3) = ao_power(i,1:3)
+         do p=1,ao_prim_num(i)
+           g_a = ao_expo_ordered_transp(p,i)
+           do m=-l,l
+             y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
+             ao_pseudo_grid(i,m,l,k,j) = ao_pseudo_grid(i,m,l,k,j) + &
+                 ao_coef_normalized_ordered_transp(p,i)*y
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+ enddo
+ deallocate(r)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size)  ]
+ implicit none
+ BEGIN_DOC
+! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
+!
+! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
+! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
+ END_DOC
+ ! l,m    : Y(l,m) parameters
+ ! c(3)   : pseudopotential center
+ ! a(3)   : Atomic Orbital center
+ ! n_a(3) : Powers of x,y,z in the Atomic Orbital
+ ! g_a    : Atomic Orbital exponent      
+ ! r      : Distance between the Atomic Orbital center and the considered point
+ double precision, external :: ylm_orb
+ integer :: n_a(3)
+ double precision :: a(3), c(3), g_a
+ integer :: i,j,k,l,m,n,p
+ double precision :: dr, Ulc
+ double precision :: y
+ double precision, allocatable :: r(:)
+
+ allocate (r(pseudo_grid_size))
+
+ dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+ r(1) = 0.d0
+ do j=2,pseudo_grid_size
+   r(j) = r(j-1) + dr
+ enddo
+
+ mo_pseudo_grid = 0.d0
+ do n=1,pseudo_grid_size
+   do k=1,nucl_num
+     do l=0,pseudo_lmax
+       do m=-l,l
+         do i=1,ao_num
+           do j=1,mo_num
+             if (dabs(ao_pseudo_grid(i,m,l,k,n)) < 1.e-12) then
+              cycle
+             endif
+             if (dabs(mo_coef(i,j)) < 1.e-8) then
+              cycle
+             endif
+             mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + &
+                ao_pseudo_grid(i,m,l,k,n)  * mo_coef(i,j)
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+ enddo
+ deallocate(r)
+
+END_PROVIDER
+
+double precision function test_pseudo_grid_ao(i,j)
+  implicit none
+  integer, intent(in) :: i,j
+  integer :: k,l,m,n
+  double precision :: r, dr,u
+  dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+
+  test_pseudo_grid_ao = 0.d0
+  r = 0.d0
+  do k=1,pseudo_grid_size
+    do n=1,nucl_num
+      do l = 0,pseudo_lmax
+        u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr
+        do m=-l,l
+          test_pseudo_grid_ao +=  ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u 
+        enddo
+      enddo
+    enddo
+    r = r+dr
+  enddo
+end

devel/qmcchem/pseudo.irp.f

@@ -0,0 +1,8 @@
+subroutine write_pseudopotential
+  implicit none
+  BEGIN_DOC
+!  Write the pseudo_potential into the EZFIO file
+  END_DOC
+  call ezfio_set_pseudo_mo_pseudo_grid(mo_pseudo_grid)
+end
+

devel/qmcchem/qmc_e_curve.irp.f

@@ -0,0 +1,113 @@
+program e_curve
+  use bitmasks
+  implicit none
+  integer :: i,j,k, kk, nab, m, l
+  double precision :: norm, E, hij, num, ci, cj
+  integer, allocatable :: iorder(:)
+  double precision , allocatable :: norm_sort(:)
+  double precision :: e_0(N_states)
+  PROVIDE mo_two_e_integrals_in_map
+
+  nab = n_det_alpha_unique+n_det_beta_unique
+  allocate ( norm_sort(0:nab), iorder(0:nab) )
+
+  double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
+  double precision, allocatable  :: u_0(:,:), v_0(:,:)
+  allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
+  allocate(u_0(N_states,N_det),v_0(N_states,N_det))
+
+
+
+  norm_sort(0) = 0.d0
+  iorder(0) = 0
+  do i=1,n_det_alpha_unique
+   norm_sort(i) = det_alpha_norm(i)
+   iorder(i) = i
+  enddo
+  
+  do i=1,n_det_beta_unique
+   norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+   iorder(i+n_det_alpha_unique) = -i
+  enddo
+  
+  call dsort(norm_sort(1),iorder(1),nab)
+
+  if (.not.read_wf) then
+    stop 'Please set read_wf to true'
+  endif
+
+  PROVIDE psi_bilinear_matrix_values nuclear_repulsion 
+  print *,  ''
+  print *,  '=============================='
+  print *,  'Energies at different cut-offs'
+  print *,  '=============================='
+  print *,  ''
+  print *,  '=========================================================='
+  print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)',  'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
+  print *,  '=========================================================='
+  double precision :: thresh
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  thresh = 1.d-10
+  do j=0,nab
+    i = iorder(j)
+    if (i<0) then
+      do k=1,n_det
+        if (psi_bilinear_matrix_columns(k) == -i) then
+          psi_bilinear_matrix_values(k,1) = 0.d0
+        endif
+      enddo
+    else
+      do k=1,n_det
+        if (psi_bilinear_matrix_rows(k) ==  i) then
+          psi_bilinear_matrix_values(k,1) = 0.d0
+        endif
+      enddo
+    endif
+    if (thresh > norm_sort(j)) then
+      cycle
+    endif
+
+    u_0 = psi_bilinear_matrix_values(1:N_det,1:N_states)
+    v_t = 0.d0
+    s_t = 0.d0
+    call dtranspose(                                                   &
+        u_0,                                                           &
+        size(u_0, 1),                                                  &
+        u_t,                                                           &
+        size(u_t, 1),                                                  &
+        N_det, N_states)
+    call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
+    call dtranspose(                                                   &
+        v_t,                                                           &
+        size(v_t, 1),                                                  &
+        v_0,                                                           &
+        size(v_0, 1),                                                  &
+        N_states, N_det)
+    
+    double precision, external :: u_dot_u, u_dot_v
+    do i=1,N_states
+      e_0(i) = u_dot_v(v_t(1,i),u_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
+    enddo
+
+    m = 0
+    do k=1,n_det
+     if (psi_bilinear_matrix_values(k,1) /= 0.d0) then
+      m = m+1
+     endif
+    enddo
+
+    if (m == 0) then
+       exit
+    endif
+    E = E_0(1) + nuclear_repulsion
+    norm = u_dot_u(u_0(1,1),N_det)
+    print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)',  thresh, m, &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, E
+    thresh = thresh * dsqrt(10.d0)
+  enddo
+  print *,  '=========================================================='
+
+  deallocate (iorder, norm_sort)
+end
+

devel/qmcchem/save_for_qmcchem.irp.f

@@ -0,0 +1,45 @@
+program save_for_qmc
+
+  integer                        :: iunit
+  logical                        :: exists
+  double precision               :: e_ref
+  
+  ! Determinants
+  read_wf = .True.
+  TOUCH read_wf
+  print *,  "N_det = ", N_det
+  call write_spindeterminants
+  
+  ! Reference Energy
+  if (do_pseudo) then
+    call write_pseudopotential
+  endif
+  call system( &
+   'mkdir -p '//trim(ezfio_filename)//'/simulation ;' // &
+   'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/simulation/.version ; ' // &
+   'mkdir -p '//trim(ezfio_filename)//'/properties ;' // &
+   'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/properties/.version ; ' // &
+   'echo T > '//trim(ezfio_filename)//'/properties/e_loc' &
+  )
+  iunit = 13
+  open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
+  call ezfio_has_fci_energy_pt2(exists)
+  if (exists) then
+    call ezfio_get_fci_energy_pt2(e_ref)
+  else
+    call ezfio_has_fci_energy(exists)
+    if (exists) then
+      call ezfio_get_fci_energy(e_ref)
+    else
+      call ezfio_has_hartree_fock_energy(exists)
+      if (exists) then
+        call ezfio_get_hartree_fock_energy(e_ref)
+      else
+        e_ref = 0.d0
+      endif
+    endif
+  endif
+  write(iunit,*) e_ref
+  close(iunit)
+
+end

devel/qmcchem/truncate_wf_qmcchem.irp.f

@@ -0,0 +1,116 @@
+program truncate
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call run
+end
+
+subroutine run
+  use bitmasks
+  implicit none
+  integer :: i,j,k, kk, nab, m, l
+  double precision :: norm, E, hij, num, ci, cj
+  integer, allocatable :: iorder(:)
+  double precision , allocatable :: norm_sort(:)
+  double precision :: e_0(N_states)
+
+  PROVIDE mo_two_e_integrals_in_map H_apply_buffer_allocated
+
+  nab = n_det_alpha_unique+n_det_beta_unique
+  allocate ( norm_sort(0:nab), iorder(0:nab) )
+
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
+  double precision, allocatable  :: u_0(:,:), v_0(:,:)
+  allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
+  allocate(u_0(N_det,N_states),v_0(N_det,N_states))
+
+  norm_sort(0) = 0.d0
+  iorder(0) = 0
+  do i=1,n_det_alpha_unique
+   norm_sort(i) = det_alpha_norm(i)
+   iorder(i) = i
+  enddo
+  
+  do i=1,n_det_beta_unique
+   norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+   iorder(i+n_det_alpha_unique) = -i
+  enddo
+  
+  call dsort(norm_sort(1),iorder(1),nab)
+
+
+  PROVIDE psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns
+  PROVIDE nuclear_repulsion 
+  print *,  ''
+  do j=0,nab
+    i = iorder(j)
+    if (i<0) then
+      !$OMP PARALLEL DO PRIVATE(k)
+      do k=1,n_det
+        if (psi_bilinear_matrix_columns(k) == -i) then
+          do l=1,N_states
+            psi_bilinear_matrix_values(k,l) = 0.d0
+          enddo
+        endif
+      enddo
+      !$OMP END PARALLEL DO
+    else
+      !$OMP PARALLEL DO PRIVATE(k)
+      do k=1,n_det
+        if (psi_bilinear_matrix_rows(k) ==  i) then
+          do l=1,N_states
+            psi_bilinear_matrix_values(k,l) = 0.d0
+          enddo
+        endif
+      enddo
+      !$OMP END PARALLEL DO
+    endif
+    if (ci_threshold <= norm_sort(j)) then
+      exit
+    endif
+  enddo
+
+  m = 0
+  do k=1,n_det
+    if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then
+    m = m+1
+    endif
+  enddo
+
+  do k=1,N_states
+    E = E_0(k) + nuclear_repulsion
+  enddo
+  print *,  'Number of determinants:', m
+  call wf_of_psi_bilinear_matrix(.True.)
+  call save_wavefunction
+
+  u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
+  v_0(1:N_det,1:N_states) = 0.d0
+  u_t(1:N_states,1:N_det) = 0.d0
+  v_t(1:N_states,1:N_det) = 0.d0
+  s_t(1:N_states,1:N_det) = 0.d0
+  call dtranspose(                                                   &
+      u_0,                                                           &
+      size(u_0, 1),                                                  &
+      u_t,                                                           &
+      size(u_t, 1),                                                  &
+      N_det, N_states)
+  print *, 'Computing H|Psi> ...'
+  call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
+  print *, 'Done'
+  call dtranspose(                                                   &
+      v_t,                                                           &
+      size(v_t, 1),                                                  &
+      v_0,                                                           &
+      size(v_0, 1),                                                  &
+      N_states, N_det)
+  
+  double precision, external :: u_dot_u, u_dot_v
+  do i=1,N_states
+    e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
+    print *,  'E(',i,') = ', e_0(i) + nuclear_repulsion
+  enddo
+
+  deallocate (iorder, norm_sort)
+end
+