diff --git a/src/AO/ao.irp.f b/src/AO/ao.irp.f new file mode 100644 index 0000000..9fd2590 --- /dev/null +++ b/src/AO/ao.irp.f @@ -0,0 +1,365 @@ + BEGIN_PROVIDER [ integer, ao_num ] +&BEGIN_PROVIDER [ integer, ao_num_8 ] + implicit none + BEGIN_DOC +! Number of atomic orbitals + END_DOC + integer, external :: mod_align + + ao_num = -1 + call get_ao_basis_ao_num(ao_num) + if (ao_num <= 0) then + call abrt(irp_here,'Number of contracted gaussians should be > 0') + endif + call iinfo(irp_here,'ao_num',ao_num) + ao_num_8 = mod_align(ao_num) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, ao_prim_num, (ao_num_8) ] + implicit none + BEGIN_DOC +! Number of primitives per atomic orbital + END_DOC + + ao_prim_num = 0 + call get_ao_basis_ao_prim_num(ao_prim_num) + integer :: i + character*(80) :: message + do i=1,ao_num + if (ao_prim_num(i) <= 0) then + write(message,'(A,I6,A)') 'Number of primitives of contraction ',i,' should be > 0' + call abrt(irp_here,message) + endif + enddo + call iset_order(ao_prim_num,ao_nucl_sort_idx,ao_num) + call iinfo(irp_here,'sum_ao_prim_num',sum(ao_prim_num)) + +END_PROVIDER + + + BEGIN_PROVIDER [ integer, ao_nucl, (ao_num_8) ] +&BEGIN_PROVIDER [ integer, ao_nucl_sort_idx, (ao_num_8) ] +&BEGIN_PROVIDER [ integer, ao_nucl_idx, (2,nucl_num) ] + implicit none + BEGIN_DOC +! Nucleus on which the atomic orbital is centered + END_DOC + + ao_nucl = -1 + call get_ao_basis_ao_nucl(ao_nucl) + + character*(80) :: message + character*(30) :: range + integer :: i,k + + write(range,'(A,I5,A)') '(1,',nucl_num,')' + + do i=1,ao_num + if ( (ao_nucl(i) <= 0) .or. (ao_nucl(i) > nucl_num) ) then + write(message,'(A,I6,A)') 'Contraction ',i,' should be centered on a nucleus in the range'//trim(range) + call abrt(irp_here,message) + endif + ao_nucl_sort_idx(i) = i + enddo + call isort(ao_nucl,ao_nucl_sort_idx,ao_num) + ao_nucl_idx(1,ao_nucl(1)) = 1 + ao_nucl_idx(2,ao_nucl(1)) = 1 + do i=2,ao_num + k = ao_nucl(i) + if (k == ao_nucl(i-1)) then + ao_nucl_idx(2,k) = i + else + ao_nucl_idx(1,k) = i + ao_nucl_idx(2,k) = i + endif + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ integer, ao_power, (ao_num,4) ] +&BEGIN_PROVIDER [ logical, ao_power_is_zero, (ao_num) ] + implicit none + BEGIN_DOC +! x,y,z powers of the atomic orbital +! +! 4 contains the sum of powers +! +! ao_power_is_zero is true where ao_power(:,4) == 0 + END_DOC + ao_power = 0 + call get_ao_basis_ao_power(ao_power) + + character*(80) :: message + integer :: i,j + do i=1,3 + do j=1,ao_num + if (ao_power(j,i) < 0) then + write(message,'(A,I1,A,I6,A)') 'Power ',i,' of contraction ',j,' should be > 0' + call abrt(irp_here,message) + endif + enddo + call iset_order(ao_power(1,i),ao_nucl_sort_idx,ao_num) + enddo + do j=1,ao_num + ao_power(j,4) = ao_power(j,1)+ao_power(j,2)+ao_power(j,3) + ao_power_is_zero(j) = ao_power(j,4) == 0 + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ integer, ao_power_transp, (4,ao_num) ] + implicit none + BEGIN_DOC +! Transposed ao_power + END_DOC + integer :: i,j + do i=1,ao_num + do j=1,4 + ao_power_transp(j,i) = ao_power(i,j) + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer , ao_class, (ao_num_8) ] + implicit none + include '../types.F' + BEGIN_DOC +! t_Slater or t_Gaussian + END_DOC + + integer :: i + do i=1,ao_num + ao_class(i) = t_Gaussian + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ integer , ao_power_max ] + implicit none + BEGIN_DOC +! Maximum power among x, y and z + END_DOC + ao_power_max = maxval(ao_power_max_nucl) +END_PROVIDER + + +BEGIN_PROVIDER [ integer , ao_power_max_nucl, (nucl_num_8,4) ] + implicit none + BEGIN_DOC +! Maximum powers of x, y and z per nucleus + END_DOC + integer :: i, j + ao_power_max_nucl = 0 + + integer :: inucl + do j=1,3 + do i=1,ao_num + inucl = ao_nucl(i) + ao_power_max_nucl(inucl,j) = max(ao_power(i,j),ao_power_max_nucl(inucl,j)) + enddo + enddo + do inucl=1,nucl_num + ao_power_max_nucl(inucl,4) = max(ao_power_max_nucl(inucl,1), & + max(ao_power_max_nucl(inucl,2),ao_power_max_nucl(inucl,3)) ) + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ integer, ao_prim_num_max ] +&BEGIN_PROVIDER [ integer, ao_prim_num_max_2 ] +&BEGIN_PROVIDER [ integer, ao_prim_num_max_8 ] +&BEGIN_PROVIDER [ integer, ao_prim_num_max_pow2 ] + implicit none + BEGIN_DOC +! Max Number of primitives per atomic orbital + END_DOC + integer, external :: mod_align + + ao_prim_num_max = maxval(ao_prim_num) + call iinfo(irp_here,'ao_prim_num_max',ao_prim_num_max) + ao_prim_num_max_8 = mod_align(ao_prim_num_max) + ao_prim_num_max_2 = mod(ao_prim_num_max,2) + if (ao_prim_num_max_2 == 0) then + ao_prim_num_max_2 = ao_prim_num_max + else + ao_prim_num_max_2 = ao_prim_num_max + 2 - ao_prim_num_max_2 + endif + ao_prim_num_max_pow2 = 16 + do while (ao_prim_num_max_pow2 < ao_prim_num_max) + ao_prim_num_max_pow2 *= 2 + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ real, ao_expo, (ao_num_8,ao_prim_num_max) ] +&BEGIN_PROVIDER [ real, ao_coef, (ao_num_8,ao_prim_num_max) ] + implicit none + BEGIN_DOC +! Exponents and coefficients of the atomic orbitals + END_DOC + + ao_expo = 0. + real, allocatable :: buf(:,:) + allocate (buf(ao_num,ao_prim_num_max)) + buf = 0. + call get_ao_basis_ao_expo(buf) + do j=1,ao_prim_num_max + do i=1,ao_num + ao_expo(i,j) = buf(i,j) + enddo + call set_order(ao_expo(1,j),ao_nucl_sort_idx,ao_num) + enddo + + integer :: i,j + do i=1,ao_num + do j=1,ao_prim_num(i) + if (ao_expo(i,j) <= 0.) then + character*(80) :: message + write(message,'(A,I6,A,I6,A)') 'Exponent ',j,' of contracted gaussian ',i,' is < 0' + call abrt(irp_here,message) + endif + enddo + enddo + + ao_coef = 0. + buf = 0. + call get_ao_basis_ao_coef(buf) + do j=1,ao_prim_num_max + do i=1,ao_num + ao_coef(i,j) = buf(i,j) + enddo + call set_order(ao_coef(1,j),ao_nucl_sort_idx,ao_num) + enddo + deallocate(buf) + + real,allocatable :: buffer(:) + integer, allocatable :: order(:) + allocate (buffer(ao_prim_num_max),order(ao_prim_num_max)) + do i=1,ao_num + do j=1,ao_prim_num(i) + buffer(j) = ao_expo(i,j) + order(j) = j + enddo + call sort(buffer,order,ao_prim_num(i)) + do j=1,ao_prim_num(i) + ao_expo(i,j) = buffer(j) + buffer(j) = ao_coef(i,j) + enddo + do j=1,ao_prim_num(i) + ao_coef(i,order(j)) = buffer(j) + enddo + enddo + deallocate(buffer,order) + +! Normalization of the AO coefficients +! ------------------------------------ + real :: norm, norm2 + real :: goverlap + integer :: pow(3), l + do i=1,ao_num + do j=1,ao_prim_num(i) + pow(1) = ao_power_transp(1,i) + pow(2) = ao_power_transp(2,i) + pow(3) = ao_power_transp(3,i) + norm = goverlap(ao_expo(i,j),ao_expo(i,j),pow) + ao_coef(i,j) = ao_coef(i,j)/sqrt(norm) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ real, ao_expo_unique, (2*sum(ao_prim_num)) ] +&BEGIN_PROVIDER [ integer, ao_expo_unique_nucl_idx, (2,nucl_num) ] +&BEGIN_PROVIDER [ integer, ao_expo_unique_idx_1, (ao_num_8) ] +&BEGIN_PROVIDER [ integer, ao_expo_unique_idx, (ao_prim_num_max_pow2,ao_num) ] + implicit none + BEGIN_DOC +! Exponents and coefficients of the atomic orbitals + END_DOC + + integer :: k, kk, kkstart, kk2, i, j + integer, allocatable :: order(:) + allocate ( order(2*sum(ao_prim_num)) ) + + ao_expo_unique_idx = 1 + kk=1 + do k=1,nucl_num + kkstart = kk + do i=1,ao_num + if (ao_nucl(i) == k) then + do j=1,ao_prim_num(i) + order(kk) = kk-kkstart+1 + ao_expo_unique(kk) = ao_expo(i,j) + kk += 1 + enddo + endif + enddo + call sort(ao_expo_unique(kkstart),order(kkstart),kk-kkstart) + order(kk) = 0 + kk = kkstart + kk2 = kk+1 + do while (order(kk2) /= 0) + if ( ao_expo_unique(kk2) /= ao_expo_unique(kk) ) then + kk += 1 + ao_expo_unique(kk) = ao_expo_unique(kk2) + endif + kk2+=1 + enddo + ao_expo_unique_nucl_idx(1,k) = kkstart + ao_expo_unique_nucl_idx(2,k) = kk + kk += 1 + enddo + deallocate(order) + + ao_expo_unique_idx = 0 + do i=1,ao_num + do j=1,ao_prim_num(i) + do k=ao_expo_unique_nucl_idx(1,ao_nucl(i)),ao_expo_unique_nucl_idx(2,ao_nucl(i)) + if (ao_expo_unique(k) == ao_expo(i,j)) then + ao_expo_unique_idx(j,i) = k + exit + endif + enddo + enddo + ao_expo_unique_idx_1(i) = ao_expo_unique_idx(1,i) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_coef_transp, (ao_prim_num_max_8,ao_num) ] + implicit none + BEGIN_DOC +! Transposed of the ao_coef matrix + END_DOC + integer :: i,j + integer, save :: ifirst = 0 + if (ifirst==0) then + ifirst = 1 + ao_coef_transp = 0. + endif + call transpose(ao_coef,ao_num_8,ao_coef_transp,ao_prim_num_max_8,ao_num,ao_prim_num_max) +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_expo_transp, (ao_prim_num_max_8,ao_num) ] + implicit none + BEGIN_DOC +! Transposed of the ao_expo matrix + END_DOC + integer :: i,j + integer, save :: ifirst = 0 + if (ifirst==0) then + ifirst = 1 + ao_expo_transp = 0. + endif + call transpose(ao_expo,ao_num_8,ao_expo_transp,ao_prim_num_max_8,ao_num,ao_prim_num_max) +END_PROVIDER + diff --git a/src/AO/ao_axis.irp.f b/src/AO/ao_axis.irp.f new file mode 100644 index 0000000..5216ded --- /dev/null +++ b/src/AO/ao_axis.irp.f @@ -0,0 +1,129 @@ +subroutine pow_l(r,a,x1,x2,x3) + implicit none + BEGIN_DOC +! Fast calculation of powers for AO calculations + END_DOC + real, intent(in) :: r + integer, intent(in) :: a + real, intent(out) :: x1,x2,x3 + if (a==0) then + x1 = 1. + x2 = 0. + x3 = 0. + return + else if (a==1) then + x1 = r + x2 = 1. + x3 = 0. + return + else if (a==2) then + x1 = r*r + x2 = r + x3 = 1. + return + endif + + select case (a) + case (3) + x2 = r*r + x3 = r + x1 = x2*r + return + case (4) + x2 = r*r*r + x3 = r*r + x1 = x3*x3 + return + case (5:) + x3 = r**(a-2) + x2 = x3*r + x1 = x2*r + return + case (:-1) + x1 = 0. + x2 = 0. + x3 = 0. + return + end select +end function + + + BEGIN_PROVIDER [ real, ao_axis_block, ((-2*simd_sp+1):ao_block_num_8) ] +&BEGIN_PROVIDER [ real, ao_axis_grad_block_x, ((-2*simd_sp+1):ao_block_num_8) ] +&BEGIN_PROVIDER [ real, ao_axis_grad_block_y, ((-2*simd_sp+1):ao_block_num_8) ] +&BEGIN_PROVIDER [ real, ao_axis_grad_block_z, ((-2*simd_sp+1):ao_block_num_8) ] +&BEGIN_PROVIDER [ real, ao_axis_lapl_block, ((-2*simd_sp+1):ao_block_num_8) ] + implicit none + include '../types.F' + + BEGIN_DOC +! Cartesian polynomial part of the atomic orbitals. (blocked) + +! Gradients of the cartesian polynomial part of the atomic orbitals. (blocked) + +! Laplacian of the cartesian atomic orbitals (blocked) + END_DOC + + + integer :: i, j, k, l, idx + real, save :: real_of_int(-1:15) + data real_of_int /0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15./ + !DIR$ ATTRIBUTES ALIGN : 64 :: real_of_int + + if (ao_block_num_8 == 0) then + return + endif + + !DIR$ VECTOR ALIGNED + do idx=1,ao_oneD_prim_non_zero_idx(0) + i=ao_oneD_prim_non_zero_idx(idx) + + real :: p10, p11, p12 ! p10 contains x**(ao_power(i,1)) + real :: p20, p21, p22 ! p11 contains x**(ao_power(i,1)-1) + real :: p30, p31, p32 ! p12 contains x**(ao_power(i,1)-2) + real :: p012, p013, p023 + + integer :: inucl + inucl = ao_nucl(i) + + real :: x, y, z + x = nucl_elec_dist_vec(1,inucl,ao_elec) + y = nucl_elec_dist_vec(2,inucl,ao_elec) + z = nucl_elec_dist_vec(3,inucl,ao_elec) + + integer :: pow1, pow2, pow3 + + pow1 = ao_power_transp(1,i) + pow2 = ao_power_transp(2,i) + pow3 = ao_power_transp(3,i) + + !DIR$ FORCEINLINE + call pow_l(x,pow1,p10, p11, p12) + !DIR$ FORCEINLINE + call pow_l(y,pow2,p20, p21, p22) + !DIR$ FORCEINLINE + call pow_l(z,pow3,p30, p31, p32) + + p012 = real_of_int(pow3) * p10*p20 + p023 = p20*p30 + p013 = real_of_int(pow2) * p10*p30 + + ao_axis_block(idx) = p023 * p10 + p023 = real_of_int(pow1) * p023 + + ao_axis_lapl_block(idx) = real_of_int(pow1-1) * p023 * p12 & + + real_of_int(pow2-1) * p013 * p22 & + + real_of_int(pow3-1) * p012 * p32 + + ao_axis_grad_block_x(idx) = p023 * p11 + ao_axis_grad_block_y(idx) = p013 * p21 + ao_axis_grad_block_z(idx) = p012 * p31 + + + enddo + + +END_PROVIDER + + + diff --git a/src/AO/ao_full.irp.f b/src/AO/ao_full.irp.f new file mode 100644 index 0000000..75df2e5 --- /dev/null +++ b/src/AO/ao_full.irp.f @@ -0,0 +1,185 @@ + BEGIN_PROVIDER [ real, ao_value_block, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_grad_block_x, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_grad_block_y, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_grad_block_z, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_lapl_block, (ao_num_8) ] +&BEGIN_PROVIDER [ integer, ao_value_non_zero_idx, ((-simd_sp+1):ao_num_8) ] + + implicit none + BEGIN_DOC +! Values of the atomic orbitals (blocked) + +! Gradients of the atomic orbitals (blocked) + +! Laplacian of the atomic orbitals (blocked) + END_DOC + + if (ao_block_num == 0) then + return + endif + + integer :: j,idx,k, kmax + ao_value_non_zero_idx(0) = ao_oneD_prim_non_zero_idx(0) + kmax = ao_oneD_prim_non_zero_idx(0) + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do idx=1,kmax + ao_value_non_zero_idx(idx) = ao_oneD_prim_non_zero_idx(idx) + ao_value_block(idx) = ao_oneD_block(idx) * ao_axis_block(idx) + ao_grad_block_x(idx) = ao_oneD_block(idx) * ao_axis_grad_block_x(idx) + ao_oneD_grad_block_x(idx) * ao_axis_block(idx) + ao_grad_block_y(idx) = ao_oneD_block(idx) * ao_axis_grad_block_y(idx) + ao_oneD_grad_block_y(idx) * ao_axis_block(idx) + ao_grad_block_z(idx) = ao_oneD_block(idx) * ao_axis_grad_block_z(idx) + ao_oneD_grad_block_z(idx) * ao_axis_block(idx) + ao_lapl_block(idx) = ao_oneD_block(idx) * ao_axis_lapl_block(idx) + ao_oneD_lapl_block(idx) * ao_axis_block(idx) + & + 2.*(ao_oneD_grad_block_x(idx) * ao_axis_grad_block_x(idx) + & + ao_oneD_grad_block_y(idx) * ao_axis_grad_block_y(idx) + & + ao_oneD_grad_block_z(idx) * ao_axis_grad_block_z(idx) ) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ logical, reduce_primitives ] + implicit none + BEGIN_DOC +! If True, remove cusp AOs from AO basis + END_DOC + reduce_primitives = do_nucl_fitcusp +! if (calc_density_fit) then +! reduce_primitives = .False. +! endif +END_PROVIDER + + +BEGIN_PROVIDER [ logical, primitives_reduced ] + implicit none + BEGIN_DOC +! Tells if the number of primitives has been reduced due to the nucl_fitcusp + END_DOC + integer, save :: first_pass = 0 + + if (reduce_primitives.and.(first_pass == 0)) then + first_pass = 1 + integer :: i,j,k,l + integer :: prim_num_old, prim_num + prim_num_old = sum(ao_prim_num) + do k=1,nucl_num + point(1) = nucl_coord(k,1) + point(2) = nucl_coord(k,2) + point(3) = nucl_coord(k,3)+ nucl_fitcusp_radius(k) + TOUCH point + PROVIDE ao_oned_prim_p + PROVIDE ao_prim_num_max + PROVIDE ao_oned_p + PROVIDE ao_power + PROVIDE ao_coef + PROVIDE ao_nucl + do i=1,ao_num + if (ao_oned_p(i) /= 0.) then + l=ao_power(i,1)+ao_power(i,2)+ao_power(i,3) + if ( (l == 0).and.(ao_nucl(i) == k) ) then + do j=1,ao_prim_num_max + if (abs(ao_coef(i,j)*ao_oneD_prim_p(i,j)/ao_oneD_p(i)) < 1.e-4) then + ao_coef(i,j) = 0. + endif + enddo + endif + endif + enddo + enddo + + do i=1,ao_num + j=1 + do while (j1) then + ao_elec = j + TOUCH ao_elec + endif + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do idx=1,ao_value_non_zero_idx(0) + i=ao_value_non_zero_idx(idx) + ao_value(i,j) = ao_value_block(idx) + enddo + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (200) + do idx=1,ao_value_non_zero_idx(0) + i=ao_value_non_zero_idx(idx) + ao_grad_x(i,j) = ao_grad_block_x(idx) + ao_grad_y(i,j) = ao_grad_block_y(idx) + ao_grad_z(i,j) = ao_grad_block_z(idx) + ao_lapl(i,j) = ao_lapl_block(idx) + enddo + enddo + ao_elec = 1 + SOFT_TOUCH ao_elec + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_value_transp, (elec_num_8,ao_num) ] + implicit none + BEGIN_DOC +! Values of the atomic orbitals + END_DOC + integer :: i,j + do i=1,ao_num + !DIR$ VECTOR ALIGNED + do j=1,elec_num + ao_value_transp(j,i) = ao_value(i,j) + enddo + enddo +END_PROVIDER + diff --git a/src/AO/ao_oneD.irp.f b/src/AO/ao_oneD.irp.f new file mode 100644 index 0000000..679b2e7 --- /dev/null +++ b/src/AO/ao_oneD.irp.f @@ -0,0 +1,352 @@ +BEGIN_PROVIDER [ logical, compute_all_aos ] + implicit none + BEGIN_DOC +! If true, compute all AOs + END_DOC + + compute_all_aos = .not. do_nucl_fitcusp +END_PROVIDER + + +BEGIN_PROVIDER [ integer, ao_elec ] + implicit none + BEGIN_DOC +! Current electron in AO evaluation + END_DOC + ao_elec = 1 +END_PROVIDER + + + BEGIN_PROVIDER [ real, ao_radius, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_radius_sorted, (ao_num_8) ] +&BEGIN_PROVIDER [ real, nucl_radius, (nucl_num_8) ] +&BEGIN_PROVIDER [ integer, ao_radius_order, (ao_num_8) ] + implicit none + BEGIN_DOC +! Radius of the AOs +! +! Radius of the nuclei +! +! Indices of the AOs per nucleus ordered by ao_radius in descending order + END_DOC + integer :: i,k, inucl + real :: to_sort(ao_num) + !DIR$ VECTOR ALIGNED + nucl_radius = 0. + !DIR$ VECTOR ALIGNED + do i=1,ao_num + ao_radius_order(i) = i + enddo + do i=1,ao_num + ao_radius(i) = 0. + do k=1,ao_prim_num(i) + ao_radius(i) = max(20./ao_expo_transp(k,i), ao_radius(i)) + enddo + to_sort(i) = 1./ao_radius(i) + inucl = ao_nucl(i) + nucl_radius(inucl) = max(nucl_radius(inucl),ao_radius(i)) + enddo + !DIR$ VECTOR ALIGNED + nucl_radius = sqrt(nucl_radius) + integer :: istart, iend + do inucl=1,nucl_num + istart = ao_nucl_idx(1,inucl) + iend = ao_nucl_idx(2,inucl) + call isort(to_sort(istart),ao_radius_order(istart),(iend-istart+1)) + enddo + ao_radius_sorted = ao_radius + call set_order(ao_radius_sorted,ao_radius_order,ao_num) + +END_PROVIDER + + BEGIN_PROVIDER [ real, ao_oneD_block, (ao_num_8) ] +&BEGIN_PROVIDER [ integer, ao_oneD_prim_non_zero_idx, ((-simd_sp+1):ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_oneD_lapl_block, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_oneD_grad_block_z, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_oneD_grad_block_y, (ao_num_8) ] +&BEGIN_PROVIDER [ real, ao_oneD_grad_block_x, (ao_num_8) ] +&BEGIN_PROVIDER [ integer, ao_block_num ] +&BEGIN_PROVIDER [ integer, ao_block_num_8 ] + implicit none + BEGIN_DOC +! Gradients and Laplacians of the primitive AOs + END_DOC + + PROVIDE nucl_elec_dist + PROVIDE ao_power + PROVIDE ao_expo_unique + PROVIDE ao_expo_transp + PROVIDE ao_coef_transp + PROVIDE compute_all_aos + PROVIDE ao_prim_num + PROVIDE ao_expo + PROVIDE nucl_fitcusp_radius + + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + !DIR$ VECTOR ALIGNED + ao_oneD_prim_non_zero_idx = 0. + ao_oneD_block = 0. + ao_oneD_lapl_block = 0. + ao_oneD_grad_block_x = 0. + ao_oneD_grad_block_y = 0. + ao_oneD_grad_block_z = 0. + endif + + real :: ao_oneD_prim_block + real :: buffer2 + + real :: b1, b2, b3, b4 + integer :: i, j, k, idx, l, m + real :: rtemp, r, r2 + + integer :: inucl, istart, iend, idx0, ii + real :: buffer(-simd_sp+1:max(ao_nucl_idx(2,nucl_num), ao_expo_unique_nucl_idx(2,nucl_num)) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer + + j = ao_elec + + real :: d + real, parameter :: rdm = 20. + idx = 1 + !DIR$ VECTOR ALIGNED + buffer = 0. + do inucl=1,nucl_num + r = nucl_elec_dist(inucl,j) + if (r > nucl_radius(inucl)) then + cycle + else + r2 = r*r + istart = ao_expo_unique_nucl_idx(1,inucl) + if (ao_expo_unique(istart)*r2 > rdm) then + cycle + endif + real :: dcomp + dcomp = rdm/r2 + b1 = -nucl_elec_dist_vec(1,inucl,j)-nucl_elec_dist_vec(1,inucl,j) + b2 = -nucl_elec_dist_vec(2,inucl,j)-nucl_elec_dist_vec(2,inucl,j) + b3 = -nucl_elec_dist_vec(3,inucl,j)-nucl_elec_dist_vec(3,inucl,j) + b4 = 4.*r2 + iend = ao_expo_unique_nucl_idx(2,inucl) + + integer :: imax, istep + imax=istart+1 + !DIR$ LOOP COUNT(20) + do while (imax<=iend .and. ao_expo_unique(imax) <= dcomp) + imax += 1 + enddo + imax = imax-1 + + buffer(istart) = exp(-r2*ao_expo_unique(istart)) + if (imax - istart == 0) then + continue + else if (imax - istart == 1) then + buffer(istart+1) = exp(-r2*ao_expo_unique(istart+1)) + else if (imax - istart == 2) then + buffer(istart+1) = exp(-r2*ao_expo_unique(istart+1)) + buffer(istart+2) = exp(-r2*ao_expo_unique(istart+2)) + else + !DIR$ VECTOR UNALIGNED + !DIR$ LOOP COUNT MIN(3) + do i=istart+1,imax + buffer(i) = exp(-r2*ao_expo_unique(i)) + enddo + endif + istart = ao_nucl_idx(1,inucl) + iend = ao_nucl_idx(2,inucl) + idx0 = idx + + if ( r > nucl_fitcusp_radius(inucl) ) then + + !DIR$ VECTOR UNALIGNED + do ii=istart,iend + if ( r2 > ao_radius_sorted(ii)) then + exit + endif + ao_oneD_prim_non_zero_idx(idx) = ao_radius_order(ii) + idx += 1 + enddo + + else if (.not.compute_all_aos) then + + !DIR$ VECTOR UNALIGNED + do ii=istart,iend + if ( r2 > ao_radius_sorted(ii)) then + exit + endif + i = ao_radius_order(ii) + if ( .not.ao_power_is_zero(i) ) then + ao_oneD_prim_non_zero_idx(idx) = i + idx += 1 + endif + enddo + + else + + !DIR$ VECTOR UNALIGNED + do ii=istart,iend + if ( r2 > ao_radius_sorted(ii)) then + exit + endif + ao_oneD_prim_non_zero_idx(idx) = ao_radius_order(ii) + idx += 1 + enddo + + endif + + + iend = idx-1 + + real :: f, g, h + integer :: lmax + i = ao_oneD_prim_non_zero_idx(idx0) + if (i==0) then + cycle + endif + lmax = ao_prim_num(i) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0) = b1*f + ao_oneD_grad_block_y(idx0) = b2*f + ao_oneD_grad_block_z(idx0) = b3*f + + select case (iend-idx0) + case (0) + continue + case (1) + i = ao_oneD_prim_non_zero_idx(idx0+1) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0+1) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0+1) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0+1) = b1*f + ao_oneD_grad_block_y(idx0+1) = b2*f + ao_oneD_grad_block_z(idx0+1) = b3*f + lmax = max(ao_prim_num(i),lmax) + case (2) + i = ao_oneD_prim_non_zero_idx(idx0+1) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0+1) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0+1) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0+1) = b1*f + ao_oneD_grad_block_y(idx0+1) = b2*f + ao_oneD_grad_block_z(idx0+1) = b3*f + lmax = max(ao_prim_num(i),lmax) + + i = ao_oneD_prim_non_zero_idx(idx0+2) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0+2) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0+2) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0+2) = b1*f + ao_oneD_grad_block_y(idx0+2) = b2*f + ao_oneD_grad_block_z(idx0+2) = b3*f + lmax = max(ao_prim_num(i),lmax) + case (3) + i = ao_oneD_prim_non_zero_idx(idx0+1) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0+1) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0+1) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0+1) = b1*f + ao_oneD_grad_block_y(idx0+1) = b2*f + ao_oneD_grad_block_z(idx0+1) = b3*f + lmax = max(ao_prim_num(i),lmax) + + i = ao_oneD_prim_non_zero_idx(idx0+2) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0+2) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0+2) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0+2) = b1*f + ao_oneD_grad_block_y(idx0+2) = b2*f + ao_oneD_grad_block_z(idx0+2) = b3*f + lmax = max(ao_prim_num(i),lmax) + + i = ao_oneD_prim_non_zero_idx(idx0+3) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx0+3) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx0+3) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx0+3) = b1*f + ao_oneD_grad_block_y(idx0+3) = b2*f + ao_oneD_grad_block_z(idx0+3) = b3*f + lmax = max(ao_prim_num(i),lmax) + + case default + !DIR$ VECTOR UNALIGNED + !DIR$ LOOP COUNT MIN(4) + do idx=idx0+1,iend + i = ao_oneD_prim_non_zero_idx(idx) + g = buffer(ao_expo_unique_idx_1(i)) + f = ao_coef(i,1)*g*ao_expo(i,1) + ao_oneD_block(idx) = ao_coef(i,1)*g + ao_oneD_lapl_block(idx) = f*(b4 * ao_expo(i,1)-6.) + ao_oneD_grad_block_x(idx) = b1*f + ao_oneD_grad_block_y(idx) = b2*f + ao_oneD_grad_block_z(idx) = b3*f + lmax = max(ao_prim_num(i),lmax) + enddo + + end select + + if (lmax == 1) then + idx = iend+1 + cycle + endif + + do idx=idx0,iend + i = ao_oneD_prim_non_zero_idx(idx) + if (ao_prim_num(i) == 1) then + cycle + endif + l = ao_prim_num(i) + + ao_oneD_block(idx) = 0. + ao_oneD_lapl_block(idx) = 0. + f = 0. + g = 0. + h = 0. + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (8) + do k=1,l + ao_oneD_prim_block = buffer(ao_expo_unique_idx(k,i))*ao_coef_transp(k,i) + g = g + ao_oneD_prim_block + buffer2 = ao_expo_transp(k,i)*ao_oneD_prim_block + f = f+buffer2 + h = h + buffer2 * (b4*ao_expo_transp(k,i)-6.) + enddo + ao_oneD_block(idx) = g + ao_oneD_lapl_block(idx) = h + ao_oneD_grad_block_x(idx) = b1*f + ao_oneD_grad_block_y(idx) = b2*f + ao_oneD_grad_block_z(idx) = b3*f + enddo + idx = iend+1 + endif + enddo + + ao_oneD_prim_non_zero_idx(0) = idx-1 + if (ao_oneD_prim_non_zero_idx(0) > ao_block_num) then + ao_block_num = ao_oneD_prim_non_zero_idx(0) + integer :: mod_align + ao_block_num_8 = mod_align(ao_block_num) + endif + + if (ao_oneD_prim_non_zero_idx(0) == 0) then + ao_oneD_prim_non_zero_idx(0) = 1 + ao_oneD_prim_non_zero_idx(1) = 1 + ao_oneD_block(idx) = 0. + ao_oneD_lapl_block(idx) = 0. + ao_oneD_grad_block_x(idx) = 0. + ao_oneD_grad_block_y(idx) = 0. + ao_oneD_grad_block_z(idx) = 0. + endif + +END_PROVIDER + + diff --git a/src/AO/ao_point.irp.f b/src/AO/ao_point.irp.f new file mode 100644 index 0000000..4e2b511 --- /dev/null +++ b/src/AO/ao_point.irp.f @@ -0,0 +1,289 @@ +BEGIN_PROVIDER [ real, ao_axis_power_p, (-2:ao_power_max,3,nucl_num) ] + implicit none + BEGIN_DOC +! Evaluation of power of x, y, z at the current point for each +! nucleus. Negative power -> 0. + END_DOC + + integer :: i,k,l + do i=1,nucl_num + do l=1,3 + ao_axis_power_p(-2,l,i) = 0. + ao_axis_power_p(-1,l,i) = 0. + ao_axis_power_p(0,l,i) = 0. + ao_axis_power_p(0,l,i) = 1. + do k=1,ao_power_max_nucl(i,l) + ao_axis_power_p(k,l,i) = point_nucl_dist_vec(i,l)*ao_axis_power_p(k-1,l,i) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_axis_p, (ao_num) ] + implicit none + include '../types.F' + BEGIN_DOC +! Cartesian polynomial part of the atomic orbitals. + END_DOC + integer :: i + + do i=1,ao_num + ao_axis_p(i) & + = ao_axis_power_p( ao_power_transp(1,i) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i) , 3 , ao_nucl(i) ) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_axis_grad_p, (ao_num,3) ] + implicit none + include '../types.F' + BEGIN_DOC +! Gradients of the cartesian polynomial part of the atomic orbitals. + END_DOC + integer :: i, l + real :: real_of_int(-1:10) + data real_of_int /0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./ + + do i=1,ao_num + ao_axis_grad_p(i,1) = real_of_int(ao_power_transp(1,i)) & + * ao_axis_power_p( ao_power_transp(1,i)-1, 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i) , 3 , ao_nucl(i) ) + enddo + + do i=1,ao_num + ao_axis_grad_p(i,2) = real_of_int(ao_power_transp(2,i)) & + * ao_axis_power_p( ao_power_transp(1,i) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i)-1, 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i) , 3 , ao_nucl(i) ) + enddo + + do i=1,ao_num + ao_axis_grad_p(i,3) = real_of_int(ao_power_transp(3,i)) & + * ao_axis_power_p( ao_power_transp(1,i) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i)-1, 3 , ao_nucl(i) ) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_axis_lapl_p, (ao_num) ] + implicit none + include '../types.F' + BEGIN_DOC +! Laplacian of the cartesian atomic orbitals + END_DOC + integer :: i, j, l + + do i=1,ao_num + real :: real_of_int(-2:10) + data real_of_int /0.,0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./ + + ao_axis_lapl_p(i) & + = real_of_int(ao_power_transp(1,i)) & + * real_of_int(ao_power_transp(1,i)-1) & + * ao_axis_power_p( ao_power_transp(1,i)-2, 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i) , 3 , ao_nucl(i) ) & + + real_of_int(ao_power_transp(2,i)) & + * real_of_int(ao_power_transp(2,i)-1) & + * ao_axis_power_p( ao_power_transp(1,i) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i)-2, 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i) , 3 , ao_nucl(i) ) & + + real_of_int(ao_power_transp(3,i)) & + * real_of_int(ao_power_transp(3,i)-1) & + * ao_axis_power_p( ao_power_transp(1,i) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(2,i) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power_transp(3,i)-2, 3 , ao_nucl(i) ) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_value_p, (ao_num) ] + implicit none + BEGIN_DOC +! Values of the atomic orbitals + END_DOC + integer :: i + do i=1,ao_num + ao_value_p(i) = ao_oneD_p(i) * ao_axis_p(i) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_grad_p, (ao_num,3) ] + implicit none + include '../types.F' + BEGIN_DOC +! Gradients of the atomic orbitals + END_DOC + + integer :: i,l + do l=1,3 + do i=1,ao_num + ao_grad_p(i,l) = ao_oneD_p(i) * ao_axis_grad_p(i,l) + enddo + do i=1,ao_num + ao_grad_p(i,l) = ao_grad_p(i,l) + ao_oneD_grad_p(i,l) * ao_axis_p(i) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_lapl_p, (ao_num) ] + implicit none + include '../types.F' + BEGIN_DOC +! Laplacian of the atomic orbitals + END_DOC + integer :: i,l + do i=1,ao_num + ao_lapl_p(i) = ao_oneD_p(i) * ao_axis_lapl_p(i) + enddo + do i=1,ao_num + ao_lapl_p(i) = ao_lapl_p(i) + ao_oneD_lapl_p(i) * ao_axis_p(i) + enddo + do l=1,3 + do i=1,ao_num + ao_lapl_p(i) = ao_lapl_p(i) + 2.*ao_oneD_grad_p(i,l) * ao_axis_grad_p(i,l) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_num,ao_prim_num_max) ] + implicit none + include '../types.F' + + BEGIN_DOC +! Exponentials of the primitive AOs + END_DOC + integer :: i, k + real :: r2, rtemp + + ! Compute alpha*r or alpha*r^2 + do i=1,ao_num + r2 = point_nucl_dist(ao_nucl(i))*point_nucl_dist(ao_nucl(i)) + do k=1,ao_prim_num_max + ao_oneD_prim_p(i,k) = r2 + enddo + enddo + + ! Compute exp(-alpha*r) or exp(-alpha*r^2) + do i=1,ao_num + do k=1,ao_prim_num(i) + ao_oneD_prim_p(i,k) = exp(-ao_oneD_prim_p(i,k)*ao_expo(i,k)) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_p, (ao_num) ] + implicit none + include '../types.F' + + BEGIN_DOC +! One-dimensional component of the AOs + END_DOC + + integer :: i, k + + do i=1,ao_num + ao_oneD_p(i) = 0. + enddo + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_p(i) = ao_oneD_p(i) + ao_coef(i,k)*ao_oneD_prim_p(i,k) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_num,ao_prim_num_max,3) ] + implicit none + include '../types.F' + BEGIN_DOC +! Gradients of the one-dimensional component of the primitive AOs + END_DOC + integer :: i, k, l + do l=1,3 + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_prim_grad_p(i,k,l) = -2.*point_nucl_dist_vec(ao_nucl(i),l)*ao_expo(i,k)*ao_oneD_prim_p(i,k) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_grad_p, (ao_num,3) ] + implicit none + include '../types.F' + BEGIN_DOC +! Gradients of the one-dimensional component of the AOs + END_DOC + integer :: i, k, l + do l=1,3 + do i=1,ao_num + ao_oneD_grad_p(i,l) = 0. + enddo + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef(i,k)*ao_oneD_prim_grad_p(i,k,l) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_num,ao_prim_num_max) ] + implicit none + include '../types.F' + BEGIN_DOC +! Laplacian of the one-dimensional component of the primitive AOs + END_DOC + integer :: i, k + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_prim_lapl_p(i,k) = ao_oneD_prim_p(i,k) * ao_expo(i,k) *& + ( 4.*ao_expo(i,k)*point_nucl_dist(ao_nucl(i))**2 - 6. ) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_lapl_p, (ao_num) ] + implicit none + include '../types.F' + BEGIN_DOC +! Laplacian of the one-dimensional component of the AOs + END_DOC + + integer :: i, k + + do i=1,ao_num + ao_oneD_lapl_p(i) = 0. + enddo + do k=1,ao_prim_num_max + do i=1,ao_num + ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef(i,k)*ao_oneD_prim_lapl_p(i,k) + enddo + enddo + +END_PROVIDER + +