mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-10-18 05:51:30 +02:00
354 lines
9.1 KiB
Fortran
354 lines
9.1 KiB
Fortran
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_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.
|
|
! AO coefficients are such that the AOs are normalized, and the primitives
|
|
! are also normalized
|
|
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
|
|
|