diff --git a/configure b/configure index 43ca9f6d..7b5fa0c9 100755 --- a/configure +++ b/configure @@ -9,7 +9,7 @@ echo "QP_ROOT="$QP_ROOT unset CC unset CCXX -TREXIO_VERSION=2.4.2 +TREXIO_VERSION=2.5.0 # Force GCC instead of ICC for dependencies export CC=gcc diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index 4373d90e..d75d0074 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -160,7 +160,6 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_shell_ang_mom(ang_mom) 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) nucl_shell_num = [] prev = None @@ -194,6 +193,10 @@ def write_ezfio(trexio_filename, filename): shell_factor = trexio.read_basis_shell_factor(trexio_file) prim_factor = trexio.read_basis_prim_factor(trexio_file) + for i,p in enumerate(prim_factor): + coefficient[i] *= prim_factor[i] + ezfio.set_ao_basis_primitives_normalized(False) + ezfio.set_basis_prim_coef(coefficient) elif basis_type.lower() == "numerical": @@ -245,13 +248,12 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_nucleus_shell_num(nucl_shell_num) shell_factor = trexio.read_basis_shell_factor(trexio_file) - prim_factor = [1.]*prim_num else: raise TypeError print(basis_type) except: - raise + basis_type = "None" print("None") ezfio.set_ao_basis_ao_cartesian(True) @@ -262,9 +264,6 @@ def write_ezfio(trexio_filename, filename): except: cartesian = True - ao_num = trexio.read_ao_num(trexio_file) - ezfio.set_ao_basis_ao_num(ao_num) - trexio_file_cart = trexio_file if basis_type.lower() == "gaussian" and not cartesian: try: @@ -278,8 +277,12 @@ def write_ezfio(trexio_filename, filename): except: pass + ao_num = trexio.read_ao_num(trexio_file_cart) + ezfio.set_ao_basis_ao_num(ao_num) + + if cartesian and basis_type.lower() == "gaussian" and shell_num > 0: - ao_shell = trexio.read_ao_shell(trexio_file) + 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) @@ -314,6 +317,7 @@ def write_ezfio(trexio_filename, filename): exponent.append(expo[i]) num_prim.append(num_prim0[i]) + print (len(coefficient), ao_num) assert (len(coefficient) == ao_num) ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) ezfio.set_ao_basis_ao_prim_num(num_prim) @@ -384,6 +388,14 @@ def write_ezfio(trexio_filename, filename): # Read coefs from temporary cartesian file created in the AO section MoMatrix = trexio.read_mo_coefficient(trexio_file_cart) + + # Renormalize MO coefs if needed + if trexio.has_ao_normalization(trexio_file_cart): + ezfio.set_ao_basis_ao_normalized(False) + norm = trexio.read_ao_normalization(trexio_file_cart) +# for j in range(mo_num): +# for i,f in enumerate(norm): +# MoMatrix[i,j] *= f ezfio.set_mo_basis_mo_coef(MoMatrix) mo_occ = [ 0. for i in range(mo_num) ] @@ -486,10 +498,10 @@ def write_ezfio(trexio_filename, filename): if trexio.has_mo_spin(trexio_file): spin = trexio.read_mo_spin(trexio_file) if max(spin) == 1: - alpha = [ i for i in range(len(spin)) if spin[i] == 0 ] - alpha = [ alpha[i] for i in range(num_alpha) ] - beta = [ i for i in range(len(spin)) if spin[i] == 1 ] - beta = [ beta[i] for i in range(num_beta) ] + tmp = [ i for i in range(len(spin)) if spin[i] == 0 ] + alpha = [ tmp[i] for i in range(num_alpha) ] + tmp = [ i for i in range(len(spin)) if spin[i] == 1 ] + beta = [ tmp[i] for i in range(num_beta) ] warnings.append("UHF orbitals orbitals read", end=' ') alpha_s = ['0']*mo_num beta_s = ['0']*mo_num diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 09604487..34853398 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -55,9 +55,6 @@ END_PROVIDER do i=1,ao_num -! powA(1) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) -! powA(2) = 0 -! powA(3) = 0 powA(1) = ao_power(i,1) powA(2) = ao_power(i,2) powA(3) = ao_power(i,3) @@ -76,16 +73,19 @@ END_PROVIDER endif ! Normalization of the contracted basis functions - if (ao_normalized) then - norm = 0.d0 - do j=1,ao_prim_num(i) - do k=1,ao_prim_num(i) - call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k) - enddo + norm = 0.d0 + do j=1,ao_prim_num(i) + do k=1,ao_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k) + enddo + enddo + ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm) + + if (.not.ao_normalized) then + do j=1,ao_prim_num(i) + ao_coef_normalized(i,j) = ao_coef_normalized(i,j) * ao_coef_normalization_factor(i) enddo - ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm) - else ao_coef_normalization_factor(i) = 1.d0 endif enddo @@ -339,3 +339,22 @@ BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] ao_l_char_space(i) = give_ao_character_space enddo END_PROVIDER + +! --- + +BEGIN_PROVIDER [ logical, use_pw ] + + implicit none + + logical :: exist + + use_pw = .false. + + call ezfio_has_ao_basis_ao_expo_pw(exist) + if(exist) then + PROVIDE ao_expo_pw_ord_transp + if(maxval(dabs(ao_expo_pw_ord_transp(4,:,:))) .gt. 1d-15) use_pw = .true. + endif + +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 8ab9fe5f..2c6b8e7e 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,4 +45,3 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)] END_PROVIDER - diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index efafd504..69b18900 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -314,7 +314,7 @@ END_PROVIDER do j=1,nq - if ( (Qmax < Dmin).or.(N+j*1_8 > ndim8) ) exit + if ( (Qmax <= Dmin).or.(N+j*1_8 > ndim8) ) exit ! i. rank = N+j 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 3734d4a0..e4bd9d1d 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -40,8 +40,11 @@ double precision function ao_two_e_integral(i, j, k, l) double precision, external :: ao_two_e_integral_erf double precision, external :: ao_two_e_integral_cgtos double precision, external :: ao_two_e_integral_schwartz_accel + double precision, external :: ao_two_e_integral_general + double precision, external :: general_primitive_integral logical, external :: do_schwartz_accel + double precision :: coef1, coef2, coef3, coef4 if(use_cgtos) then @@ -58,83 +61,44 @@ double precision function ao_two_e_integral(i, j, k, l) else - dim1 = n_pt_max_integrals - num_i = ao_nucl(i) num_j = ao_nucl(j) num_k = ao_nucl(k) num_l = ao_nucl(l) ao_two_e_integral = 0.d0 - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo + if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then - double precision :: coef1, coef2, coef3, coef4 - double precision :: p_inv,q_inv - double precision :: general_primitive_integral + ao_two_e_integral = ao_two_e_integral_general(i,j,k,l,general_primitive_integral) - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = general_primitive_integral(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p + else - else + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + enddo + double precision :: ERI - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p + do p = 1, ao_prim_num(i) + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) + do r = 1, ao_prim_num(k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + integral = ERI( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& + I_power(1),J_power(1),K_power(1),L_power(1), & + I_power(2),J_power(2),K_power(2),L_power(2), & + I_power(3),J_power(3),K_power(3),L_power(3)) + ao_two_e_integral = ao_two_e_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p endif @@ -144,6 +108,76 @@ end ! --- +double precision function ao_two_e_integral_general(i, j, k, l, op) + + BEGIN_DOC + ! integral of the AO basis or (ij|kl) + ! i(r1) j(r1) 1/r12 k(r2) l(r2) + END_DOC + + implicit none + include 'utils/constants.include.F' + + integer, intent(in) :: i, j, k, l + double precision, external :: op ! Operator function + + integer :: p, q, r, s + integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) + integer :: iorder_p(3), iorder_q(3) + double precision :: I_center(3), J_center(3), K_center(3), L_center(3) + double precision :: integral + double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp + double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq + + dim1 = n_pt_max_integrals + + num_i = ao_nucl(i) + num_j = ao_nucl(j) + num_k = ao_nucl(k) + num_l = ao_nucl(l) + ao_two_e_integral_general = 0.d0 + + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + I_center(p) = nucl_coord(num_i,p) + J_center(p) = nucl_coord(num_j,p) + K_center(p) = nucl_coord(num_k,p) + L_center(p) = nucl_coord(num_l,p) + enddo + + double precision :: coef1, coef2, coef3, coef4 + double precision :: p_inv,q_inv + + do p = 1, ao_prim_num(i) + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) + call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & + I_power,J_power,I_center,J_center,dim1) + p_inv = 1.d0/pp + do r = 1, ao_prim_num(k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & + K_power,L_power,K_center,L_center,dim1) + q_inv = 1.d0/qq + integral = op(dim1, & + P_new,P_center,fact_p,pp,p_inv,iorder_p, & + Q_new,Q_center,fact_q,qq,q_inv,iorder_q) + ao_two_e_integral_general = ao_two_e_integral_general + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + +end + double precision function ao_two_e_integral_schwartz_accel(i,j,k,l) implicit none BEGIN_DOC @@ -512,7 +546,7 @@ double precision function general_primitive_integral(dim, & double precision :: a,b,c,d,e,f,accu,pq,const double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2 integer :: n_pt_tmp,n_pt_out, iorder - double precision :: d1(0:max_dim),d_poly(0:max_dim),rint,d1_screened(0:max_dim) + double precision :: d1(0:max_dim),d_poly(0:max_dim),d1_screened(0:max_dim) general_primitive_integral = 0.d0 diff --git a/src/mo_one_e_ints/mo_one_e_ints.irp.f b/src/mo_one_e_ints/mo_one_e_ints.irp.f index fba55beb..19cd2159 100644 --- a/src/mo_one_e_ints/mo_one_e_ints.irp.f +++ b/src/mo_one_e_ints/mo_one_e_ints.irp.f @@ -21,3 +21,14 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)] call nullify_small_elements(mo_num,mo_num,mo_one_e_integrals,size(mo_one_e_integrals,1),1.d-15) END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_one_e_integrals_from_mo, (ao_num, ao_num)] + implicit none + BEGIN_DOC +! Integrals of the one e hamiltonian obtained from the integrals on the MO basis +! +! WARNING : this is equal to ao_one_e_integrals only if the AO and MO basis have the same number of functions + END_DOC + call mo_to_ao(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_num) +END_PROVIDER