Computing alpha_0 correction

devel/sr_correction/EZFIO.cfg

@@ -0,0 +1,12 @@
+[energy_mu]
+type: double precision
+doc: Energy E(mu)
+size: (determinants.n_states)
+interface: ezfio
+
+[energy]
+type: double precision
+doc: Energy E(mu) + alpha_0 <Psi(mu)|\bar{W}|Psi(mu)>
+size: (determinants.n_states)
+interface: ezfio
+

devel/sr_correction/NEED

@@ -0,0 +1,3 @@
+davidson_undressed
+two_body_rdm
+mo_two_e_ints

devel/sr_correction/README.rst

@@ -0,0 +1,5 @@
+=============
+sr_correction
+=============
+
+Correction of the energy obtained with a Long-Range Hamiltonian

devel/sr_correction/alpha.irp.f

@@ -0,0 +1,18 @@
+BEGIN_PROVIDER [ double precision, alpha_coef, (0:1) ]
+ implicit none
+ BEGIN_DOC
+ ! SavCar-JCP-23
+ END_DOC
+
+ double precision :: num, den
+
+ num = 0.319820d0 + mu_erf * (1.063846d0 + mu_erf)
+ den = 0.487806d0 + mu_erf * (1.375439d0 + mu_erf)
+ alpha_coef(0) = num/den
+
+ num = 0.113074d0 + mu_erf * (0.638308d0 + mu_erf)
+ den = 0.122652d0 + mu_erf * (0.674813d0 + mu_erf)
+ alpha_coef(1) = num/den
+END_PROVIDER
+
+

devel/sr_correction/energy_mu.irp.f

@@ -0,0 +1,80 @@
+BEGIN_PROVIDER [ double precision, W_mu, (mo_num, mo_num, mo_num, mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! $<ij|W(mu)|kl>$ in MO basis
+ END_DOC
+
+ integer :: i,j,k,l
+ do l=1,mo_num
+   do k=1,mo_num
+     do j=1,mo_num
+       do i=1,mo_num
+         double precision, external :: mo_two_e_integral_erf
+         W_mu(i,j,k,l) = mo_two_e_integral_erf(i,j,k,l)
+       enddo
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, W_bar_mu, (mo_num, mo_num, mo_num, mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! $<ij|\bar{W}(mu)|kl>$ in MO basis
+ END_DOC
+
+ integer :: i,j,k,l
+ do l=1,mo_num
+   do k=1,mo_num
+     do j=1,mo_num
+       do i=1,mo_num
+         double precision, external :: mo_two_e_integral
+         W_bar_mu(i,j,k,l) = mo_two_e_integral(i,j,k,l) - W_mu(i,j,k,l)
+       enddo
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, energy_mu ]
+ implicit none
+ BEGIN_DOC
+ ! E(mu)
+ END_DOC
+
+ double precision :: one_e, two_e
+
+ integer :: k,l
+ double precision, external :: ddot
+
+ one_e = 0.d0
+ two_e = 0.d0
+ do l=1,mo_num
+   one_e += ddot(mo_num, one_e_dm_mo(1,l), 1, mo_one_e_integrals(1,l), 1)
+   do k=1,mo_num
+     two_e += 0.5d0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_mu(1,1,k,l), 1)
+   enddo
+ enddo
+ energy_mu = one_e + two_e + nuclear_repulsion
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, correction_alpha_0 ]
+ implicit none
+ BEGIN_DOC
+ ! alpha_0(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>
+ END_DOC
+
+ double precision :: one_e, two_e
+
+ integer :: k,l
+ double precision, external :: ddot
+
+ correction_alpha_0 = 0.d0
+ do l=1,mo_num
+   do k=1,mo_num
+     correction_alpha_0 += 0.5d0 * ddot (mo_num*mo_num, two_e_dm_mo(1,1,k,l), 1, W_bar_mu(1,1,k,l), 1)
+   enddo
+ enddo
+ correction_alpha_0 = correction_alpha_0*alpha_coef(0)
+END_PROVIDER
+

devel/sr_correction/sr_correction.irp.f

@@ -0,0 +1,14 @@
+program sr_correction_opt_psi
+ implicit none
+ read_wf = .True.
+ use_only_lr = .False.
+ SOFT_TOUCH read_wf use_only_lr
+ call run
+end
+
+subroutine run
+ implicit none
+ print *, 'E(mu)', energy_mu
+ print *, 'correction', correction_alpha_0
+ print *, 'E = ', energy_mu + correction_alpha_0
+end

stable/qmcchem/pot_ao_pseudo_ints.irp.f

@@ -89,13 +89,10 @@ BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
      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
+           do j=1,mo_num
-              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