Numeric integrals for test. Analytic integrals OK

This commit is contained in:
Anthony Scemama 2009-05-14 15:25:18 +02:00
parent af8b3d3dfc
commit 1284eab958
3 changed files with 148 additions and 28 deletions

View File

@ -1,7 +1,6 @@
IRPF90 = irpf90 # -a #-d
FC = ifort
FCFLAGS= -O2 #-xT -fast-transcendentals
IRPF90 = irpf90 #-a -d
FC = ifort
FCFLAGS= -O3 -xT
SRC=
OBJ=

View File

@ -1,9 +1,42 @@
program debug
implicit none
double precision :: eplf_integral, ao_overlap
PROVIDE ao_prim_num_max
integer :: i,j
integer :: k
read(*,*) i,j
print *, eplf_integral(i,j,eplf_gamma,point)
print *, ao_overlap(i,j)
print *, ''
do k=1,nucl_num
print *, nucl_coord(k,:)
enddo
print *, ''
print *, 'AO ', i
print *, 'prim num:', ao_prim_num(i)
print *, 'powers :', ao_power(i,:)
print *, 'center :', ao_nucl(i)
print *, 'expo / coef'
do k=1,ao_prim_num(i)
print *, ao_expo(k,i), ao_coef(k,i)
enddo
print *, ''
print *, 'AO ', j
print *, 'prim num:', ao_prim_num(j)
print *, 'powers :', ao_power(j,:)
print *, 'center :', ao_nucl(j)
print *, 'expo / coef'
do k=1,ao_prim_num(j)
print *, ao_expo(k,j), ao_coef(k,j)
enddo
double precision :: ao_overlap, ao_overlap_numeric
print *, ''
print *, 'Overlap integral :', ao_overlap(i,j)
print *, 'Overlap integral N :', ao_overlap_numeric(i,j)
double precision :: eplf_integral, eplf_integral_numeric
print *, ''
print *, 'EPLF gamma : ', eplf_gamma
print *, 'EPLF integral :', eplf_integral(i,j,eplf_gamma,point)
print *, 'EPLF integral N :', eplf_integral_numeric(i,j,eplf_gamma,point)
end

View File

@ -1,3 +1,112 @@
BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Overlap matrix between the Atomic Orbitals
END_DOC
integer :: i, j
double precision :: ao_overlap
do j=1,ao_num
do i=j,ao_num
ao_overlap_matrix(i,j) = ao_overlap(i,j)
enddo
enddo
do j=1,ao_num
do i=1,j-1
ao_overlap_matrix(i,j) = ao_overlap(j,i)
enddo
enddo
END_PROVIDER
double precision function primitive_overlap_oneD_numeric(a,xa,i,b,xb,j)
implicit none
include 'constants.F'
real, intent(in) :: a,b ! Exponents
real, intent(in) :: xa,xb ! Centers
integer, intent(in) :: i,j ! Powers of xa and xb
integer,parameter :: Npoints=1000
real :: x, xmin, xmax, dx
ASSERT (a>0.)
ASSERT (b>0.)
ASSERT (i>=0)
ASSERT (j>=0)
xmin = min(xa,xb) - 10.
xmax = max(xa,xb) + 10.
dx = (xmax-xmin)/real(Npoints)
real :: dtemp
dtemp = 0.
x = xmin
integer :: k
do k=1,Npoints
dtemp = dtemp + &
(x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2))
x = x+dx
enddo
primitive_overlap_oneD_numeric = dtemp*dx
end function
double precision function ao_overlap_numeric(i,j)
implicit none
integer, intent(in) :: i, j
integer :: p,q,k
double precision :: integral(ao_prim_num_max,ao_prim_num_max)
double precision :: primitive_overlap_oneD_numeric
ASSERT(i>0)
ASSERT(j>0)
ASSERT(i<=ao_num)
ASSERT(j<=ao_num)
do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i)
integral(p,q) = &
primitive_overlap_oneD_numeric ( &
ao_expo(p,i), &
nucl_coord(ao_nucl(i),1), &
ao_power(i,1), &
ao_expo(q,j), &
nucl_coord(ao_nucl(j),1), &
ao_power(j,1) ) * &
primitive_overlap_oneD_numeric ( &
ao_expo(p,i), &
nucl_coord(ao_nucl(i),2), &
ao_power(i,2), &
ao_expo(q,j), &
nucl_coord(ao_nucl(j),2), &
ao_power(j,2) ) * &
primitive_overlap_oneD_numeric ( &
ao_expo(p,i), &
nucl_coord(ao_nucl(i),3), &
ao_power(i,3), &
ao_expo(q,j), &
nucl_coord(ao_nucl(j),3), &
ao_power(j,3) )
enddo
enddo
do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i)
integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
enddo
enddo
ao_overlap_numeric = 0.
do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i)
ao_overlap_numeric = ao_overlap_numeric + integral(p,q)
enddo
enddo
end function
subroutine gaussian_product(a,xa,b,xb,k,p,xp)
implicit none
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
@ -139,24 +248,3 @@ double precision function ao_overlap(i,j)
end function
BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Overlap matrix between the Atomic Orbitals
END_DOC
integer :: i, j
double precision :: ao_overlap
do j=1,ao_num
do i=j,ao_num
ao_overlap_matrix(i,j) = ao_overlap(i,j)
enddo
enddo
do j=1,ao_num
do i=1,j-1
ao_overlap_matrix(i,j) = ao_overlap(j,i)
enddo
enddo
END_PROVIDER