diff --git a/plugins/local/basis_correction/weak_corr_func.irp.f b/plugins/local/basis_correction/weak_corr_func.irp.f index 6af5e49d..a31ff3af 100644 --- a/plugins/local/basis_correction/weak_corr_func.irp.f +++ b/plugins/local/basis_correction/weak_corr_func.irp.f @@ -1,13 +1,13 @@ - BEGIN_PROVIDER [double precision, ecmd_lda_mu_of_r, (N_states)] + BEGIN_PROVIDER [double precision, ecmd_lda_mu_of_r, (N_states)] BEGIN_DOC -! ecmd_lda_mu_of_r = multi-determinantal Ecmd within the LDA approximation with mu(r) , +! ecmd_lda_mu_of_r = multi-determinantal Ecmd within the LDA approximation with mu(r) , ! ! see equation 40 in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714 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,31 +15,52 @@ ecmd_lda_mu_of_r = 0.d0 call wall_time(wall0) - 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_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) - rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate) - ! 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 + 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 - 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) + mu = mu_of_r_prov(ipoint,istate) + weight = final_weight_at_r_vector(ipoint) + rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate) + ! 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 + 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)] BEGIN_DOC ! ecmd_pbe_ueg_mu_of_r = multi-determinantal Ecmd within the PBE-UEG approximation with mu(r) , - ! + ! ! see Eqs. 13-14b in Phys.Chem.Lett.2019, 10, 2931 2937; https://pubs.acs.org/doi/10.1021/acs.jpclett.9b01176 ! ! Based on the PBE-on-top functional (see Eqs. 26, 27 of J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.5082638) @@ -53,7 +74,7 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)] double precision :: g0_UEG_mu_inf ecmd_pbe_ueg_mu_of_r = 0.d0 - + print*,'Providing ecmd_pbe_ueg_mu_of_r ...' call wall_time(wall0) do istate = 1, N_states @@ -69,7 +90,7 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)] grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate) ! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0 - on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b) + on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b) ! The form of interpolated (mu=0 ---> mu=infinity) functional originally introduced in JCP, 150, 084103 1-10 (2019) call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_PBE) diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index 00c1a40c..156c4344 100755 --- a/scripts/qp_import_trexio.py +++ b/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") diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index b19be1e2..bdfa2096 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/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, & diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 1cb7617e..a04eae69 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/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 - do ib=0,n_pt_tmp - do ic = 0,n_Iz - d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib) + if (ic > ib) then + do ib=0,n_pt_tmp + d1(ib:) = d1(ib:) + Iz_pol(:) * d_poly(ib) enddo - enddo + else + do ic = 0,n_Iz + d1(ic:) = d1(ic:) + Iz_pol(ic) * d_poly(:) + 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 diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 4ba0964a..66bedb07 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/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