Acceleration

This commit is contained in:
Anthony Scemama 2009-06-15 11:33:16 +02:00
parent dc47a48108
commit da570c030f
6 changed files with 76 additions and 91 deletions

View File

@ -1,16 +1,17 @@
# MPI-ifort
IRPF90 = irpf90 -DMPI #-a -d
FC = mpif90 -static-intel -static-libgcc -xT
FCFLAGS= -O3 -g
FCFLAGS= -O3
# Gfortran
#FC = gfortran -ffree-line-length-none -static-libgcc -static
#FCFLAGS= -O3 -ffast-math -L ~/QCIO/lib
# Mono
#IRPF90 = irpf90
#FC = ifort -static-intel -static-libgcc
#FCFLAGS= -O3 -axP
IRPF90 = irpf90
FC = ifort -static-intel -static-libgcc
FCFLAGS= -O3 -axP
FCFLAGS= -O3 -xT -g
SRC=
OBJ=

View File

@ -1,38 +1,16 @@
double precision function Boys(x,n)
recursive double precision function Boys(x,n) result(res)
implicit none
include 'constants.F'
real, intent(in) :: x
integer, intent(in) :: n
integer :: k
real, parameter :: thr = 6.
integer ,parameter :: Nmax = 20
double precision :: fact,fact2
ASSERT (x > 0.)
if (n == 0) then
if (x > thr) then
Boys = 0.5d0*sqrt(pi/x)
else
Boys = 1.d0/dble(2*n+1)
do k=1,Nmax
Boys = Boys + (-x)**k/dble(fact(k)*(2*n+2*k+1))
enddo
endif
res = sqrt(pi/(4.*x))*erf(sqrt(x))
else
if (x > thr) then
Boys = fact2(2*n-1)*0.5d0**(n+1)*sqrt(pi/x**(2*n+1))
else
Boys = 1.d0/dble(2*n+1)
do k=1,Nmax
Boys = Boys + (-x)**k/dble(fact(k)*(2*n+2*k+1))
enddo
endif
res = (dble(2*n-1) * Boys(x,(n-1)) - exp(-x) )/(2.*x)
endif
end function
@ -96,3 +74,42 @@ double precision function fact(n)
end function
double precision function rintgauss(n)
implicit none
include 'constants.F'
integer :: n
rintgauss = sqrt(pi)
if ( n == 0 ) then
return
else if ( n == 1 ) then
rintgauss = 0.
else if ( mod(n,2) == 1) then
rintgauss = 0.
else
double precision :: fact2
rintgauss = rintgauss/(2.**(n/2))
rintgauss = rintgauss * fact2(n-1)
endif
end function
double precision function goverlap(gamA,gamB,nA)
implicit none
real :: gamA, gamB
integer :: nA(3)
double precision :: gamtot
gamtot = gamA+gamB
goverlap=1.0
integer :: l
double precision :: rintgauss
do l=1,3
goverlap = goverlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
enddo
end function

View File

@ -123,10 +123,10 @@ END_PROVIDER
!$OMP END CRITICAL (qcio_critical)
double precision :: norm, norm2
double precision :: overlap
double precision :: goverlap
do i=1,ao_num
do j=1,ao_prim_num(i)
norm = overlap(ao_expo(j,i),ao_expo(j,i),ao_power(i,:))
norm = goverlap(ao_expo(j,i),ao_expo(j,i),ao_power(i,:))
norm = sqrt(norm)
ao_coef(j,i) = buffer(j,i)/norm
enddo
@ -134,43 +134,3 @@ END_PROVIDER
END_PROVIDER
double precision function rintgauss(n)
implicit none
integer :: n
double precision :: pi
pi = acos(-1.)
rintgauss = sqrt(pi)
if ( n == 0 ) then
return
else if ( n == 1 ) then
rintgauss = 0.
else if ( mod(n,2) == 1) then
rintgauss = 0.
else
double precision :: fact2
rintgauss = rintgauss/(2.**(n/2))
rintgauss = rintgauss * fact2(n-1)
endif
end function
double precision function overlap(gamA,gamB,nA)
implicit none
real :: gamA, gamB
integer :: nA(3)
double precision :: gamtot
gamtot = gamA+gamB
overlap=1.0
integer :: l
double precision :: rintgauss
do l=1,3
overlap = overlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
enddo
end function

View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [ real, eplf_gamma ]
BEGIN_DOC
! Value of the gaussian for the EPLF
END_DOC
eplf_gamma = 10.
eplf_gamma = 10000.
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
@ -27,27 +27,33 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
END_DOC
integer :: i, j, k, l
double precision :: t
PROVIDE ao_eplf_integral_matrix
PROVIDE mo_coef
do i=1,mo_num
do j=i,mo_num
mo_eplf_integral_matrix(j,i) = 0.
enddo
enddo
do i=1,mo_num
do k=1,ao_num
if (mo_coef(k,i) /= 0.) then
do j=i,mo_num
do l=1,ao_num
mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(l,k)
enddo
enddo
endif
do k=1,ao_num
if (abs(mo_coef(k,i)) /= 0.) then
do l=1,ao_num
t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
do j=i,mo_num
mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
t*mo_coef_transp(j,l)
enddo
endif
enddo
endif
enddo
do j=i,mo_num
enddo
do i=1,mo_num
do j=i+1,mo_num
mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
enddo
enddo
@ -132,7 +138,7 @@ double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,
real, intent(in) :: a,b,gmma ! Exponents
real, intent(in) :: xa,xb,xr ! Centers
integer, intent(in) :: i,j ! Powers of xa and xb
integer,parameter :: Npoints=1000
integer,parameter :: Npoints=10000
real :: x, xmin, xmax, dx
ASSERT (a>0.)
@ -165,12 +171,13 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
integer :: p,q,k
double precision :: integral
double precision :: ao_eplf_integral_primitive_oneD_numeric
real :: gmma, center(3)
real :: gmma, center(3), c
ao_eplf_integral_numeric = 0.
do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i)
c = ao_coef(p,i)*ao_coef(q,j)
integral = &
ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(p,i), &
@ -199,7 +206,7 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
ao_power(j,3), &
gmma, &
center(3))
ao_eplf_integral_numeric = ao_eplf_integral_numeric + integral*ao_coef(p,i)*ao_coef(q,j)
ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral
enddo
enddo

View File

@ -27,7 +27,7 @@ double precision function primitive_overlap_oneD_numeric(a,xa,i,b,xb,j)
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
integer,parameter :: Npoints=10000
real :: x, xmin, xmax, dx
ASSERT (a>0.)

View File

@ -1,7 +1,7 @@
program debug
PROVIDE ao_prim_num_max
read(*,*) eplf_gamma
TOUCH eplf_gamma
!read(*,*) eplf_gamma
!TOUCH eplf_gamma
call run()
end