diff --git a/etc/openmp.rc b/etc/openmp.rc index 7f71c9b8..a8ced5a0 100644 --- a/etc/openmp.rc +++ b/etc/openmp.rc @@ -1 +1 @@ -export OMP_NESTED=True +#export OMP_NESTED=True diff --git a/external/irpf90 b/external/irpf90 index 4ab1b175..43160c60 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6 +Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b diff --git a/plugins/local/basis_correction/weak_corr_func.irp.f b/plugins/local/basis_correction/weak_corr_func.irp.f index 6af5e49d..99ea504b 100644 --- a/plugins/local/basis_correction/weak_corr_func.irp.f +++ b/plugins/local/basis_correction/weak_corr_func.irp.f @@ -1,45 +1,73 @@ - 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 + double precision, external :: g0_UEG_mu_inf dospin = .true. ! JT dospin have to be set to true for open shell print*,'Providing ecmd_lda_mu_of_r ...' 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_DUMMY")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_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) + rho = rho_a + rho_b +! p2 = rho_a*rho_b +! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0 + p2 = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b) + rho_a = 0.5d0*(rho + dsqrt(rho*rho - 4.d0*p2)) + rho_b = 0.5d0*(rho - dsqrt(rho*rho - 4.d0*p2)) +! rho_a = 0.5d0*rho +! rho_b = 0.5d0*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) @@ -50,10 +78,10 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)] double precision :: weight integer :: ipoint,istate double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top - double precision :: g0_UEG_mu_inf + double precision, external :: 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 +97,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 a515efba..156c4344 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -210,6 +210,7 @@ def write_ezfio(trexio_filename, filename): nucl_index = trexio.read_basis_nucleus_index(trexio_file) exponent = [1.]*prim_num coefficient = [1.]*prim_num + prim_factor = [1.]*prim_num shell_index = [i for i in range(shell_num)] ao_shell = trexio.read_ao_shell(trexio_file) @@ -221,6 +222,9 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ]) ezfio.set_basis_prim_expo(exponent) ezfio.set_basis_prim_coef(coefficient) + ezfio.set_basis_prim_normalization_factor(prim_factor) + ezfio.set_basis_primitives_normalized(True) + ezfio.set_basis_ao_normalized(False) nucl_shell_num = [] prev = None @@ -271,12 +275,11 @@ def write_ezfio(trexio_filename, filename): if basis_type.lower() == "gaussian" and not cartesian: try: import trexio_tools - fd, tmp = tempfile.mkstemp() - os.close(fd) + tmp = "cartesian_"+trexio_filename retcode = subprocess.call(["trexio", "convert-to", "-t", "cartesian", "-o", tmp, trexio_filename]) trexio_file_cart = trexio.File(tmp,mode='r',back_end=trexio.TREXIO_AUTO) cartesian = trexio.read_ao_cartesian(trexio_file_cart) - os.unlink(tmp) + ezfio.set_trexio_trexio_file(tmp) except: pass @@ -284,7 +287,7 @@ def write_ezfio(trexio_filename, filename): ezfio.set_ao_basis_ao_num(ao_num) - if cartesian and basis_type.lower() == "gaussian" and shell_num > 0: + if cartesian and basis_type.lower() in ["gaussian", "numerical"] and shell_num > 0: ao_shell = trexio.read_ao_shell(trexio_file_cart) at = [ nucl_index[i]+1 for i in ao_shell ] ezfio.set_ao_basis_ao_nucl(at) @@ -521,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_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 440cc865..02eedf53 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -20,7 +20,35 @@ BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ] ao_shell(k) = i enddo enddo +END_PROVIDER +BEGIN_PROVIDER [ integer, ao_sphe_num ] + implicit none + BEGIN_DOC + ! Number of spherical AOs + END_DOC + integer :: n, i + ao_sphe_num=0 + do i=1,shell_num + n = shell_ang_mom(i) + ao_sphe_num += 2*n+1 + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer, ao_sphe_shell, (ao_sphe_num) ] + implicit none + BEGIN_DOC + ! Index of the shell to which the AO corresponds + END_DOC + integer :: i, j, k, n + k=0 + do i=1,shell_num + n = shell_ang_mom(i) + do j=-n,n + k = k+1 + ao_sphe_shell(k) = i + enddo + enddo END_PROVIDER BEGIN_PROVIDER [ integer, ao_first_of_shell, (shell_num) ] @@ -133,6 +161,16 @@ END_PROVIDER enddo enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_sphe_coef_normalization_factor, (ao_sphe_num) ] + implicit none + BEGIN_DOC + ! Normalization factor in spherical AO basis + END_DOC + + ao_sphe_coef_normalization_factor(:) = 1.d0 + END_PROVIDER BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num,ao_prim_num_max) ] diff --git a/src/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f index 336161f8..53a60413 100644 --- a/src/ao_basis/spherical_to_cartesian.irp.f +++ b/src/ao_basis/spherical_to_cartesian.irp.f @@ -4,7 +4,8 @@ ! First index is the index of the cartesian AO, obtained by ao_power_index ! Second index is the index of the spherical AO -BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_0, (1) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=0 @@ -12,10 +13,12 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ] cart_to_sphe_0 = 0.d0 cart_to_sphe_0 ( 1, 1) = 1.0d0 + cart_to_sphe_norm_0 (1) = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_1, (3) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=1 @@ -25,10 +28,14 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ] cart_to_sphe_1 ( 3, 1) = 1.0d0 cart_to_sphe_1 ( 1, 2) = 1.0d0 cart_to_sphe_1 ( 2, 3) = 1.0d0 + cart_to_sphe_norm_1 (1) = 1.d0 + cart_to_sphe_norm_1 (2) = 1.d0 + cart_to_sphe_norm_1 (3) = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_2, (6) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=2 @@ -43,10 +50,14 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ] cart_to_sphe_2 ( 1, 4) = 0.86602540378443864676d0 cart_to_sphe_2 ( 4, 4) = -0.86602540378443864676d0 cart_to_sphe_2 ( 2, 5) = 1.0d0 + + cart_to_sphe_norm_2 = (/ 1.0d0, 1.7320508075688772d0, 1.7320508075688772d0, 1.0d0, & + 1.7320508075688772d0, 1.0d0 /) END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_3, (10) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=3 @@ -69,10 +80,15 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ] cart_to_sphe_3 ( 4, 6) = -1.0606601717798212866d0 cart_to_sphe_3 ( 2, 7) = 1.0606601717798212866d0 cart_to_sphe_3 ( 7, 7) = -0.790569415042094833d0 + + cart_to_sphe_norm_3 = (/ 1.0d0, 2.23606797749979d0, 2.23606797749979d0, & + 2.23606797749979d0, 3.872983346207417d0, 2.23606797749979d0, 1.0d0, 2.23606797749979d0, & + 2.23606797749979d0, 1.d00 /) END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_4, (15) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=4 @@ -107,10 +123,18 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ] cart_to_sphe_4 (11, 8) = 0.73950997288745200532d0 cart_to_sphe_4 ( 2, 9) = 1.1180339887498948482d0 cart_to_sphe_4 ( 7, 9) = -1.1180339887498948482d0 + + cart_to_sphe_norm_4 = (/ 1.0d0, 2.6457513110645907d0, 2.6457513110645907d0, & + 3.4156502553198664d0, 5.916079783099616d0, 3.415650255319866d0, & + 2.6457513110645907d0, 5.916079783099616d0, 5.916079783099616d0, & + 2.6457513110645907d0, 1.0d0, 2.6457513110645907d0, 3.415650255319866d0, & + 2.6457513110645907d0, 1.d00 /) + END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_5, (21) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=5 @@ -163,10 +187,18 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ] cart_to_sphe_5 ( 2,11) = 1.169267933366856683d0 cart_to_sphe_5 ( 7,11) = -1.5309310892394863114d0 cart_to_sphe_5 (16,11) = 0.7015607600201140098d0 + + cart_to_sphe_norm_5 = (/ 1.0d0, 3.0d0, 3.0d0, 4.58257569495584d0, & + 7.937253933193773d0, 4.58257569495584d0, 4.58257569495584d0, & + 10.246950765959598d0, 10.246950765959598d0, 4.582575694955841d0, 3.0d0, & + 7.937253933193773d0, 10.246950765959598d0, 7.937253933193773d0, 3.0d0, 1.0d0, & + 3.0d0, 4.58257569495584d0, 4.582575694955841d0, 3.0d0, 1.d00 /) + END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_6, (28) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=6 @@ -243,10 +275,22 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ] cart_to_sphe_6 ( 2,13) = 1.2151388809514737933d0 cart_to_sphe_6 ( 7,13) = -1.9764235376052370825d0 cart_to_sphe_6 (16,13) = 1.2151388809514737933d0 + + cart_to_sphe_norm_6 = (/ 1.0d0, 3.3166247903554003d0, 3.3166247903554003d0, & + 5.744562646538029d0, 9.949874371066201d0, 5.744562646538029d0, & + 6.797058187186571d0, 15.198684153570666d0, 15.198684153570664d0, & + 6.797058187186572d0, 5.744562646538029d0, 15.198684153570666d0, & + 19.621416870348583d0, 15.198684153570666d0, 5.744562646538029d0, & + 3.3166247903554003d0, 9.949874371066201d0, 15.198684153570664d0, & + 15.198684153570666d0, 9.9498743710662d0, 3.3166247903554003d0, 1.0d0, & + 3.3166247903554003d0, 5.744562646538029d0, 6.797058187186572d0, & + 5.744562646538029d0, 3.3166247903554003d0, 1.d00 /) + END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_7, (36) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=7 @@ -355,10 +399,25 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ] cart_to_sphe_7 ( 7,15) = -2.4456993503903949804d0 cart_to_sphe_7 (16,15) = 1.96875d0 cart_to_sphe_7 (29,15) = -0.64725984928774934788d0 + + cart_to_sphe_norm_7 = (/ 1.0d0, 3.6055512754639896d0, 3.605551275463989d0, & + 6.904105059069327d0, 11.958260743101398d0, 6.904105059069326d0, & + 9.26282894152753d0, 20.712315177207984d0, 20.71231517720798d0, & + 9.26282894152753d0, 9.26282894152753d0, 24.507141816213494d0, & + 31.63858403911275d0, 24.507141816213494d0, 9.262828941527529d0, & + 6.904105059069327d0, 20.712315177207984d0, 31.63858403911275d0, & + 31.63858403911275d0, 20.71231517720798d0, 6.904105059069327d0, & + 3.6055512754639896d0, 11.958260743101398d0, 20.71231517720798d0, & + 24.507141816213494d0, 20.71231517720798d0, 11.958260743101398d0, & + 3.6055512754639896d0, 1.0d0, 3.605551275463989d0, 6.904105059069326d0, & + 9.26282894152753d0, 9.262828941527529d0, 6.904105059069327d0, & + 3.6055512754639896d0, 1.d00 /) + END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_8, (45) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=8 @@ -506,10 +565,28 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ] cart_to_sphe_8 ( 7,17) = -2.9348392204684739765d0 cart_to_sphe_8 (16,17) = 2.9348392204684739765d0 cart_to_sphe_8 (29,17) = -1.2945196985754986958d0 + + cart_to_sphe_norm_8 = (/ 1.0d0, 3.872983346207417d0, 3.872983346207417d0, & + 8.062257748298551d0, 13.964240043768942d0, 8.06225774829855d0, & + 11.958260743101398d0, 26.739483914241877d0, 26.739483914241877d0, & + 11.958260743101398d0, 13.55939315961975d0, 35.874782229304195d0, & + 46.31414470763765d0, 35.874782229304195d0, 13.55939315961975d0, & + 11.958260743101398d0, 35.874782229304195d0, 54.79963503528103d0, & + 54.79963503528103d0, 35.874782229304195d0, 11.958260743101398d0, & + 8.062257748298551d0, 26.739483914241877d0, 46.31414470763765d0, & + 54.79963503528103d0, 46.314144707637645d0, 26.739483914241877d0, & + 8.06225774829855d0, 3.872983346207417d0, 13.964240043768942d0, & + 26.739483914241877d0, 35.874782229304195d0, 35.874782229304195d0, & + 26.739483914241877d0, 13.96424004376894d0, 3.8729833462074166d0, 1.0d0, & + 3.872983346207417d0, 8.06225774829855d0, 11.958260743101398d0, & + 13.55939315961975d0, 11.958260743101398d0, 8.06225774829855d0, & + 3.8729833462074166d0, 1.d0 /) + END_PROVIDER -BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ] + BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ] +&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_9, (55) ] implicit none BEGIN_DOC ! Spherical -> Cartesian Transformation matrix for l=9 @@ -703,5 +780,28 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ] cart_to_sphe_9 (16,19) = 4.1179360680974030877d0 cart_to_sphe_9 (29,19) = -2.3781845426185916576d0 cart_to_sphe_9 (46,19) = 0.60904939217552380708d0 + + cart_to_sphe_norm_9 = (/ 1.0d0, 4.1231056256176615d0, 4.1231056256176615d0, & + 9.219544457292889d0, 15.968719422671313d0, 9.219544457292889d0, & + 14.86606874731851d0, 33.24154027718933d0, 33.24154027718933d0, & + 14.866068747318508d0, 18.635603405463275d0, 49.30517214248421d0, & + 63.652703529910404d0, 49.30517214248421d0, 18.635603405463275d0, & + 18.635603405463275d0, 55.90681021638982d0, 85.39906322671229d0, & + 85.39906322671229d0, 55.90681021638983d0, 18.635603405463275d0, & + 14.86606874731851d0, 49.30517214248421d0, 85.39906322671229d0, & + 101.04553429023969d0, 85.3990632267123d0, 49.30517214248421d0, & + 14.866068747318508d0, 9.219544457292889d0, 33.24154027718933d0, & + 63.652703529910404d0, 85.39906322671229d0, 85.3990632267123d0, & + 63.65270352991039d0, 33.24154027718933d0, 9.219544457292887d0, & + 4.1231056256176615d0, 15.968719422671313d0, 33.24154027718933d0, & + 49.30517214248421d0, 55.90681021638983d0, 49.30517214248421d0, & + 33.24154027718933d0, 15.968719422671313d0, 4.1231056256176615d0, 1.0d0, & + 4.1231056256176615d0, 9.219544457292889d0, 14.866068747318508d0, & + 18.635603405463275d0, 18.635603405463275d0, 14.866068747318508d0, & + 9.219544457292887d0, 4.1231056256176615d0, 1.d0 /) + END_PROVIDER + + + diff --git a/src/ao_one_e_ints/ao_one_e_ints.irp.f b/src/ao_one_e_ints/ao_one_e_ints.irp.f index 2c6b8e7e..8b0edfbc 100644 --- a/src/ao_one_e_ints/ao_one_e_ints.irp.f +++ b/src/ao_one_e_ints/ao_one_e_ints.irp.f @@ -45,3 +45,19 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)] END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals,(ao_sphe_num,ao_sphe_num)] +&BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals_diag,(ao_sphe_num)] + implicit none + integer :: i,j,n,l + BEGIN_DOC + ! One-electron Hamiltonian in the spherical |AO| basis. + END_DOC + + ao_sphe_one_e_integrals = ao_sphe_integrals_n_e + ao_sphe_kinetic_integrals + + do j = 1, ao_num + ao_sphe_one_e_integrals_diag(j) = ao_sphe_one_e_integrals(j,j) + enddo + +END_PROVIDER + 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 eff7e7be..523e49f7 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -1,35 +1,42 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_num)] -&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ] +&BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_num)] implicit none BEGIN_DOC ! Coefficients to go from cartesian to spherical coordinates in the current ! basis set +! +! S_cart^-1 END_DOC integer :: i integer, external :: ao_power_index integer :: ibegin,j,k - integer :: prev + integer :: prev, ao_sphe_count prev = 0 ao_cart_to_sphe_coef(:,:) = 0.d0 + ao_cart_to_sphe_normalization(:) = 1.d0 ! Assume order provided by ao_power_index i = 1 - ao_cart_to_sphe_num = 0 + ao_sphe_count = 0 do while (i <= ao_num) select case ( ao_l(i) ) case (0) - ao_cart_to_sphe_num += 1 - ao_cart_to_sphe_coef(i,ao_cart_to_sphe_num) = 1.d0 + ao_sphe_count += 1 + ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 + ao_cart_to_sphe_normalization(i) = 1.d0 i += 1 BEGIN_TEMPLATE case ($SHELL) if (ao_power(i,1) == $SHELL) then do k=1,size(cart_to_sphe_$SHELL,2) do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_coef(i+j-1,ao_cart_to_sphe_num+k) = cart_to_sphe_$SHELL(j,k) + ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) enddo enddo + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j) + enddo i += size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_num += size(cart_to_sphe_$SHELL,2) + ao_sphe_count += size(cart_to_sphe_$SHELL,2) endif SUBST [ SHELL ] 1;; @@ -47,22 +54,25 @@ end select enddo + if (ao_sphe_count /= ao_sphe_num) then + call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num") + endif END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ] +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ] implicit none BEGIN_DOC - ! |AO| overlap matrix in the spherical basis set + ! T^T . S . T END_DOC double precision, allocatable :: S(:,:) - allocate (S(ao_cart_to_sphe_num,ao_num)) + allocate (S(ao_sphe_num,ao_num)) - call dgemm('T','N',ao_cart_to_sphe_num,ao_num,ao_num, 1.d0, & + 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_cart_to_sphe_num,ao_cart_to_sphe_num,ao_num, 1.d0, & + 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)) @@ -71,15 +81,33 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_num) ] +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ] implicit none BEGIN_DOC ! Inverse of :c:data:`ao_cart_to_sphe_coef` END_DOC - call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),& - ao_num,ao_cart_to_sphe_num, & - ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1), lin_dep_cutoff) + ! Normalize + integer :: m,k + double precision, allocatable :: S(:,:), R(:,:), Rinv(:,:), Sinv(:,:) + + k = size(ao_cart_to_sphe_coef,1) + m = size(ao_cart_to_sphe_coef,2) + + allocate(S(k,k), R(k,m), Rinv(m,k), Sinv(k,k)) + + R(:,:) = ao_cart_to_sphe_coef(:,:) + + call dgemm('N','T', m, m, k, 1.d0, R, k, R, k, 0.d0, S, m) + call get_pseudo_inverse(S, k, k, m, Sinv, k, 1.d-12) + call dgemm('T','N', m, m, k, 1.d0, R, k, Sinv, k, 0.d0, Rinv, m) + + + integer :: i + do i=1,ao_num + ao_cart_to_sphe_inv(:,i) = Rinv(:,i) + enddo + END_PROVIDER @@ -120,17 +148,17 @@ END_PROVIDER double precision, allocatable :: S(:,:) - allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num)) + allocate(S(ao_sphe_num,ao_sphe_num)) S = 0.d0 - do i=1,ao_cart_to_sphe_num + do i=1,ao_sphe_num S(i,i) = 1.d0 enddo - ao_ortho_canonical_num = ao_cart_to_sphe_num + ao_ortho_canonical_num = ao_sphe_num call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & - ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff) + 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_cart_to_sphe_num, 1.d0, & + call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_sphe_num, 1.d0, & ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & S, size(S,1), & 0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1)) @@ -167,3 +195,4 @@ BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonica enddo enddo END_PROVIDER + diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index e5341f6a..17eb7cbc 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -308,3 +308,26 @@ BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ] END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_sphe_overlap, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! |AO| overlap matrix in the spherical basis set + END_DOC + double precision, allocatable :: tmp(:,:) + allocate (tmp(ao_sphe_num,ao_num)) + + call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & + ao_overlap,size(ao_overlap,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & + tmp, size(tmp,1), & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & + ao_sphe_overlap,size(ao_sphe_overlap,1)) + + deallocate(tmp) + +END_PROVIDER + diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f index 49eb56ad..1e3e684d 100644 --- a/src/ao_one_e_ints/kin_ao_ints.irp.f +++ b/src/ao_one_e_ints/kin_ao_ints.irp.f @@ -190,3 +190,25 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)] endif END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_sphe_kinetic_integrals, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! |AO| kinetic inntegrals matrix in the spherical basis set + END_DOC + double precision, allocatable :: tmp(:,:) + allocate (tmp(ao_sphe_num,ao_num)) + + call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & + ao_kinetic_integrals,size(ao_kinetic_integrals,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & + tmp, size(tmp,1), & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & + ao_sphe_kinetic_integrals,size(ao_sphe_kinetic_integrals,1)) + + deallocate(tmp) + +END_PROVIDER diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 776b5ec0..b4b4aa02 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -609,3 +609,25 @@ double precision function V_r(n,alpha) end + +BEGIN_PROVIDER [ double precision, ao_sphe_integrals_n_e, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! |AO| VneVne inntegrals matrix in the spherical basis set + END_DOC + double precision, allocatable :: tmp(:,:) + allocate (tmp(ao_sphe_num,ao_num)) + + call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & + ao_integrals_n_e,size(ao_integrals_n_e,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & + tmp, size(tmp,1), & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & + ao_sphe_integrals_n_e,size(ao_sphe_integrals_n_e,1)) + + deallocate(tmp) + +END_PROVIDER diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index e75ca056..851f26d8 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -296,3 +296,67 @@ END_PROVIDER enddo END_PROVIDER +BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_local, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! |AO| pseudo_integrals_local matrix in the spherical basis set + END_DOC + double precision, allocatable :: tmp(:,:) + allocate (tmp(ao_sphe_num,ao_num)) + + call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & + ao_pseudo_integrals_local,size(ao_pseudo_integrals_local,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & + tmp, size(tmp,1), & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & + ao_sphe_pseudo_integrals_local,size(ao_sphe_pseudo_integrals_local,1)) + + deallocate(tmp) + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_non_local, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! |AO| pseudo_integrals_non_local matrix in the spherical basis set + END_DOC + double precision, allocatable :: tmp(:,:) + allocate (tmp(ao_sphe_num,ao_num)) + + call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & + ao_pseudo_integrals_non_local,size(ao_pseudo_integrals_non_local,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & + tmp, size(tmp,1), & + ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & + ao_sphe_pseudo_integrals_non_local,size(ao_sphe_pseudo_integrals_non_local,1)) + + deallocate(tmp) + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals, (ao_sphe_num,ao_sphe_num)] + implicit none + BEGIN_DOC + ! Pseudo-potential integrals in the |AO| basis set. + END_DOC + + ao_sphe_pseudo_integrals = 0.d0 + if (do_pseudo) then + if (pseudo_klocmax > 0) then + ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_local + endif + if (pseudo_kmax > 0) then + ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_non_local + endif + endif + +END_PROVIDER + 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/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 1ed859ee..4da8ac9b 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -346,3 +346,37 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) deallocate(T) end + +BEGIN_PROVIDER [ double precision, mo_sphe_coef, (ao_sphe_num, mo_num) ] + implicit none + BEGIN_DOC + ! MO coefficients in the basis of spherical harmonics AOs. + END_DOC + double precision, allocatable, dimension (:,:) :: C0, S, tmp + allocate(C0(ao_sphe_num,mo_num), S(mo_num,mo_num), tmp(mo_num,ao_sphe_num)) + + call dgemm('T','N',ao_sphe_num,mo_num,ao_num, 1.d0, & + ao_cart_to_sphe_coef,ao_num, & + mo_coef,size(mo_coef,1), 0.d0, & + C0, size(C0,1)) + + ! C0^T S S^0 + call dgemm('T','N',mo_num,ao_sphe_num,ao_sphe_num, 1.d0, & + C0,size(C0,1), & + ao_sphe_overlap,size(ao_sphe_overlap,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','N',mo_num,mo_num,ao_sphe_num, 1.d0, & + tmp, size(tmp,1), & + C0,size(C0,1), 0.d0, & + S, size(S,1)) + + integer :: m + m = ao_sphe_num + mo_sphe_coef = C0 + call ortho_lowdin(S,size(S,1),mo_num,mo_sphe_coef,ao_sphe_num,m,1.d-10) + + + deallocate (tmp, S, C0) +END_PROVIDER + diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 4ba0964a..7a73121b 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= 20250211) then + rc = trexio_write_ao_normalization(f(1), ao_coef_normalization_factor) + else + + allocate(factor(ao_num)) + do i=1,ao_num + l = ao_first_of_shell(ao_shell(i)) + factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0)) + enddo + rc = trexio_write_ao_normalization(f(1), factor) + deallocate(factor) + endif + + call trexio_assert(rc, TREXIO_SUCCESS) - rc = trexio_write_ao_shell(f(1), ao_shell) - call trexio_assert(rc, TREXIO_SUCCESS) - if (ezfio_convention >= 20250211) then - rc = trexio_write_ao_normalization(f(1), ao_coef_normalization_factor) else - integer :: pow0(3), powA(3), nz - double precision :: normA, norm0, C_A(3), overlap_x, overlap_z, overlap_y, c - nz=100 + rc = trexio_write_ao_cartesian(f(1), 0) + call trexio_assert(rc, TREXIO_SUCCESS) - C_A(1) = 0.d0 - C_A(2) = 0.d0 - C_A(3) = 0.d0 + rc = trexio_write_ao_num(f(1), ao_sphe_num) + call trexio_assert(rc, TREXIO_SUCCESS) + + rc = trexio_write_ao_shell(f(1), ao_sphe_shell) + call trexio_assert(rc, TREXIO_SUCCESS) + + rc = trexio_write_ao_normalization(f(1), ao_sphe_coef_normalization_factor) + call trexio_assert(rc, TREXIO_SUCCESS) - allocate(factor(ao_num)) - do i=1,ao_num - l = ao_first_of_shell(ao_shell(i)) - factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0)) - enddo - rc = trexio_write_ao_normalization(f(1), factor) - deallocate(factor) endif - call trexio_assert(rc, TREXIO_SUCCESS) endif @@ -341,23 +356,45 @@ subroutine export_trexio(update,full_path) ! ------------------ if (export_ao_one_e_ints) then - print *, 'AO one-e integrals' - rc = trexio_write_ao_1e_int_overlap(f(1),ao_overlap) - call trexio_assert(rc, TREXIO_SUCCESS) + double precision, allocatable :: tmp_ao(:,:) + if (export_cartesian) then + print *, 'AO one-e integrals (cartesian)' - rc = trexio_write_ao_1e_int_kinetic(f(1),ao_kinetic_integrals) - call trexio_assert(rc, TREXIO_SUCCESS) - - rc = trexio_write_ao_1e_int_potential_n_e(f(1),ao_integrals_n_e) - call trexio_assert(rc, TREXIO_SUCCESS) - - if (do_pseudo) then - rc = trexio_write_ao_1e_int_ecp(f(1), ao_pseudo_integrals_local + ao_pseudo_integrals_non_local) + rc = trexio_write_ao_1e_int_overlap(f(1),ao_overlap) call trexio_assert(rc, TREXIO_SUCCESS) - endif - rc = trexio_write_ao_1e_int_core_hamiltonian(f(1),ao_one_e_integrals) + rc = trexio_write_ao_1e_int_kinetic(f(1),ao_kinetic_integrals) + call trexio_assert(rc, TREXIO_SUCCESS) + + rc = trexio_write_ao_1e_int_potential_n_e(f(1),ao_integrals_n_e) + call trexio_assert(rc, TREXIO_SUCCESS) + + if (do_pseudo) then + rc = trexio_write_ao_1e_int_ecp(f(1), ao_pseudo_integrals) + call trexio_assert(rc, TREXIO_SUCCESS) + endif + + rc = trexio_write_ao_1e_int_core_hamiltonian(f(1),ao_one_e_integrals) + else + print *, 'AO one-e integrals (spherical)' + + rc = trexio_write_ao_1e_int_overlap(f(1),ao_sphe_overlap) + call trexio_assert(rc, TREXIO_SUCCESS) + + rc = trexio_write_ao_1e_int_kinetic(f(1),ao_sphe_kinetic_integrals) + call trexio_assert(rc, TREXIO_SUCCESS) + + rc = trexio_write_ao_1e_int_potential_n_e(f(1),ao_sphe_integrals_n_e) + call trexio_assert(rc, TREXIO_SUCCESS) + + if (do_pseudo) then + rc = trexio_write_ao_1e_int_ecp(f(1), ao_sphe_pseudo_integrals) + call trexio_assert(rc, TREXIO_SUCCESS) + endif + + rc = trexio_write_ao_1e_int_core_hamiltonian(f(1),ao_sphe_one_e_integrals) + endif call trexio_assert(rc, TREXIO_SUCCESS) end if @@ -465,8 +502,13 @@ subroutine export_trexio(update,full_path) call trexio_assert(rc, TREXIO_SUCCESS) enddo - rc = trexio_write_mo_coefficient(f(1), mo_coef) - call trexio_assert(rc, TREXIO_SUCCESS) + if (export_cartesian) then + rc = trexio_write_mo_coefficient(f(1), mo_coef) + call trexio_assert(rc, TREXIO_SUCCESS) + else + rc = trexio_write_mo_coefficient(f(1), mo_sphe_coef) + call trexio_assert(rc, TREXIO_SUCCESS) + endif if ( (trim(mo_label) == 'Canonical').and. & (export_mo_two_e_ints_cholesky.or.export_mo_two_e_ints) ) then diff --git a/src/trexio/import_trexio_integrals.irp.f b/src/trexio/import_trexio_integrals.irp.f index d85a0283..0fda2ff0 100644 --- a/src/trexio/import_trexio_integrals.irp.f +++ b/src/trexio/import_trexio_integrals.irp.f @@ -119,6 +119,28 @@ subroutine run(f) call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e('Read') endif + ! Some codes only provide ao_1e_int_core_hamiltonian rather than + ! kinetic and nuclear potentials separately, so we need to work + ! around that. This is needed for non-GTO basis sets since some QP + ! functions will try to calculate these matrices from the nonexisting + ! GTO basis if they are not set. + if (trexio_has_ao_1e_int_core_hamiltonian(f) == TREXIO_SUCCESS .and. & + trexio_has_ao_1e_int_potential_n_e(f) /= TREXIO_SUCCESS .and. & + trexio_has_ao_1e_int_kinetic(f) /= TREXIO_SUCCESS) then + rc = trexio_read_ao_1e_int_core_hamiltonian(f, A) + if (rc /= TREXIO_SUCCESS) then + print *, irp_here + print *, 'Error reading AO core Hamiltonian.' + call trexio_assert(rc, TREXIO_SUCCESS) + stop -1 + endif + call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A) + call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e('Read') + A=0.d0 + call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A) + call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic('Read') + endif + deallocate(A,B) ! AO 2e integrals diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 4e7ca87d..629a998b 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1392,15 +1392,6 @@ subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff) call dgemm('T', 'T', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1)) -! C = 0.d0 -! do i=1,m -! do j=1,n -! do k=1,n_svd -! C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) -! enddo -! enddo -! enddo - deallocate(U,D,Vt,work,A_tmp) end