LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-18 11:15:30 +02:00
Code Releases Activity

new files

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-03-13 16:18:06 +01:00
parent 3059085d57
commit ca96c59429
7 changed files with 560 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

135
src/MCQC/AOtoMO_oooooo.f90 Normal file
View File

@ -0,0 +1,135 @@
subroutine AOtoMO_oooooo(nBas,nO,cO,O,oooOooo)
! AO to MO transformation of three-electron integrals for the block oooooo
! Semi-direct O(N^7) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: scr1(:,:,:,:,:,:),scr2(:,:,:,:,:,:)
integer :: mu,nu,la,si,ka,ta,i,j,k,l,m,n
! Output variables
double precision,intent(out) :: oooOooo(nO,nO,nO,nO,nO,nO)
! Memory allocation
allocate(scr1(nBas,nBas,nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas,nBas,nBas))
write(*,*)
write(*,'(A42)') '------------------------------------------'
write(*,'(A42)') ' AO to MO transformation for oooooo block '
write(*,'(A42)') '------------------------------------------'
write(*,*)
scr1 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do ta=1,nBas
do n=1,nO
scr1(mu,nu,la,si,ka,n) = scr1(mu,nu,la,si,ka,n) + O(mu,nu,la,si,ka,ta)*cO(ta,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do i=1,nO
do n=1,nO
scr2(i,nu,la,si,ka,n) = scr2(i,nu,la,si,ka,n) + cO(mu,i)*scr1(mu,nu,la,si,ka,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do i=1,nO
do m=1,nO
do n=1,nO
scr1(i,nu,la,si,m,n) = scr1(i,nu,la,si,m,n) + scr2(i,nu,la,si,m,n)*cO(ka,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do i=1,nO
do j=1,nO
do m=1,nO
do n=1,nO
scr2(i,j,la,si,m,n) = scr2(i,j,la,si,m,n) + cO(nu,j)*scr1(i,nu,la,si,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do la=1,nBas
do si=1,nBas
do i=1,nO
do j=1,nO
do l=1,nO
do m=1,nO
do n=1,nO
scr1(i,j,la,l,m,n) = scr1(i,j,la,l,m,n) + cO(si,l)*scr2(i,j,la,si,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
oooOooo = 0d0
do si=1,nBas
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
do m=1,nO
do n=1,nO
oooOooo(i,j,k,l,m,n) = oooOooo(i,j,k,l,m,n) + cO(la,k)*scr1(i,j,la,k,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(scr1,scr2)
end subroutine AOtoMO_oooooo

1
src/MCQC/MP2F12.f90
View File

@ -55,6 +55,7 @@ subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,EcMP2,c,EcMP2F12)
! Compute the three-electron part of the MP2-F12 energy
allocate(oooFCooo(nO,nO,nO,nO,nO,nO))
call AOtoMO_oooooo(nBas,nO,cO,FC,oooFCooo)
E3a = 0d0
E3b = 0d0

29
src/MCQC/NormCoeff.f90 Normal file
View File

@ -0,0 +1,29 @@
function NormCoeff(alpha,a)
implicit none
! Input variables
double precision,intent(in) :: alpha
integer,intent(in) :: a(3)
! local variable
double precision :: pi,dfa(3),dfac
integer :: atot
! Output variable
double precision NormCoeff
pi = 4d0*atan(1d0)
atot = a(1) + a(2) + a(3)
dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1)))
dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2)))
dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3)))
NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot
NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3))
NormCoeff = sqrt(NormCoeff)
end function NormCoeff

170
src/MCQC/elements.f90 Normal file
View File

@ -0,0 +1,170 @@
function element_number(element_name)
implicit none
integer,parameter :: nelement_max = 103
character(len=2),intent(in) :: element_name
integer :: element_number
character(len=2),parameter :: element_list(nelement_max) = &
(/' H', 'He', & ! 2
'Li','Be', ' B',' C',' N',' O',' F','Ne', & ! 10
'Na','Mg', 'Al','Si',' P',' S','Cl','Ar', & ! 18
' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', & ! 36
'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', & ! 54
'Cs','Ba', & ! 56
'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & ! 70
'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', & ! 86
'Fr','Ra', & ! 88
'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No', & ! 102
'Lr' & ! 103
/)
!=====
integer :: ielement
!=====
ielement=1
do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) )
if( ielement == nelement_max ) then
write(*,'(a,a)') ' Input symbol ',element_name
write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements'
write(*,*) '!!! element symbol not understood !!!'
stop
endif
ielement = ielement + 1
enddo
element_number = ielement
end function element_number
function element_core(zval,zatom)
implicit none
double precision,intent(in) :: zval
double precision,intent(in) :: zatom
integer :: element_core
!=====
!
! If zval /= zatom, this is certainly an effective core potential
! and no core states should be frozen.
if( ABS(zval - zatom) > 1d0-3 ) then
element_core = 0
else
if( zval <= 4.00001d0 ) then ! up to Be
element_core = 0
else if( zval <= 12.00001d0 ) then ! up to Mg
element_core = 1
else if( zval <= 30.00001d0 ) then ! up to Ca
element_core = 5
else if( zval <= 48.00001d0 ) then ! up to Sr
element_core = 9
else
write(*,*) '!!! not imlemented in element_core !!!'
stop
endif
endif
end function element_core
function element_covalent_radius(zatom)
! Return covalent radius of an atom
implicit none
include 'parameters.h'
integer,intent(in) :: zatom
double precision :: element_covalent_radius
!
! Data from Cambridge Structural Database
! http://en.wikipedia.org/wiki/Covalent_radius
!
! Values are first given in picometer
! They will be converted in bohr later on
select case(zatom)
case( 1)
element_covalent_radius = 31.
case( 2)
element_covalent_radius = 28.
case( 3)
element_covalent_radius = 128.
case( 4)
element_covalent_radius = 96.
case( 5)
element_covalent_radius = 84.
case( 6)
element_covalent_radius = 73.
case( 7)
element_covalent_radius = 71.
case( 8)
element_covalent_radius = 66.
case( 9)
element_covalent_radius = 57.
case(10) ! Ne.
element_covalent_radius = 58.
case(11)
element_covalent_radius = 166.
case(12)
element_covalent_radius = 141.
case(13)
element_covalent_radius = 121.
case(14)
element_covalent_radius = 111.
case(15)
element_covalent_radius = 107.
case(16)
element_covalent_radius = 105.
case(17)
element_covalent_radius = 102.
case(18) ! Ar.
element_covalent_radius = 106.
case(19)
element_covalent_radius = 203.
case(20)
element_covalent_radius = 176.
case(21)
element_covalent_radius = 170.
case(22)
element_covalent_radius = 160.
case(23)
element_covalent_radius = 153.
case(24)
element_covalent_radius = 139.
case(25)
element_covalent_radius = 145.
case(26)
element_covalent_radius = 145.
case(27)
element_covalent_radius = 140.
case(28)
element_covalent_radius = 124.
case(29)
element_covalent_radius = 132.
case(30)
element_covalent_radius = 122.
case(31)
element_covalent_radius = 120.
case(32)
element_covalent_radius = 119.
case(34)
element_covalent_radius = 120.
case(35)
element_covalent_radius = 120.
case(36) ! Kr.
element_covalent_radius = 116.
case default
write(*,*) '!!! covalent radius not available !!!'
stop
end select
! pm to bohr conversion
element_covalent_radius = element_covalent_radius*pmtoau
end function element_covalent_radius

118
src/MCQC/read_basis.f90 Normal file
View File

@ -0,0 +1,118 @@
subroutine read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
! Read basis set information
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nAt,nO
double precision,intent(in) :: rAt(nAt,ncart)
! Local variables
integer :: nShAt,iAt,iShell
integer :: i,j,k
character :: shelltype
! Output variables
integer,intent(out) :: nShell,nBas
double precision,intent(out) :: CenterShell(maxShell,ncart)
integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
integer,intent(out) :: nV
!------------------------------------------------------------------------
! Primary basis set information
!------------------------------------------------------------------------
! Open file with basis set specification
open(unit=2,file='input/basis')
! Read basis information
write(*,'(A28)') 'Gaussian basis set'
write(*,'(A28)') '------------------'
nShell = 0
do i=1,nAt
read(2,*) iAt,nShAt
write(*,'(A28,1X,I16)') 'Atom n. ',iAt
write(*,'(A28,1X,I16)') 'number of shells ',nShAt
write(*,'(A28)') '------------------'
! Basis function centers
do j=1,nShAt
nShell = nShell + 1
do k=1,ncart
CenterShell(nShell,k) = rAt(iAt,k)
enddo
! Shell type and contraction degree
read(2,*) shelltype,KShell(nShell)
if(shelltype == "S") then
TotAngMomShell(nShell) = 0
write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
elseif(shelltype == "P") then
TotAngMomShell(nShell) = 1
write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
elseif(shelltype == "D") then
TotAngMomShell(nShell) = 2
write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
elseif(shelltype == "F") then
TotAngMomShell(nShell) = 3
write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
elseif(shelltype == "G") then
TotAngMomShell(nShell) = 4
write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
elseif(shelltype == "H") then
TotAngMomShell(nShell) = 5
write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
elseif(shelltype == "I") then
TotAngMomShell(nShell) = 6
write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
endif
! Read exponents and contraction coefficients
write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
do k=1,Kshell(nShell)
read(2,*) ExpShell(nShell,k),DShell(nShell,k)
write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
enddo
enddo
write(*,'(A28)') '------------------'
enddo
! Total number of shells
write(*,'(A28,1X,I16)') 'Number of shells',nShell
write(*,'(A28)') '------------------'
write(*,*)
! Close file with basis set specification
close(unit=2)
! Calculate number of basis functions
nBas = 0
do iShell=1,nShell
nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2
enddo
write(*,'(A28)') '------------------'
write(*,'(A28,1X,I16)') 'Number of basis functions',NBas
write(*,'(A28)') '------------------'
write(*,*)
! Number of virtual orbitals
nV = nBas - nO
end subroutine read_basis

68
src/MCQC/read_geometry.f90 Normal file
View File

@ -0,0 +1,68 @@
subroutine read_geometry(nAt,ZNuc,rA,ENuc)
! Read molecular geometry
implicit none
include 'parameters.h'
! Ouput variables
integer,intent(in) :: nAt
! Local variables
integer :: i,j
double precision :: RAB
character(len=2) :: El
integer,external :: element_number
! Ouput variables
double precision,intent(out) :: ZNuc(NAt),rA(nAt,ncart),ENuc
! Open file with geometry specification
open(unit=1,file='input/molecule')
! Read geometry
read(1,*)
read(1,*)
read(1,*)
do i=1,nAt
read(1,*) El,rA(i,1),rA(i,2),rA(i,3)
ZNuc(i) = element_number(El)
enddo
! Compute nuclear repulsion energy
ENuc = 0
do i=1,nAt-1
do j=i+1,nAt
RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2
ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB)
enddo
enddo
! Close file with geometry specification
close(unit=1)
! Print geometry
write(*,'(A28)') '------------------'
write(*,'(A28)') 'Molecular geometry'
write(*,'(A28)') '------------------'
do i=1,NAt
write(*,'(A28,1X,I16)') 'Atom n. ',i
write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i)
write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,ncart)
enddo
write(*,*)
write(*,'(A28)') '------------------'
write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
write(*,'(A28)') '------------------'
write(*,*)
end subroutine read_geometry

39
src/MCQC/read_molecule.f90 Normal file
View File

@ -0,0 +1,39 @@
subroutine read_molecule(nAt,nEl,nO,nCore,nRyd)
! Read number of atoms nAt and number of electrons nEl
implicit none
include 'parameters.h'
! Input variables
integer,intent(out) :: nAt,nEl,nO,nCore,nRyd
! Open file with geometry specification
open(unit=1,file='input/molecule')
! Read number of atoms and number of electrons
read(1,*)
read(1,*) nAt,nEl,nCore,nRyd
nO = nEl/2
! Print results
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of atoms',nAt
write(*,'(A28)') '----------------------'
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of electrons',nEl
write(*,'(A28)') '----------------------'
write(*,*)
! Close file with geometry specification
close(unit=1)
end subroutine read_molecule