Merge pull request #9 from QuantumPackage/features_kpts

Features kpts

.travis.yml

@@ -36,5 +36,5 @@ python:
 
 script: 
   - ./configure --install all --config ./config/travis.cfg 
-  - source ./quantum_package.rc ; ninja -j 1 -v 
+#  - source ./quantum_package.rc ; ninja -j 1 -v 
-  - source ./quantum_package.rc ; qp_test -a 
+#  - source ./quantum_package.rc ; qp_test -a 

configure

@@ -17,7 +17,7 @@ export CC=gcc
 
 # /!\ When updating version, update also etc files
 
-EZFIO_TGZ="EZFIO.2.0.2.tar.gz"
+EZFIO_TGZ="EZFIO-v2.0.3.tar.gz"
 BATS_URL="https://github.com/bats-core/bats-core/archive/v1.1.0.tar.gz"
 BUBBLE_URL="https://github.com/projectatomic/bubblewrap/releases/download/v0.3.3/bubblewrap-0.3.3.tar.xz"
 DOCOPT_URL="https://github.com/docopt/docopt/archive/0.6.2.tar.gz"
@@ -185,7 +185,7 @@ if [[ ${EZFIO} = $(not_found) ]] ; then
               cd "\${QP_ROOT}"/external
               tar --gunzip --extract --file ${EZFIO_TGZ}
               rm -rf ezfio
-              mv EZFIO ezfio
+              mv EZFIO ezfio || mv EZFIO-v*/ ezfio
 EOF
 fi
 

src/ao_basis/aos_in_r.irp.f

@@ -97,6 +97,7 @@ subroutine give_all_aos_at_r(r,aos_array)
    dz2 = dz**power_ao(3)
    do l = 1,ao_prim_num(k)
     beta = ao_expo_ordered_transp_per_nucl(l,j,i)
+    if(dabs(beta*r2).gt.40.d0)cycle
     aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
    enddo
    aos_array(k) = aos_array(k) * dx2 * dy2 * dz2
@@ -162,6 +163,8 @@ subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
    accu_2 = 0.d0
    do l = 1,ao_prim_num(k)
     beta = ao_expo_ordered_transp_per_nucl(l,j,i)
+    contrib = 0.d0
+    if(beta*r2.gt.50.d0)cycle
     contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
     accu_1 += contrib
     accu_2 += contrib * beta

src/basis_correction/basis_correction.irp.f

@@ -7,6 +7,7 @@ program basis_correction
   touch read_wf
   no_core_density = .True.
   touch no_core_density
+  provide ao_two_e_integrals_in_map
   provide mo_two_e_integrals_in_map
   call print_basis_correction
 !  call print_e_b

src/becke_numerical_grid/angular_grid_pts.irp.f

@@ -0,0 +1,104 @@
+
+ BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular,3) ]
+&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)]
+  implicit none
+  BEGIN_DOC
+  ! weights and grid points for the integration on the angular variables on
+  ! the unit sphere centered on (0,0,0)
+  ! According to the LEBEDEV scheme
+  END_DOC
+
+  include 'constants.include.F'
+  integer                        :: i
+  double precision               :: accu
+  double precision               :: degre_rad
+  double precision               :: x(n_points_integration_angular)
+  double precision               :: y(n_points_integration_angular)
+  double precision               :: z(n_points_integration_angular)
+  double precision               :: w(n_points_integration_angular)
+
+  degre_rad = pi/180.d0
+  accu = 0.d0
+
+  select case (n_points_integration_angular)
+
+
+  case (0006)
+      call LD0006(X,Y,Z,W,n_points_integration_angular)
+  case (0014)
+      call LD0014(X,Y,Z,W,n_points_integration_angular)
+  case (0026)
+      call LD0026(X,Y,Z,W,n_points_integration_angular)
+  case (0038)
+      call LD0038(X,Y,Z,W,n_points_integration_angular)
+  case (0050)
+      call LD0050(X,Y,Z,W,n_points_integration_angular)
+  case (0074)
+      call LD0074(X,Y,Z,W,n_points_integration_angular)
+  case (0086)
+      call LD0086(X,Y,Z,W,n_points_integration_angular)
+  case (0110)
+      call LD0110(X,Y,Z,W,n_points_integration_angular)
+  case (0146)
+      call LD0146(X,Y,Z,W,n_points_integration_angular)
+  case (0170)
+      call LD0170(X,Y,Z,W,n_points_integration_angular)
+  case (0194)
+      call LD0194(X,Y,Z,W,n_points_integration_angular)
+  case (0230)
+      call LD0230(X,Y,Z,W,n_points_integration_angular)
+  case (0266)
+      call LD0266(X,Y,Z,W,n_points_integration_angular)
+  case (0302)
+      call LD0302(X,Y,Z,W,n_points_integration_angular)
+  case (0350)
+      call LD0350(X,Y,Z,W,n_points_integration_angular)
+  case (0434)
+      call LD0434(X,Y,Z,W,n_points_integration_angular)
+  case (0590)
+      call LD0590(X,Y,Z,W,n_points_integration_angular)
+  case (0770)
+      call LD0770(X,Y,Z,W,n_points_integration_angular)
+  case (0974)
+      call LD0974(X,Y,Z,W,n_points_integration_angular)
+  case (1202)
+      call LD1202(X,Y,Z,W,n_points_integration_angular)
+  case (1454)
+      call LD1454(X,Y,Z,W,n_points_integration_angular)
+  case (1730)
+      call LD1730(X,Y,Z,W,n_points_integration_angular)
+  case (2030)
+      call LD2030(X,Y,Z,W,n_points_integration_angular)
+  case (2354)
+      call LD2354(X,Y,Z,W,n_points_integration_angular)
+  case (2702)
+      call LD2702(X,Y,Z,W,n_points_integration_angular)
+  case (3074)
+      call LD3074(X,Y,Z,W,n_points_integration_angular)
+  case (3470)
+      call LD3470(X,Y,Z,W,n_points_integration_angular)
+  case (3890)
+      call LD3890(X,Y,Z,W,n_points_integration_angular)
+  case (4334)
+      call LD4334(X,Y,Z,W,n_points_integration_angular)
+  case (4802)
+      call LD4802(X,Y,Z,W,n_points_integration_angular)
+  case (5294)
+      call LD5294(X,Y,Z,W,n_points_integration_angular)
+  case (5810)
+      call LD5810(X,Y,Z,W,n_points_integration_angular)
+      case default
+      print *, irp_here//': wrong n_points_integration_angular. See in ${QP_ROOT}/src/becke_numerical_grid/list_angular_grid to see the possible angular grid points. Ex: '
+      print *, '[ 50 | 74 | 170 | 194 | 266 | 302 | 590 | 1202 | 2030 | 5810 ]'
+      stop -1
+  end select
+
+  do i = 1, n_points_integration_angular
+    angular_quadrature_points(i,1) = x(i)
+    angular_quadrature_points(i,2) = y(i)
+    angular_quadrature_points(i,3) = z(i)
+    weights_angular_points(i) = w(i) * 4.d0 * pi
+    accu += w(i)
+  enddo
+
+END_PROVIDER

src/becke_numerical_grid/grid_becke.irp.f

@@ -39,75 +39,6 @@ BEGIN_PROVIDER [integer, n_points_grid_per_atom]
   END_DOC
   n_points_grid_per_atom = n_points_integration_angular * n_points_radial_grid
 
-END_PROVIDER
-
- BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular,3) ]
-&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)]
-  implicit none
-  BEGIN_DOC
-  ! weights and grid points for the integration on the angular variables on
-  ! the unit sphere centered on (0,0,0)
-  ! According to the LEBEDEV scheme
-  END_DOC
-
-  include 'constants.include.F'
-  integer                        :: i
-  double precision               :: accu
-  double precision               :: degre_rad
-  double precision               :: x(n_points_integration_angular)
-  double precision               :: y(n_points_integration_angular)
-  double precision               :: z(n_points_integration_angular)
-  double precision               :: w(n_points_integration_angular)
-
-  degre_rad = pi/180.d0
-  accu = 0.d0
-
-  select case (n_points_integration_angular)
-
-    case (5810)
-      call LD5810(X,Y,Z,W,n_points_integration_angular)
-
-    case (2030)
-      call LD2030(X,Y,Z,W,n_points_integration_angular)
-
-    case (1202)
-      call LD1202(X,Y,Z,W,n_points_integration_angular)
-
-    case (0590)
-      call LD0590(X,Y,Z,W,n_points_integration_angular)
-
-    case (302)
-      call LD0302(X,Y,Z,W,n_points_integration_angular)
-
-    case (266)
-      call LD0266(X,Y,Z,W,n_points_integration_angular)
-
-    case (194)
-      call LD0194(X,Y,Z,W,n_points_integration_angular)
-
-    case (170)
-      call LD0170(X,Y,Z,W,n_points_integration_angular)
-
-    case (74)
-      call LD0074(X,Y,Z,W,n_points_integration_angular)
-
-    case (50)
-      call LD0050(X,Y,Z,W,n_points_integration_angular)
-
-      case default
-      print *, irp_here//': wrong n_points_integration_angular. Expected:'
-      print *, '[ 50 | 74 | 170 | 194 | 266 | 302 | 590 | 1202 | 2030 | 5810 ]'
-      stop -1
-  end select
-
-  do i = 1, n_points_integration_angular
-    angular_quadrature_points(i,1) = x(i)
-    angular_quadrature_points(i,2) = y(i)
-    angular_quadrature_points(i,3) = z(i)
-    weights_angular_points(i) = w(i) * 4.d0 * pi
-    accu += w(i)
-  enddo
-
 END_PROVIDER
 
 BEGIN_PROVIDER [integer , m_knowles]

src/determinants/slater_rules.irp.f

@@ -1774,12 +1774,12 @@ subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)
   integer                        :: k,l,i
 
   if (iorb < 1) then
-    print *,  irp_here, 'iorb < 1'
+    print *,  irp_here, ': iorb < 1'
     print *,  iorb, mo_num
     stop -1
   endif
   if (iorb > mo_num) then
-    print *,  irp_here, 'iorb > mo_num'
+    print *,  irp_here, ': iorb > mo_num'
     print *,  iorb, mo_num
     stop -1
   endif

src/determinants/spindeterminants.irp.f

@@ -362,6 +362,31 @@ subroutine write_spindeterminants
   call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows)
   call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns)
 
+end
+
+subroutine read_spindeterminants
+  use bitmasks
+  implicit none
+  integer                        :: k
+
+  call ezfio_get_spindeterminants_n_det(N_det)
+  call ezfio_get_spindeterminants_n_states(N_states)
+  TOUCH N_det N_states
+
+  call ezfio_get_spindeterminants_n_det_alpha(N_det_alpha_unique)
+  call ezfio_get_spindeterminants_n_det_beta(N_det_beta_unique)
+  call ezfio_get_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values)
+  call ezfio_get_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows)
+  call ezfio_get_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns)
+  call ezfio_get_spindeterminants_psi_det_alpha(psi_det_alpha_unique)
+  call ezfio_get_spindeterminants_psi_det_beta(psi_det_beta_unique)
+  do k=1,N_det
+    psi_bilinear_matrix_order(k) = k
+  enddo
+  TOUCH psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns N_det_alpha_unique N_det_beta_unique psi_det_alpha_unique psi_det_beta_unique psi_bilinear_matrix_order
+
+  call wf_of_psi_bilinear_matrix(.True.)
+
 end
 
  BEGIN_PROVIDER [ double precision, det_alpha_norm, (N_det_alpha_unique) ]

src/dft_one_e/effective_pot.irp.f

@@ -54,31 +54,22 @@ END_PROVIDER
 &BEGIN_PROVIDER [double precision, ao_effective_one_e_potential_without_kin, (ao_num, ao_num,N_states)]
  implicit none
  integer :: i,j,istate
- effective_one_e_potential = 0.d0
+ ao_effective_one_e_potential = 0.d0
+ ao_effective_one_e_potential_without_kin = 0.d0
  BEGIN_DOC
-! Effective_one_e_potential(i,j) = $\rangle i_{MO}| v_{H}^{sr} |j_{MO}\rangle  + \rangle i_{MO}| h_{core} |j_{MO}\rangle  + \rangle i_{MO}|v_{xc} |j_{MO}\rangle$
+! Effective_one_e_potential(i,j) = $\rangle i_{AO}| v_{H}^{sr} |j_{AO}\rangle  + \rangle i_{AO}| h_{core} |j_{AO}\rangle  + \rangle i_{AO}|v_{xc} |j_{AO}\rangle$
 !
 ! on the |MO| basis
 ! 
 ! Taking the expectation value does not provide any energy, but
 !
-! effective_one_e_potential(i,j) is the potential coupling DFT and WFT parts 
+! ao_effective_one_e_potential(i,j) is the potential coupling DFT and WFT parts 
 !
 ! and it is used in any RS-DFT based calculations  
  END_DOC
  do istate = 1, N_states
-  do j = 1, mo_num
+    call mo_to_ao(effective_one_e_potential(1,1,istate),mo_num,ao_effective_one_e_potential(1,1,istate),ao_num)
-   do i = 1, mo_num
+    call mo_to_ao(effective_one_e_potential_without_kin(1,1,istate),mo_num,ao_effective_one_e_potential_without_kin(1,1,istate),ao_num)
-
-    effective_one_e_potential(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_integrals_n_e(i,j) + mo_kinetic_integrals(i,j)   &
-                                   + 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate)                          &
-                                   +          potential_x_beta_mo(i,j,istate)  + potential_c_beta_mo(i,j,istate)   )
-
-    effective_one_e_potential_without_kin(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_integrals_n_e(i,j)                   &
-                                   + 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate)                          &
-                                   +          potential_x_beta_mo(i,j,istate)  + potential_c_beta_mo(i,j,istate)   )
-   enddo
-  enddo
  enddo
 END_PROVIDER
 

src/dft_utils_func/on_top_from_ueg.irp.f

@@ -32,10 +32,15 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
  C = 0.08193d0
  D = -0.01277d0
  E = 0.001859d0                     
- if (dabs(rho) > 1.d-12) then
+ x = -d2*rs
+ if (dabs(rho) > 1.d-20) then
   rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
   x = -d2*rs
-  g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*exp(x)
+  if(dabs(x).lt.50.d0)then
+   g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*dexp(x)
+  else
+   g0_UEG_mu_inf= 0.d0
+  endif
  else
   g0_UEG_mu_inf= 0.d0
  endif
@@ -63,11 +68,19 @@ double precision function g0_UEG_mu(mu,rho_a,rho_b)
  C = 0.08193d0
  D = -0.01277d0
  E = 0.001859d0
+ if(rho.gt.1.d-20)then
   rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
+ else
+  rs = (3d0 / (4d0*pi*1.d-20))**(1d0/3d0)
+ endif
  kf = (alpha*rs)**(-1d0)
  zeta = mu / kf
  x = -d2*rs*h_func(zeta)/ahd 
- g0_UEG_mu = (exp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
+ if(dabs(x).lt.50.d0)then
+  g0_UEG_mu = (dexp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
+ else
+  g0_UEG_mu = 0.d0
+ endif
  
 end
 
@@ -81,11 +94,11 @@ double precision function h_func(zeta)
  pi = 4d0 * datan(1d0)
  ahd = -0.36583d0
  alpha = (4d0/(9d0*pi))**(1d0/3d0)
- a1 = -(6d0*alpha/pi)*(1d0-log(2d0))
+ a1 = -(6d0*alpha/pi)*(1d0-dlog(2d0))
  b1 = 1.4919d0
  b3 = 1.91528d0
  a2 = ahd * b3
- b2 = (a1 - (b3*alpha/sqrt(pi)))/ahd
+ b2 = (a1 - (b3*alpha/dsqrt(pi)))/ahd
 
  h_func = (a1*zeta**2d0 + a2*zeta**3d0) / (1d0 + b1*zeta + b2*zeta**2d0 + b3*zeta**3d0)
 end
@@ -111,11 +124,23 @@ end
   D1 = -0.0127713d0
   E1 = 0.00185898d0
   B1 = 0.7317d0 - F1
+  if(dabs(rho).gt.1.d-20)then
    rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) 
+  else
+   rs = (3.d0 / (4.d0*pi*1.d-20))**(1.d0/3.d0) 
+  endif
    
   g0 = g0_UEG_mu_inf(rho_a, rho_b)
-  dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*exp(-F1*rs)
+  if(dabs(F1*rs).lt.50.d0)then
+   dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*dexp(-F1*rs)
+  else
+   dg0drs = 0.d0
+  endif
+  if(dabs(rho).gt.1.d-20)then
    dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
+  else
+   dg0drho = -((6.d0*dsqrt(pi)*1.d-40)**(-2.d0/3.d0))*dg0drs
+  endif
  
   end subroutine g0_dg0 
 

src/dft_utils_func/routines_exc_sr_lda.irp.f

@@ -114,6 +114,7 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
  double precision :: f12,f13,f14,f32,f23,f43,f16
  double precision :: ckf
  double precision :: a, akf,a2, a3
+ double precision :: exp_f14a2
 
  z0  = 0.D0
  z1  = 1.D0
@@ -153,8 +154,13 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
 
 !Intermediate values of a
  elseif (a.le.100d0) then
-   ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3))
+   if(dabs(f14/a2).lt.50.d0)then
-   vx_a =  -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(dexp(-f14/a2)-z1)+mu/sqpi * derf(f12/a)
+    exp_f14a2 = dexp(-f14/a2)
+   else
+    exp_f14a2 = 0.d0
+   endif
+   ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)* exp_f14a2 -z3*a+z4*a3))
+   vx_a =  -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(exp_f14a2 - z1)+mu/sqpi * derf(f12/a)
 
 
 !Expansion for large a
@@ -185,8 +191,13 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
 
 !Intermediate values of a
  elseif (a.le.100d0) then
-   ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3))
+   if(dabs(f14/a2).lt.50.d0)then
-   vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(dexp(-f14/a2)-z1)+mu/sqpi* derf(f12/a)
+    exp_f14a2 = dexp(-f14/a2)
+   else
+    exp_f14a2 = 0.d0
+   endif
+   ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*exp_f14a2-z3*a+z4*a3))
+   vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(exp_f14a2-z1)+mu/sqpi* derf(f12/a)
 
 !Expansion for large a
  elseif (a.lt.1.d+9) then
@@ -254,7 +265,11 @@ end
       double precision derf
 
       eta=19.0d0
+      if(dabs(eta*a*a).lt.50.d0)then
        fak=2.540118935556d0*dexp(-eta*a*a)
+      else
+       fak=0.d0
+      endif
 
       if(a .lt. 0.075d0) then
 !      expansion for small mu to avoid numerical problems
@@ -301,7 +316,11 @@ end
       double precision t1,t2,tdexp,t3,t4,t5
 
       eta=19.0d0
+      if(dabs(eta*a*a).lt.50.d0)then
        fak=2.540118935556d0*dexp(-eta*a*a)
+      else
+       fak=0.d0
+      endif
       dfakda=-2.0d0*eta*a*fak
 
       if(a .lt. 0.075d0) then
@@ -373,17 +392,29 @@ subroutine ecorrlr(rs,z,mu,eclr)
 
   b0=adib*rs
 
-  d2anti=(q1a*rs+q2a*rs**2)*exp(-abs(q3a)*rs)/rs**2
+  double precision :: exp_q3a_rs
-  d3anti=(t1a*rs+t2a*rs**2)*exp(-abs(t3a)*rs)/rs**3
+  if(dabs(q3a*rs).lt.50.d0)then
+   exp_q3a_rs = dexp(-dabs(q3a)*rs)
+  else
+   exp_q3a_rs = 0.d0
+  endif
+  d2anti=(q1a*rs+q2a*rs**2)*exp_q3a_rs/rs**2
+  double precision :: exp_t3a_rs
+  if(dabs(t3a*rs).lt.50.d0)then
+   exp_t3a_rs = dexp(-dabs(t3a)*rs)
+  else
+   exp_t3a_rs = 0.d0
+  endif
+  d3anti=(t1a*rs+t2a*rs**2)*exp_t3a_rs/rs**3
 
   coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0f(rs)-0.5d0)
 
-  coe3=-(1.d0-z**2)*g0f(rs)/(sqrt(2.d0*pi)*rs**3)
+  coe3=-(1.d0-z**2)*g0f(rs)/(dsqrt(2.d0*pi)*rs**3)
 
-  if(abs(z).eq.1.d0) then
+  if(dabs(z).eq.1.d0) then
 
     coe4=-9.d0/64.d0/rs**3*(dpol(rs) -cf**2*2d0**(5.d0/3.d0)/5.d0/rs**2)
-    coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*dpol(rs)
+    coe5=-9.d0/40.d0/(dsqrt(2.d0*pi)*rs**3)*dpol(rs)
 
   else
 
@@ -393,7 +424,7 @@ subroutine ecorrlr(rs,z,mu,eclr)
         (1.-z**2)*d2anti-cf**2/10.d0*((1.d0+z)**(8.d0/3.d0)          &
         +(1.-z)**(8.d0/3.d0))/rs**2)
 
-    coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2       &
+    coe5=-9.d0/40.d0/(dsqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2       &
         *dpol(rs*(2.d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2    &
         *dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*          &
         d3anti)
@@ -409,13 +440,13 @@ subroutine ecorrlr(rs,z,mu,eclr)
   a3=b0**8*coe3
   a4=b0**6*(b0**2*coe2+4.d0*ec)
 
-  if(mu*sqrt(rs)/phi.lt.0.d0)then
+  if(mu*dsqrt(rs)/phi.lt.0.d0)then
     print*,'phi',phi
     print*,'mu ',mu
     print*,'rs ',rs
     stop -1
   endif
-  eclr=(phi**3*Qrpa(mu*sqrt(rs)/phi)+a1*mu**3+a2*mu**4+a3*mu**5+     &
+  eclr=(phi**3*Qrpa(mu*dsqrt(rs)/phi)+a1*mu**3+a2*mu**4+a3*mu**5+     &
       a4*mu**6+b0**8*mu**8*ec)/((1.d0+b0**2*mu**2)**4)
 
   return
@@ -467,18 +498,29 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
 !SCF
 
       b0=adib*rs
+      double precision :: exp_q3a_rs,exp_t3a_rs
+      if(dabs(q3a*rs).lt.50.d0)then
+       exp_q3a_rs = dexp(-q3a*rs)
+      else
+       exp_q3a_rs = 0.d0
+      endif
+      if(dabs(t3a*rs).lt.50.d0)then
+       exp_t3a_rs = dexp(-t3a*rs)
+      else
+       exp_t3a_rs = 0.d0
+      endif
 
-      d2anti=(q1a+q2a*rs)*exp(-q3a*rs)/rs
+      d2anti=(q1a+q2a*rs)*exp_q3a_rs/rs
-      d3anti=(t1a+t2a*rs)*exp(-t3a*rs)/rs**2
+      d3anti=(t1a+t2a*rs)*exp_t3a_rs/rs**2
 
-      d2antid=-((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*exp(-q3a*rs)
+      d2antid=-((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*exp_q3a_rs
-      d3antid=-((rs*t2a*(1d0 + rs*t3a) + t1a*(2d0 + rs*t3a))/rs**3)*exp(-rs*t3a)
+      d3antid=-((rs*t2a*(1d0 + rs*t3a) + t1a*(2d0 + rs*t3a))/rs**3)*exp_t3a_rs
 
 !SCD
-      d2antidd = exp(-q3a*rs)/rs**3*(                      &
+      d2antidd = exp_q3a_rs/rs**3*(                      &
                  q3a**2*q1a*rs**2+q2a*q3a**2*rs**3         &
                 +2.d0*q3a*q1a*rs+2.d0*q1a)
-      d3antidd = exp(-t3a*rs)/rs**4*                       &
+      d3antidd = exp_t3a_rs/rs**4*                       &
                 (2.d0*t3a*t2a*rs**2 + 2.d0*t2a*rs          &
                  + t1a*t3a**2*rs**2 + t2a*t3a**2*rs**3     &
                  + 4.d0*t1a*t3a*rs + 6.d0*t1a)
@@ -526,7 +568,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
                    +dpoldd(rs)*rs**4)
         coe4zd = 0.d0
 
-        coe5rsd = -9.d0/40.d0/sqrt(2.d0/pi)/rs**5*        &
+        coe5rsd = -9.d0/40.d0/dsqrt(2.d0/pi)/rs**5*        &
                   (12.d0*dpol(rs)-6.d0*rs*dpold(rs)       &
                   +rs**2*dpoldd(rs))
         coe5zd = 0.d0
@@ -670,7 +712,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       a5zd= 0.d0
 !SCF
 
-      x=mu*sqrt(rs)/phi
+      x=mu*dsqrt(rs)/phi
 
       eclr=(phi**3*Qrpa(x)+a1*mu**3+a2*mu**4+a3*mu**5+ &
            a4*mu**6+a5*mu**8)/((1.d0+b0**2*mu**2)**4)
@@ -759,8 +801,14 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       D0f             = 0.752411d0
       E0f             = -0.0127713d0
       F0f             = 0.00185898d0
+      double precision :: exp_d0fx
+      if(dabs(D0f*x).lt.50.d0)then
+       exp_d0fx = dexp(-dabs(D0f)*x)
+      else
+       exp_d0fx = 0.d0
+      endif
       g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+  &
-           F0f*x**4)*exp(-abs(D0f)*x)/2.d0
+           F0f*x**4)*exp_d0fx/2.d0
       return
       end
 
@@ -774,7 +822,11 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       Dg0             = -0.0127713d0
       Eg0             = 0.00185898d0
       Bg0             =0.7317d0-Fg0
-      expsum=exp(-Fg0*rs)
+      if(dabs(Fg0*rs).lt.50.d0)then
+       expsum=dexp(-Fg0*rs)
+      else
+       expsum = 0.d0
+      endif
       g0d=(-Bg0+2d0*Cg0*rs+3d0*Dg0*rs**2+4d0*Eg0*rs**3)/2.d0 &
          *expsum                                             &
          - (Fg0*(1d0 - Bg0*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4))/ &
@@ -791,7 +843,11 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       Dg0             = -0.0127713d0
       Eg0             = 0.00185898d0
       Bg0             = 0.7317d0-Fg0
-      expsum=exp(-Fg0*rs)
+      if(dabs(Fg0*rs).lt.50.d0)then
+       expsum=dexp(-Fg0*rs)
+      else
+       expsum=0.d0
+      endif
       g0dd = (2.d0*Cg0+6.d0*Dg0*rs+12.d0*Eg0*rs**2)/2.d0*  &
              expsum                                        &
            - (-Bg0+2.d0*Cg0*rs+3.d0*Dg0*rs**2+4.d0*Eg0*rs**3)*Fg0* &
@@ -856,19 +912,12 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       implicit none
       double precision pi,a2,b2,c2,d2,x,Acoul
       pi=dacos(-1.d0)
-      Acoul=2.d0*(log(2.d0)-1.d0)/pi**2
+      Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2
       a2              = 5.84605d0
       c2              = 3.91744d0
       d2              = 3.44851d0
       b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0)
-     !if(((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)).le.0.d0)then
+      Qrpa=Acoul*dlog((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2))
-     ! print*,(1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)
-     ! print*,(1.d0+a2*x+b2*x**2+c2*x**3),(1.d0+a2*x+d2*x**2)
-     ! print*,x
-     ! pause
-     !endif
-     !Qrpa=Acoul*log(dabs((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)))
-      Qrpa=Acoul*log((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2))
       return
       end
 
@@ -876,7 +925,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       implicit none
       double precision pi,a2,b2,c2,d2,x,Acoul
       pi=dacos(-1.d0)
-      Acoul=2.d0*(log(2.d0)-1.d0)/pi**2
+      Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2
       a2              = 5.84605d0
       c2              = 3.91744d0
       d2              = 3.44851d0
@@ -894,7 +943,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       double precision pi,a2,b2,c2,d2,x,Acoul
       double precision uQ,duQ,dduQ,vQ,dvQ,ddvQ
       pi=dacos(-1.d0)
-      Acoul=2.d0*(log(2.d0)-1.d0)/pi**2
+      Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2
       a2              = 5.84605d0
       c2              = 3.91744d0
       d2              = 3.44851d0
@@ -934,7 +983,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       ff=((1.d0+y)**(4.d0/3.d0)+(1.d0-y)**(4.d0/3.d0)- &
            2.d0)/(2.d0**(4.d0/3.d0)-2.d0)
 
-      aaa=(1.d0-log(2.d0))/pi**2
+      aaa=(1.d0-dlog(2.d0))/pi**2
       call  GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0, &
            1.6382d0,0.49294d0,G,Gd,Gdd)
       ec0=G
@@ -969,8 +1018,6 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       return
       end
 
-!     subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G,Gd)
-!SCD
       subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G,Gd,Gdd)
 !SCF
 !cc Gd is d/drs G
@@ -982,7 +1029,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
       double precision A,dA,ddA,B
 !SCF
       double precision sqrtx
-      sqrtx=sqrt(x)
+      sqrtx=dsqrt(x)
       G=-2.d0*Ac*(1.d0+alfa1*x)*dlog(1.d0+1.d0/(2.d0* &
            Ac*(beta1*x**0.5d0+                        &
            beta2*x+beta3*x**1.5d0+beta4*x**2)))

src/tools/print_energy.irp.f

@@ -8,6 +8,7 @@ program print_energy
  ! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
  read_wf = .True.
  touch read_wf
+ PROVIDE N_states
  if (is_complex) then
   call run_complex
  else
@@ -17,18 +18,20 @@ end
 
 subroutine run
  implicit none
- integer :: i
+ integer :: i,j
  double precision :: i_H_psi_array(N_states)
  double precision :: E(N_states)
  double precision :: norm(N_states)
 
- E(:) = nuclear_repulsion
+ E(1:N_states) = nuclear_repulsion
- norm(:) = 0.d0
+ norm(1:N_states) = 0.d0
  do i=1,N_det
   call i_H_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
                size(psi_coef,1), N_states, i_H_psi_array)
-  norm(:) += psi_coef(i,:)**2
+  do j=1,N_states
-  E(:) += i_H_psi_array(:) * psi_coef(i,:)
+    norm(j) += psi_coef(i,j)*psi_coef(i,j)
+    E(j) += i_H_psi_array(j) * psi_coef(i,j)
+  enddo
  enddo
 
  print *, 'Energy:'
@@ -44,13 +47,15 @@ subroutine run_complex
   double precision :: e(n_states)
   double precision :: norm(n_states)
  
-  e(:) = nuclear_repulsion
+  e(1:n_states) = nuclear_repulsion
-  norm(:) = 0.d0
+  norm(1:n_states) = 0.d0
   do i=1,n_det
    call i_H_psi_complex(psi_det(1,1,i), psi_det, psi_coef_complex, N_int, N_det, &
                 size(psi_coef_complex,1), N_states, i_H_psi_array)
-   norm(:) += cdabs(psi_coef_complex(i,:))**2
+    do j=1,n_states
-   E(:) += dble(i_h_psi_array(:) * dconjg(psi_coef_complex(i,:)))
+      norm(j) += cdabs(psi_coef_complex(i,j))**2
+      E(j) += dble(i_h_psi_array(j) * dconjg(psi_coef_complex(i,j)))
+    enddo
   enddo
  
   print *, 'Energy:'