Fixed guess

plugins/local/basis_correction/weak_corr_func.irp.f

@@ -7,7 +7,7 @@
  END_DOC
  implicit none
  integer :: ipoint,istate
- double precision :: rho_a, rho_b, ec
+ double precision :: rho_a, rho_b, ec, rho, p2
  double precision :: wall0,wall1,weight,mu
  logical :: dospin
  dospin = .true. ! JT dospin have to be set to true for open shell
@@ -15,6 +15,26 @@
 
  ecmd_lda_mu_of_r = 0.d0
  call wall_time(wall0)
+ if (mu_of_r_potential.EQ."proj")then
+  do istate = 1, N_states
+    do ipoint = 1, n_points_final_grid
+    ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
+    mu     = mu_of_r_prov(ipoint,istate)
+    weight = final_weight_at_r_vector(ipoint)
+    rho = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
+    p2 = !TODO
+    rho_a = 0.5d0*(rho - dsqrt(-p2 + rho*rho))
+    rho_b = 0.5d0*(rho + dsqrt(-p2 + rho*rho))
+    ! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
+    call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
+    if(isnan(ec))then
+      print*,'ec is nan'
+      stop
+    endif
+    ecmd_lda_mu_of_r(istate) += weight * ec
+    enddo
+  enddo
+ else
   do istate = 1, N_states
     do ipoint = 1, n_points_final_grid
     ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
@@ -31,9 +51,10 @@
     ecmd_lda_mu_of_r(istate) += weight * ec
     enddo
   enddo
+ endif
  call wall_time(wall1)
  print*,'Time for ecmd_lda_mu_of_r :',wall1-wall0
- END_PROVIDER 
+END_PROVIDER
 
 
 BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]

scripts/qp_import_trexio.py

@@ -524,11 +524,11 @@ def write_ezfio(trexio_filename, filename):
     alpha = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(alpha_s) ][::-1]
     beta  = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(beta_s ) ][::-1]
     ezfio.set_determinants_bit_kind(8)
-    ezfio.set_determinants_n_int(1+mo_num//64)
+    ezfio.set_determinants_n_int(1+(mo_num-1)//64)
     ezfio.set_determinants_n_det(1)
     ezfio.set_determinants_n_states(1)
-    ezfio.set_determinants_psi_det(alpha+beta)
     ezfio.set_determinants_psi_coef([[1.0]])
+    ezfio.set_determinants_psi_det(alpha+beta)
     print("OK")
 
 

src/ao_one_e_ints/ao_ortho_canonical.irp.f

@@ -59,6 +59,28 @@
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! T^T . S . T
+ END_DOC
+ double precision, allocatable :: S(:,:)
+ allocate (S(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
+   ao_overlap,size(ao_overlap,1), 0.d0, &
+   S, size(S,1))
+
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
+   S, size(S,1), &
+   ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, &
+   ao_cart_to_sphe_overlap,size(ao_cart_to_sphe_overlap,1))
+
+ deallocate(S)
+
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ]
  implicit none
  BEGIN_DOC
@@ -82,7 +104,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ]
 
   integer :: i
   do i=1,ao_num
-    ao_cart_to_sphe_inv(:,i) = Rinv(:,i)  !/ ao_cart_to_sphe_normalization(i)
+    ao_cart_to_sphe_inv(:,i) = Rinv(:,i)
   enddo
 
 END_PROVIDER
@@ -132,7 +154,7 @@ END_PROVIDER
     enddo
 
     ao_ortho_canonical_num = ao_sphe_num
-    call ortho_canonical(ao_sphe_overlap, size(ao_sphe_overlap,1), &
+    call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
       ao_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff)
 
     call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_sphe_num, 1.d0, &

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -649,11 +649,20 @@ double precision function general_primitive_integral(dim,            &
 !  call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
   if (ior(n_pt_tmp,n_Iz) >= 0) then
     ! Bottleneck here
+    if (ic > ib) then
       do ib=0,n_pt_tmp
+          d1(ib:) = d1(ib:) + Iz_pol(:) * d_poly(ib)
+      enddo
+    else
       do ic = 0,n_Iz
-        d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib)
+        d1(ic:) = d1(ic:) + Iz_pol(ic) * d_poly(:)
-      enddo
       enddo
+    endif
+!    do ib=0,n_pt_tmp
+!      do ic = 0,n_Iz
+!        d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib)
+!      enddo
+!    enddo
 
     do n_pt_out = n_pt_tmp+n_Iz, 0, -1
       if (d1(n_pt_out) /= 0.d0) exit

src/scf_utils/roothaan_hall_scf.irp.f

@@ -222,7 +222,7 @@ END_DOC
   endif
 
 
-  ! Identify degenerate MOs and force them to be on the axes
+  ! Identify degenerate MOs and combine them to force them to be on the axes
   allocate(S(ao_num,ao_num))
   i=1
   do while (i<mo_num)
@@ -239,6 +239,10 @@ END_DOC
     endif
     i = j
   enddo
+  TOUCH mo_coef
+  call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), &
+       size(Fock_matrix_mo,2),mo_label,1,.true.)
+  call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
 
   if(do_mom)then
      call reorder_mo_max_overlap