eplf/src/ao.irp.f

BEGIN_PROVIDER [ integer, ao_num ]
 implicit none

 BEGIN_DOC
! Number of atomic orbitals
 END_DOC

 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

END_PROVIDER

BEGIN_PROVIDER [ integer, ao_prim_num, (ao_num) ]
 implicit none

 BEGIN_DOC
! Number of primitives per atomic orbital
 END_DOC

 ao_prim_num = -1
 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

END_PROVIDER

BEGIN_PROVIDER [ integer, ao_nucl, (ao_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
 write(range,'(A,I5,A)') '(1,',nucl_num,')'
 integer :: i
 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
 enddo

END_PROVIDER

BEGIN_PROVIDER [ integer, ao_power, (ao_num,3) ]
 implicit none

 BEGIN_DOC
! x,y,z powers of the atomic orbital
 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
 enddo
END_PROVIDER


BEGIN_PROVIDER [ integer , ao_power_max ]
 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,3) ]
 implicit none
 BEGIN_DOC
! Maximum powers of x, y and z per nucleus
 END_DOC
 integer :: i, j
 do j=1,3
  do i=1,nucl_num
   ao_power_max_nucl(i,j) = 0
  enddo
 enddo

 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
END_PROVIDER


BEGIN_PROVIDER [ integer, ao_prim_num_max ]
 implicit none

 BEGIN_DOC
! Max Number of primitives per atomic orbital
 END_DOC

 ao_prim_num_max = maxval(ao_prim_num)

END_PROVIDER

 BEGIN_PROVIDER [ real, ao_expo_transp, (ao_num,ao_prim_num_max) ]
&BEGIN_PROVIDER [ real, ao_coef_transp, (ao_num,ao_prim_num_max) ]
 implicit none
 BEGIN_DOC
! Exponents and coefficients of the atomic orbitals
 END_DOC

 ao_expo_transp = 0.
 call get_ao_basis_ao_expo(ao_expo_transp)

 integer :: i,j
 do i=1,ao_num
  do j=1,ao_prim_num(i)
   if (ao_expo_transp(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_transp = 0.
 call get_ao_basis_ao_coef(ao_coef_transp)

! Normalization of the AO coefficients
! ------------------------------------
 double precision :: norm, norm2
 double precision :: goverlap
 integer :: pow(3), l
 do i=1,ao_num
  do j=1,ao_prim_num(i)
   pow(1) = ao_power(i,1)
   pow(2) = ao_power(i,2)
   pow(3) = ao_power(i,3)
   norm = goverlap(ao_expo_transp(i,j),ao_expo_transp(i,j),pow)
   ao_coef_transp(i,j) = ao_coef_transp(i,j)/sqrt(norm)
  enddo
 enddo

END_PROVIDER

 BEGIN_PROVIDER [ real, ao_expo, (ao_prim_num_max,ao_num) ]
&BEGIN_PROVIDER [ real, ao_coef, (ao_prim_num_max,ao_num) ]
 implicit none
 BEGIN_DOC
! Exponents and coefficients of the atomic orbitals
 END_DOC

 integer :: i,j
 do i=1,ao_num
  do j=1,ao_prim_num(i)
   ao_coef(j,i) = ao_coef_transp(i,j)
   ao_expo(j,i) = ao_expo_transp(i,j)
  enddo
 enddo

END_PROVIDER