LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-06-02 11:25:18 +02:00
Code Releases Activity
4981bbe9e2
qmcchem/src/AO/ao.irp.f

352 lines
9.0 KiB
Fortran
Raw Blame History

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
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
View Git Blame Copy Permalink