LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-12-21 11:53:30 +01:00
Code Releases Activity

Added AO

Browse Source
This commit is contained in:
Anthony Scemama 2015-12-19 03:32:49 +01:00
parent ea185b405f
commit 9580050b80
5 changed files with 1320 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

365
src/AO/ao.irp.f Normal file
View File

@ -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

129
src/AO/ao_axis.irp.f Normal file
View File

@ -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

185
src/AO/ao_full.irp.f Normal file
View File

@ -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 (j<ao_prim_num(i))
if (ao_coef(i,j) == 0.) then
do k=j,ao_prim_num(i)-1
ao_coef(i,k) = ao_coef(i,k+1)
ao_expo(i,k) = ao_expo(i,k+1)
enddo
ao_expo(i,ao_prim_num(i)) = 0.
ao_coef(i,ao_prim_num(i)) = 0.
ao_prim_num(i) = ao_prim_num(i)-1
else
j = j+1
endif
enddo
enddo
TOUCH ao_expo ao_coef ao_prim_num
prim_num = sum(ao_prim_num)
character*(80) :: message
write(message,'(A8,I4,A)') 'Removed ',prim_num_old-prim_num,' primitive gaussians'
call info(irp_here,message)
call iinfo(irp_here,"sum_ao_prim_num",sum(ao_prim_num))
primitives_reduced = .True.
FREE ao_oneD_p ao_oneD_prim_p
else
primitives_reduced = .False.
endif
END_PROVIDER
BEGIN_PROVIDER [ real, ao_value, (ao_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, ao_grad_x , (ao_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, ao_grad_y , (ao_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, ao_grad_z , (ao_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, ao_lapl , (ao_num_8,elec_num) ]
implicit none
BEGIN_DOC
! Values, Gradients and Laplacians of the atomic orbitals
END_DOC
integer :: i,j,idx,k
!DIR$ VECTOR ALIGNED
ao_value = 0.
!DIR$ VECTOR ALIGNED
ao_grad_x = 0.
!DIR$ VECTOR ALIGNED
ao_grad_y = 0.
!DIR$ VECTOR ALIGNED
ao_grad_z = 0.
!DIR$ VECTOR ALIGNED
ao_lapl = 0.
do j=1,elec_num
if (j>1) 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

352
src/AO/ao_oneD.irp.f Normal file
View File

@ -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

289
src/AO/ao_point.irp.f Normal file
View File

@ -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