From 14d59648af826ab7ec2c938cd3e52efcaf0a823a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 24 Apr 2015 21:45:18 +0200 Subject: [PATCH] Changed ao_coef to ao_coef_normalized_ordered --- ocaml/.merlin | 4 + src/AOs/ASSUMPTIONS.rst | 7 - src/AOs/README.rst | 90 +++---- src/AOs/ao_overlap.irp.f | 20 +- src/AOs/aos.irp.f | 299 +++++++++++---------- src/Bielec_integrals/ao_bi_integrals.irp.f | 58 ++-- src/Dets/H_apply.irp.f | 35 +++ src/Dets/README.rst | 5 +- src/Dets/diagonalize_CI_SC2.irp.f | 3 +- src/Dets/filter_connected.irp.f | 4 +- src/Dets/save_for_casino.irp.f | 2 +- src/MOs/ASSUMPTIONS.rst | 4 + src/MOs/README.rst | 2 + src/Molden/print_mo.irp.f | 4 +- src/MonoInts/kin_ao_ints.irp.f | 10 +- src/MonoInts/pot_ao_ints.irp.f | 16 +- src/MonoInts/spread_dipole_ao.irp.f | 30 +-- src/Properties/delta_rho.irp.f | 16 +- 18 files changed, 323 insertions(+), 286 deletions(-) create mode 100644 ocaml/.merlin diff --git a/ocaml/.merlin b/ocaml/.merlin new file mode 100644 index 00000000..3683fed6 --- /dev/null +++ b/ocaml/.merlin @@ -0,0 +1,4 @@ +PKG core ZMQ cryptokit +B _build/ + + diff --git a/src/AOs/ASSUMPTIONS.rst b/src/AOs/ASSUMPTIONS.rst index 0c8e0ccc..cede1de9 100644 --- a/src/AOs/ASSUMPTIONS.rst +++ b/src/AOs/ASSUMPTIONS.rst @@ -1,8 +1 @@ -* The atomic orbitals are normalized: - - .. math:: - - \int \left(\chi_i({\bf r}) \right)^2 dr = 1 - * The AO coefficients in the EZFIO files are not necessarily normalized and are normalized after reading -* The AO coefficients and exponents are ordered in increasing order of exponents diff --git a/src/AOs/README.rst b/src/AOs/README.rst index f9f81f5f..5978e16f 100644 --- a/src/AOs/README.rst +++ b/src/AOs/README.rst @@ -17,21 +17,19 @@ The AO coefficients are normalized as: {\tilde c}_{ki} = \frac{c_{ki}}{ \int \left( (x-X_A)^a (y-Y_A)^b (z-Z_A)^c e^{-\gamma_{ki} |{\bf r} - {\bf R}_A|^2} \right)^2} dr +Warning: ``ao_coef`` contains the AO coefficients given in input. These do not +include the normalization constant of the AO. The ``ao_coef_normalized`` includes +this normalization factor. + +The AOs are also sorted by increasing exponent to accelerate the calculation of +the two electron integrals. + Assumptions =========== .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* The atomic orbitals are normalized: - - .. math:: - - \int \left(\chi_i({\bf r}) \right)^2 dr = 1 - -* The AO coefficients in the EZFIO files are not necessarily normalized and are normalized after reading -* The AO coefficients and exponents are ordered in increasing order of exponents - Needed Modules ============== @@ -70,45 +68,41 @@ Documentation Overlap between atomic basis functions: :math:`\int \chi_i(r) \chi_j(r) dr)` -`ao_coef `_ - Coefficients, exponents and powers of x,y and z +`ao_coef `_ + AO Coefficients, read from input. Those should not be used directly, as + the MOs are expressed on the basis of **normalized** AOs. -`ao_coef_transp `_ - Transposed ao_coef and ao_expo +`ao_coef_normalized `_ + Coefficients including the AO normalization -`ao_coef_unnormalized `_ - Coefficients, exponents and powers of x,y and z as in the EZFIO file - ao_coef(i,j) = coefficient of the jth primitive on the ith ao +`ao_coef_normalized_ordered `_ + Sorted primitives to accelerate 4 index MO transformation + +`ao_coef_normalized_ordered_transp `_ + Transposed ao_coef_normalized_ordered + +`ao_expo `_ + AO Exponents read from input + +`ao_expo_ordered `_ + Sorted primitives to accelerate 4 index MO transformation + +`ao_expo_ordered_transp `_ + Transposed ao_expo_ordered + +`ao_l `_ ao_l = l value of the AO: a+b+c in x^a y^b z^c -`ao_expo `_ - Coefficients, exponents and powers of x,y and z - -`ao_expo_transp `_ - Transposed ao_coef and ao_expo - -`ao_expo_unsorted `_ - Coefficients, exponents and powers of x,y and z as in the EZFIO file - ao_coef(i,j) = coefficient of the jth primitive on the ith ao +`ao_l_char `_ ao_l = l value of the AO: a+b+c in x^a y^b z^c -`ao_l `_ - Coefficients, exponents and powers of x,y and z as in the EZFIO file - ao_coef(i,j) = coefficient of the jth primitive on the ith ao - ao_l = l value of the AO: a+b+c in x^a y^b z^c - -`ao_l_char `_ - Coefficients, exponents and powers of x,y and z as in the EZFIO file - ao_coef(i,j) = coefficient of the jth primitive on the ith ao - ao_l = l value of the AO: a+b+c in x^a y^b z^c - -`ao_l_char_space `_ +`ao_l_char_space `_ Undocumented -`ao_md5 `_ +`ao_md5 `_ MD5 key characteristic of the AO basis -`ao_nucl `_ +`ao_nucl `_ Index of the nuclei on which the ao is centered `ao_num `_ @@ -118,35 +112,35 @@ Documentation Number of atomic orbitals `ao_power `_ - Coefficients, exponents and powers of x,y and z + Powers of x,y and z read from input -`ao_prim_num `_ +`ao_prim_num `_ Number of primitives per atomic orbital -`ao_prim_num_max `_ +`ao_prim_num_max `_ Undocumented -`ao_prim_num_max_align `_ +`ao_prim_num_max_align `_ Undocumented -`l_to_charater `_ +`l_to_charater `_ character corresponding to the "L" value of an AO orbital -`n_aos_max `_ +`n_aos_max `_ Number of AOs per atom -`nucl_aos `_ +`nucl_aos `_ List of AOs attached on each atom -`nucl_list_shell_aos `_ +`nucl_list_shell_aos `_ Index of the shell type Aos and of the corresponding Aos Per convention, for P,D,F and G AOs, we take the index of the AO with the the corresponding power in the "X" axis -`nucl_n_aos `_ +`nucl_n_aos `_ Number of AOs per atom -`nucl_num_shell_aos `_ +`nucl_num_shell_aos `_ Index of the shell type Aos and of the corresponding Aos Per convention, for P,D,F and G AOs, we take the index of the AO with the the corresponding power in the "X" axis diff --git a/src/AOs/ao_overlap.irp.f b/src/AOs/ao_overlap.irp.f index f77924d3..737f03f7 100644 --- a/src/AOs/ao_overlap.irp.f +++ b/src/AOs/ao_overlap.irp.f @@ -21,8 +21,8 @@ !$OMP overlap_x,overlap_y, overlap_z, overlap, & !$OMP alpha, beta,i,j,c) & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_transp,ao_nucl, & - !$OMP ao_expo_transp,dim1) + !$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & + !$OMP ao_expo_ordered_transp,dim1) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(2) = nucl_coord( ao_nucl(j), 2 ) @@ -44,12 +44,12 @@ power_B(2) = ao_power( i, 2 ) power_B(3) = ao_power( i, 3 ) do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) !DEC$ VECTOR ALIGNED do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - c = ao_coef_transp(n,j) * ao_coef_transp(l,i) + c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) ao_overlap(i,j) += c * overlap ao_overlap_x(i,j) += c * overlap_x ao_overlap_y(i,j) += c * overlap_y @@ -84,8 +84,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ] !$OMP overlap_x,overlap_y, overlap_z, overlap, & !$OMP alpha, beta,i,j,dx) & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_overlap_abs,ao_num,ao_coef_transp,ao_nucl, & - !$OMP ao_expo_transp,dim1,lower_exp_val) + !$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & + !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(2) = nucl_coord( ao_nucl(j), 2 ) @@ -104,14 +104,14 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ] power_B(2) = ao_power( i, 2 ) power_B(3) = ao_power( i, 3 ) do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) !DEC$ VECTOR ALIGNED do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1) call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1) call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1) - ao_overlap_abs(i,j) += abs(ao_coef_transp(n,j) * ao_coef_transp(l,i)) * overlap_x * overlap_y * overlap_z + ao_overlap_abs(i,j) += abs(ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z enddo enddo enddo diff --git a/src/AOs/aos.irp.f b/src/AOs/aos.irp.f index c75694fb..b9c82327 100644 --- a/src/AOs/aos.irp.f +++ b/src/AOs/aos.irp.f @@ -1,152 +1,171 @@ BEGIN_PROVIDER [ integer, ao_num ] &BEGIN_PROVIDER [ integer, ao_num_align ] - implicit none - - BEGIN_DOC -! Number of atomic orbitals - END_DOC - - ao_num = -1 - PROVIDE ezfio_filename - call ezfio_get_ao_basis_ao_num(ao_num) - if (ao_num <= 0) then - stop 'Number of contracted gaussians should be > 0' - endif - integer :: align_double - ao_num_align = align_double(ao_num) + implicit none + + BEGIN_DOC + ! Number of atomic orbitals + END_DOC + + ao_num = -1 + PROVIDE ezfio_filename + call ezfio_get_ao_basis_ao_num(ao_num) + if (ao_num <= 0) then + stop 'Number of contracted gaussians should be > 0' + endif + integer :: align_double + ao_num_align = align_double(ao_num) +END_PROVIDER + +BEGIN_PROVIDER [ integer, ao_power, (ao_num_align,3) ] + implicit none + BEGIN_DOC + ! Powers of x,y and z read from input + END_DOC + PROVIDE ezfio_filename + + integer :: i,j,k + integer, allocatable :: ibuffer(:,:) + allocate ( ibuffer(ao_num,3) ) + ibuffer = 0 + call ezfio_get_ao_basis_ao_power(ibuffer) + ao_power = 0 + do j = 1, 3 + do i = 1, ao_num + ao_power(i,j) = ibuffer(i,j) + enddo + enddo + deallocate(ibuffer) + END_PROVIDER - BEGIN_PROVIDER [ integer, ao_power, (ao_num_align,3) ] -&BEGIN_PROVIDER [ double precision, ao_expo, (ao_num_align,ao_prim_num_max) ] -&BEGIN_PROVIDER [ double precision, ao_coef, (ao_num_align,ao_prim_num_max) ] - implicit none - - BEGIN_DOC -! Coefficients, exponents and powers of x,y and z - END_DOC - PROVIDE ezfio_filename - - double precision, allocatable :: buffer(:,:) - allocate ( buffer(ao_num,ao_prim_num_max) ) - integer :: ibuffer(ao_num,3) - integer :: i,j,k - character*(128) :: give_ao_character_space - ibuffer = 0 - call ezfio_get_ao_basis_ao_power(ibuffer) - ao_power = 0 - do j = 1, 3 - do i = 1, ao_num - ao_power(i,j) = ibuffer(i,j) +BEGIN_PROVIDER [ double precision, ao_expo, (ao_num_align,ao_prim_num_max) ] + implicit none + BEGIN_DOC + ! AO Exponents read from input + END_DOC + PROVIDE ezfio_filename + + double precision, allocatable :: buffer(:,:) + allocate ( buffer(ao_num,ao_prim_num_max) ) + integer :: i,j,k + ao_expo = 0.d0 + buffer = 0.d0 + call ezfio_get_ao_basis_ao_expo(buffer) + do j = 1, ao_prim_num_max + do i = 1, ao_num + ao_expo(i,j) = buffer(i,j) + enddo enddo - enddo - ao_expo = 0.d0 - buffer = 0.d0 - call ezfio_get_ao_basis_ao_expo(buffer) - do j = 1, ao_prim_num_max - do i = 1, ao_num - ao_expo(i,j) = buffer(i,j) - enddo - enddo + deallocate(buffer) +END_PROVIDER - ao_coef = 0.d0 - buffer = 0.d0 - call ezfio_get_ao_basis_ao_coef(buffer) - do j = 1, ao_prim_num_max - do i = 1, ao_num - ao_coef(i,j) = buffer(i,j) +BEGIN_PROVIDER [ double precision, ao_coef, (ao_num_align,ao_prim_num_max) ] + implicit none + BEGIN_DOC + ! AO Coefficients, read from input. Those should not be used directly, as + ! the MOs are expressed on the basis of **normalized** AOs. + END_DOC + PROVIDE ezfio_filename + + double precision, allocatable :: buffer(:,:) + allocate ( buffer(ao_num,ao_prim_num_max) ) + integer :: i,j,k + ao_coef = 0.d0 + buffer = 0.d0 + call ezfio_get_ao_basis_ao_coef(buffer) + do j = 1, ao_prim_num_max + do i = 1, ao_num + ao_coef(i,j) = buffer(i,j) + enddo enddo - enddo + deallocate(buffer) +END_PROVIDER - deallocate(buffer) +BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num_align,ao_prim_num_max) ] + implicit none + BEGIN_DOC + ! Coefficients including the AO normalization + END_DOC + double precision :: norm, norm2,overlap_x,overlap_y,overlap_z,C_A(3) + integer :: l, powA(3), nz + integer :: i,j + nz=100 + C_A(1) = 0.d0 + C_A(2) = 0.d0 + C_A(3) = 0.d0 + do i=1,ao_num + powA(1) = ao_power(i,1) + powA(2) = ao_power(i,2) + powA(3) = ao_power(i,3) + do j=1,ao_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm) + enddo + enddo +END_PROVIDER -! Normalization of the AO coefficients -! ------------------------------------ - double precision :: norm, norm2,overlap_x,overlap_y,overlap_z,C_A(3) - integer :: l, powA(3), nz - nz=100 - C_A(1) = 0.d0 - C_A(2) = 0.d0 - C_A(3) = 0.d0 - do i=1,ao_num - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) - do j=1,ao_prim_num(i) - call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) - ao_coef(i,j) = ao_coef(i,j)/sqrt(norm) - enddo - enddo - - ! Sorting of the exponents for efficient integral calculations - integer :: iorder(ao_prim_num_max) - double precision :: d(ao_prim_num_max,2) - do i=1,ao_num - do j=1,ao_prim_num(i) - iorder(j) = j - d(j,1) = ao_expo(i,j) - d(j,2) = ao_coef(i,j) - enddo - call dsort(d(1,1),iorder,ao_prim_num(i)) - call dset_order(d(1,2),iorder,ao_prim_num(i)) - do j=1,ao_prim_num(i) - ao_expo(i,j) = d(j,1) - ao_coef(i,j) = d(j,2) - enddo - enddo + BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num_align,ao_prim_num_max) ] +&BEGIN_PROVIDER [ double precision, ao_expo_ordered, (ao_num_align,ao_prim_num_max) ] + implicit none + BEGIN_DOC + ! Sorted primitives to accelerate 4 index MO transformation + END_DOC + + integer :: iorder(ao_prim_num_max) + double precision :: d(ao_prim_num_max,2) + integer :: i,j + do i=1,ao_num + do j=1,ao_prim_num(i) + iorder(j) = j + d(j,1) = ao_expo(i,j) + d(j,2) = ao_coef_normalized(i,j) + enddo + call dsort(d(1,1),iorder,ao_prim_num(i)) + call dset_order(d(1,2),iorder,ao_prim_num(i)) + do j=1,ao_prim_num(i) + ao_expo_ordered(i,j) = d(j,1) + ao_coef_normalized_ordered(i,j) = d(j,2) + enddo + enddo END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ] - implicit none - BEGIN_DOC -! Transposed ao_coef and ao_expo - END_DOC - integer :: i,j - do j=1, ao_num - do i=1, ao_prim_num_max - ao_coef_transp(i,j) = ao_coef(j,i) - ao_expo_transp(i,j) = ao_expo(j,i) +BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered_transp, (ao_prim_num_max_align,ao_num) ] + implicit none + BEGIN_DOC + ! Transposed ao_coef_normalized_ordered + END_DOC + integer :: i,j + do j=1, ao_num + do i=1, ao_prim_num_max + ao_coef_normalized_ordered_transp(i,j) = ao_coef_normalized_ordered(j,i) + enddo enddo - enddo - - + END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_coef_unnormalized, (ao_num_align,ao_prim_num_max) ] -&BEGIN_PROVIDER [ double precision, ao_expo_unsorted, (ao_num_align,ao_prim_num_max) ] -&BEGIN_PROVIDER [ integer, ao_l, (ao_num) ] +BEGIN_PROVIDER [ double precision, ao_expo_ordered_transp, (ao_prim_num_max_align,ao_num) ] + implicit none + BEGIN_DOC + ! Transposed ao_expo_ordered + END_DOC + integer :: i,j + do j=1, ao_num + do i=1, ao_prim_num_max + ao_expo_ordered_transp(i,j) = ao_expo_ordered(j,i) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ integer, ao_l, (ao_num) ] &BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ] implicit none - BEGIN_DOC -! Coefficients, exponents and powers of x,y and z as in the EZFIO file -! ao_coef(i,j) = coefficient of the jth primitive on the ith ao ! ao_l = l value of the AO: a+b+c in x^a y^b z^c END_DOC - PROVIDE ezfio_filename - - double precision, allocatable :: buffer(:,:) - allocate ( buffer(ao_num,ao_prim_num_max) ) - integer :: i,j,k - character*(128) :: give_ao_character_space - buffer = 0.d0 - call ezfio_get_ao_basis_ao_expo(buffer) - do j = 1, ao_prim_num_max - do i = 1, ao_num - ao_expo_unsorted(i,j) = buffer(i,j) - enddo - enddo - - buffer = 0.d0 - call ezfio_get_ao_basis_ao_coef(buffer) - do j = 1, ao_prim_num_max - do i = 1, ao_num - ao_coef_unnormalized(i,j) = buffer(i,j) - enddo - enddo - deallocate(buffer) - + integer :: i do i=1,ao_num ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) ao_l_char(i) = l_to_charater(ao_l(i)) @@ -154,23 +173,6 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ] - implicit none - BEGIN_DOC -! Transposed ao_coef and ao_expo - END_DOC - integer :: i,j - do j=1, ao_num - do i=1, ao_prim_num_max - ao_coef_transp(i,j) = ao_coef(j,i) - ao_expo_transp(i,j) = ao_expo(j,i) - enddo - enddo - - -END_PROVIDER - BEGIN_PROVIDER [ integer, ao_prim_num, (ao_num_align) ] implicit none @@ -303,10 +305,10 @@ END_PROVIDER enddo enddo - END_PROVIDER +END_PROVIDER - BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] +BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] implicit none integer :: i character*(4) :: give_ao_character_space @@ -397,6 +399,7 @@ END_PROVIDER ao_l_char_space(i) = give_ao_character_space enddo END_PROVIDER + BEGIN_PROVIDER [ character*(32), ao_md5 ] BEGIN_DOC ! MD5 key characteristic of the AO basis diff --git a/src/Bielec_integrals/ao_bi_integrals.irp.f b/src/Bielec_integrals/ao_bi_integrals.irp.f index 0da76021..17cf4afc 100644 --- a/src/Bielec_integrals/ao_bi_integrals.irp.f +++ b/src/Bielec_integrals/ao_bi_integrals.irp.f @@ -42,24 +42,24 @@ double precision function ao_bielec_integral(i,j,k,l) do p = 1, ao_prim_num(i) double precision :: coef1 - coef1 = ao_coef_transp(p,i) + coef1 = ao_coef_normalized_ordered_transp(p,i) do q = 1, ao_prim_num(j) double precision :: coef2 - coef2 = coef1*ao_coef_transp(q,j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) double precision :: p_inv,q_inv call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_transp(p,i),ao_expo_transp(q,j), & + 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) double precision :: coef3 - coef3 = coef2*ao_coef_transp(r,k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) do s = 1, ao_prim_num(l) double precision :: coef4 - coef4 = coef3*ao_coef_transp(s,l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) double precision :: general_primitive_integral call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_transp(r,k),ao_expo_transp(s,l), & + 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, & @@ -82,15 +82,15 @@ double precision function ao_bielec_integral(i,j,k,l) double precision :: ERI do p = 1, ao_prim_num(i) - coef1 = ao_coef_transp(p,i) + coef1 = ao_coef_normalized_ordered_transp(p,i) do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_transp(q,j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_transp(r,k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_transp(s,l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) integral = ERI( & - ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(r,k),ao_expo_transp(s,l),& + 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)) @@ -149,12 +149,12 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) schwartz_kl(0,0) = 0.d0 do r = 1, ao_prim_num(k) - coef1 = ao_coef_transp(r,k)*ao_coef_transp(r,k) + coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) schwartz_kl(0,r) = 0.d0 do s = 1, ao_prim_num(l) - coef2 = coef1 * ao_coef_transp(s,l) * ao_coef_transp(s,l) + coef2 = coef1 * ao_coef_normalized_ordered_transp(s,l) * ao_coef_normalized_ordered_transp(s,l) call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_transp(r,k),ao_expo_transp(s,l), & + 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 schwartz_kl(s,r) = general_primitive_integral(dim1, & @@ -168,13 +168,13 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) do p = 1, ao_prim_num(i) double precision :: coef1 - coef1 = ao_coef_transp(p,i) + coef1 = ao_coef_normalized_ordered_transp(p,i) do q = 1, ao_prim_num(j) double precision :: coef2 - coef2 = coef1*ao_coef_transp(q,j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) double precision :: p_inv,q_inv call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_transp(p,i),ao_expo_transp(q,j), & + 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 schwartz_ij = general_primitive_integral(dim1, & @@ -189,16 +189,16 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) cycle endif double precision :: coef3 - coef3 = coef2*ao_coef_transp(r,k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) do s = 1, ao_prim_num(l) double precision :: coef4 if (schwartz_kl(s,r)*schwartz_ij < thresh) then cycle endif - coef4 = coef3*ao_coef_transp(s,l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) double precision :: general_primitive_integral call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_transp(r,k),ao_expo_transp(s,l), & + 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, & @@ -222,12 +222,12 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) schwartz_kl(0,0) = 0.d0 do r = 1, ao_prim_num(k) - coef1 = ao_coef_transp(r,k)*ao_coef_transp(r,k) + coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) schwartz_kl(0,r) = 0.d0 do s = 1, ao_prim_num(l) - coef2 = coef1*ao_coef_transp(s,l)*ao_coef_transp(s,l) + coef2 = coef1*ao_coef_normalized_ordered_transp(s,l)*ao_coef_normalized_ordered_transp(s,l) schwartz_kl(s,r) = ERI( & - ao_expo_transp(r,k),ao_expo_transp(s,l),ao_expo_transp(r,k),ao_expo_transp(s,l),& + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& K_power(1),L_power(1),K_power(1),L_power(1), & K_power(2),L_power(2),K_power(2),L_power(2), & K_power(3),L_power(3),K_power(3),L_power(3)) * & @@ -238,11 +238,11 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) enddo do p = 1, ao_prim_num(i) - coef1 = ao_coef_transp(p,i) + coef1 = ao_coef_normalized_ordered_transp(p,i) do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_transp(q,j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) schwartz_ij = ERI( & - ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(p,i),ao_expo_transp(q,j),& + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),& I_power(1),J_power(1),I_power(1),J_power(1), & I_power(2),J_power(2),I_power(2),J_power(2), & I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2 @@ -253,14 +253,14 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) if (schwartz_kl(0,r)*schwartz_ij < thresh) then cycle endif - coef3 = coef2*ao_coef_transp(r,k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) do s = 1, ao_prim_num(l) if (schwartz_kl(s,r)*schwartz_ij < thresh) then cycle endif - coef4 = coef3*ao_coef_transp(s,l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) integral = ERI( & - ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(r,k),ao_expo_transp(s,l),& + 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)) diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index 801d00a5..d230c765 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -181,8 +181,43 @@ subroutine copy_H_apply_buffer_to_wf call normalize(psi_coef,N_det) SOFT_TOUCH N_det psi_det psi_coef + call debug_unicity_of_determinants end +subroutine debug_unicity_of_determinants + implicit none + BEGIN_DOC +! This subroutine checks that there are no repetitions in the wave function + END_DOC + logical :: same, failed + integer :: i,k + print *, "======= DEBUG UNICITY =========" + failed = .False. + do i=2,N_det + same = .True. + do k=1,N_int + if ( psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,i-1) ) then + same = .False. + exit + endif + if ( psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,i-1) ) then + same = .False. + exit + endif + enddo + if (same) then + failed = .True. + call debug_det(psi_det_sorted_bit(1,1,i)) + endif + enddo + + if (failed) then + print *, '======= Determinants not unique : Failed ! =========' + stop + else + print *, '======= Determinants are unique : OK ! =========' + endif +end subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) use bitmasks diff --git a/src/Dets/README.rst b/src/Dets/README.rst index e9077510..cbc1c352 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -54,7 +54,10 @@ Documentation after calling this function. After calling this subroutine, N_det, psi_det and psi_coef need to be touched -`fill_h_apply_buffer_no_selection `_ +`debug_unicity_of_determinants `_ + This subroutine checks that there are no repetitions in the wave function + +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD `h_apply_buffer_allocated `_ diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index 86ba72b9..81931a62 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -35,8 +35,7 @@ END_PROVIDER do i=1,N_det CI_SC2_eigenvectors(i,j) = psi_coef(i,j) enddo -! TODO : check comment -! CI_SC2_electronic_energy(j) = CI_electronic_energy(j) + CI_SC2_electronic_energy(j) = CI_electronic_energy(j) enddo call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f index 93a6ee7b..c3d333ad 100644 --- a/src/Dets/filter_connected.irp.f +++ b/src/Dets/filter_connected.irp.f @@ -235,8 +235,8 @@ subroutine filter_connected_davidson(key1,key2,Nint,sze,idx) else if (Nint==3) then - !DIR$ LOOP COUNT (1000) i = idx(0) + !DIR$ LOOP COUNT (1000) do j_int=1,N_con_int itmp = det_connections(j_int,i) do while (itmp /= 0_8) @@ -261,8 +261,8 @@ subroutine filter_connected_davidson(key1,key2,Nint,sze,idx) else - !DIR$ LOOP COUNT (1000) i = idx(0) + !DIR$ LOOP COUNT (1000) do j_int=1,N_con_int itmp = det_connections(j_int,i) do while (itmp /= 0_8) diff --git a/src/Dets/save_for_casino.irp.f b/src/Dets/save_for_casino.irp.f index 631f79bd..35c0c3a7 100644 --- a/src/Dets/save_for_casino.irp.f +++ b/src/Dets/save_for_casino.irp.f @@ -161,7 +161,7 @@ subroutine save_casino if (ao_l(i) == ao_power(i,1)) then do j=1,ao_prim_num(i) l+=1 - rtmp(l) = ao_coef(i,ao_prim_num(i)-j+1) + rtmp(l) = ao_coef_normalized(i,ao_prim_num(i)) enddo endif enddo diff --git a/src/MOs/ASSUMPTIONS.rst b/src/MOs/ASSUMPTIONS.rst index e69de29b..8514ef73 100644 --- a/src/MOs/ASSUMPTIONS.rst +++ b/src/MOs/ASSUMPTIONS.rst @@ -0,0 +1,4 @@ +ASSUMPTONS +========== + +* The AO basis functions are normalized. diff --git a/src/MOs/README.rst b/src/MOs/README.rst index d7a15219..d7a869b4 100644 --- a/src/MOs/README.rst +++ b/src/MOs/README.rst @@ -8,6 +8,8 @@ Molecular orbitals are expressed as \phi_k({\bf r}) = \sum_i C_{ik} \chi_k({\bf r}) +where :math:`\chi_k` are *normalized* atomic basis set. + The current set of molecular orbitals has a label ``mo_label``. When the orbitals are modified, the label should also be updated to keep everything consistent. diff --git a/src/Molden/print_mo.irp.f b/src/Molden/print_mo.irp.f index 9cec8fbd..b147fe50 100644 --- a/src/Molden/print_mo.irp.f +++ b/src/Molden/print_mo.irp.f @@ -92,9 +92,9 @@ subroutine write_Ao_basis(i_unit_output) do k = 1, ao_prim_num(i_ao) i_prim +=1 if(i_prim.lt.100)then - write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k) + write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo(i_ao,k),ao_coef(i_ao,k) else - write(i_unit_output,'(4X,I3,3X,A1,5X,I3,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k) + write(i_unit_output,'(4X,I3,3X,A1,5X,I3,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo(i_ao,k),ao_coef(i_ao,k) endif enddo write(i_unit_output,*)'' diff --git a/src/MonoInts/kin_ao_ints.irp.f b/src/MonoInts/kin_ao_ints.irp.f index 444c3880..10b065b4 100644 --- a/src/MonoInts/kin_ao_ints.irp.f +++ b/src/MonoInts/kin_ao_ints.irp.f @@ -36,8 +36,8 @@ !$OMP alpha, beta,i,j,c,d_a_2,d_2,deriv_tmp, & !$OMP overlap_x0,overlap_y0,overlap_z0) & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_transp,ao_nucl, & - !$OMP ao_expo_transp,dim1) + !$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & + !$OMP ao_expo_ordered_transp,dim1) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(2) = nucl_coord( ao_nucl(j), 2 ) @@ -58,12 +58,12 @@ power_B(2) = ao_power( i, 2 ) power_B(3) = ao_power( i, 3 ) do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) !DEC$ VECTOR ALIGNED do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1) - c = ao_coef_transp(n,j) * ao_coef_transp(l,i) + c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) ! if (abs(c) < 1.d-8) then ! cycle ! endif diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index f430ace9..32e98a2c 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -19,7 +19,7 @@ !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, & !$OMP num_A,num_B,Z,c,n_pt_in) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge) n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (guided) @@ -40,9 +40,9 @@ B_center(2) = nucl_coord(num_B,2) B_center(3) = nucl_coord(num_B,3) do l=1,ao_prim_num(j) - alpha = ao_expo_transp(l,j) + alpha = ao_expo_ordered_transp(l,j) do m=1,ao_prim_num(i) - beta = ao_expo_transp(m,i) + beta = ao_expo_ordered_transp(m,i) c = 0.d0 do k = 1, nucl_num double precision :: Z,c @@ -53,7 +53,7 @@ c = c+Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) enddo ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - & - ao_coef_transp(l,j)*ao_coef_transp(m,i)*c + ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c enddo enddo enddo @@ -90,7 +90,7 @@ END_PROVIDER !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,l,m,alpha,beta,A_center,B_center,power_A,power_B, & !$OMP num_A,num_B,c,n_pt_in) & - !$OMP SHARED (k,ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & + !$OMP SHARED (k,ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & !$OMP n_pt_max_integrals,ao_nucl_elec_integral_per_atom,nucl_num,C_center) n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (guided) @@ -114,11 +114,11 @@ END_PROVIDER B_center(3) = nucl_coord(num_B,3) c = 0.d0 do l=1,ao_prim_num(j) - alpha = ao_expo_transp(l,j) + alpha = ao_expo_ordered_transp(l,j) do m=1,ao_prim_num(i) - beta = ao_expo_transp(m,i) + beta = ao_expo_ordered_transp(m,i) c = c + NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) & - * ao_coef_transp(l,j)*ao_coef_transp(m,i) + * ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i) enddo enddo ao_nucl_elec_integral_per_atom(i,j,k) = -c diff --git a/src/MonoInts/spread_dipole_ao.irp.f b/src/MonoInts/spread_dipole_ao.irp.f index c0d7c88e..d7aa738a 100644 --- a/src/MonoInts/spread_dipole_ao.irp.f +++ b/src/MonoInts/spread_dipole_ao.irp.f @@ -26,8 +26,8 @@ !$OMP overlap_x,overlap_y, overlap_z, overlap, & !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_transp,ao_nucl, & - !$OMP ao_expo_transp,dim1,lower_exp_val) + !$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & + !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(2) = nucl_coord( ao_nucl(j), 2 ) @@ -48,11 +48,11 @@ accu_y = 0.d0 accu_z = 0.d0 do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) !DEC$ VECTOR ALIGNED do l = 1, ao_prim_num(i) - c = ao_coef_transp(n,j)*ao_coef_transp(l,i) - beta = ao_expo_transp(l,i) + c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) accu_x += c*(tmp*overlap_y*overlap_z) @@ -100,8 +100,8 @@ !$OMP overlap_x,overlap_y, overlap_z, overlap, & !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_transp,ao_nucl, & - !$OMP ao_expo_transp,dim1,lower_exp_val) + !$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & + !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(2) = nucl_coord( ao_nucl(j), 2 ) @@ -122,11 +122,11 @@ accu_y = 0.d0 accu_z = 0.d0 do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) !DEC$ VECTOR ALIGNED do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) - c = ao_coef_transp(l,i)*ao_coef_transp(n,j) + beta = ao_expo_ordered_transp(l,i) + c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) @@ -174,8 +174,8 @@ !$OMP overlap_x,overlap_y, overlap_z, overlap, & !$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_transp,ao_nucl, & - !$OMP ao_expo_transp,dim1,lower_exp_val) + !$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & + !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(2) = nucl_coord( ao_nucl(j), 2 ) @@ -196,12 +196,12 @@ accu_y = 0.d0 accu_z = 0.d0 do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) !DEC$ VECTOR ALIGNED do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - c = ao_coef_transp(l,i) * ao_coef_transp(n,j) + c = ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(n,j) i_component = 1 call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) accu_x += c*(tmp*overlap_y*overlap_z) diff --git a/src/Properties/delta_rho.irp.f b/src/Properties/delta_rho.irp.f index 3cfe136c..69894c38 100644 --- a/src/Properties/delta_rho.irp.f +++ b/src/Properties/delta_rho.irp.f @@ -79,7 +79,7 @@ BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_all_points, (ao_num_a !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(i,j,n,l,A_center,power_A,B_center,power_B,accu_z, & !$OMP overlap_x,overlap_y,overlap_z,overlap,c,alpha,beta) & - !$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_transp,ao_coef_transp, & + !$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_ordered_transp,ao_coef_normalized_ordered_transp, & !$OMP ao_integrated_delta_rho_all_points,N_z_pts,dim1,i_z,z,delta_z) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) @@ -98,12 +98,12 @@ BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_all_points, (ao_num_a accu_z = 0.d0 do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - c = ao_coef_transp(n,j) * ao_coef_transp(l,i) + c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) accu_z += c* overlap_x * overlap_y * SABpartial(z,z+delta_z,A_center,B_center,power_A,power_B,alpha,beta) enddo enddo @@ -147,7 +147,7 @@ BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_one_point, (ao_num_al !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(i,j,n,l,A_center,power_A,B_center,power_B,accu_z, & !$OMP overlap_x,overlap_y,overlap_z,overlap,c,alpha,beta) & - !$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_transp,ao_coef_transp, & + !$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_ordered_transp,ao_coef_normalized_ordered_transp, & !$OMP ao_integrated_delta_rho_one_point,dim1,z,delta_z) do j=1,ao_num A_center(1) = nucl_coord( ao_nucl(j), 1 ) @@ -166,12 +166,12 @@ BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_one_point, (ao_num_al accu_z = 0.d0 do n = 1,ao_prim_num(j) - alpha = ao_expo_transp(n,j) + alpha = ao_expo_ordered_transp(n,j) do l = 1, ao_prim_num(i) - beta = ao_expo_transp(l,i) + beta = ao_expo_ordered_transp(l,i) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - c = ao_coef_transp(n,j) * ao_coef_transp(l,i) + c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) accu_z += c* overlap_x * overlap_y * SABpartial(z,z+delta_z,A_center,B_center,power_A,power_B,alpha,beta) enddo enddo