From ca96c594294530015b1c6a5bd4ea218bb467eb17 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 13 Mar 2019 16:18:06 +0100 Subject: [PATCH] new files --- src/MCQC/AOtoMO_oooooo.f90 | 135 +++++++++++++++++++++++++++++ src/MCQC/MP2F12.f90 | 1 + src/MCQC/NormCoeff.f90 | 29 +++++++ src/MCQC/elements.f90 | 170 +++++++++++++++++++++++++++++++++++++ src/MCQC/read_basis.f90 | 118 +++++++++++++++++++++++++ src/MCQC/read_geometry.f90 | 68 +++++++++++++++ src/MCQC/read_molecule.f90 | 39 +++++++++ 7 files changed, 560 insertions(+) create mode 100644 src/MCQC/AOtoMO_oooooo.f90 create mode 100644 src/MCQC/NormCoeff.f90 create mode 100644 src/MCQC/elements.f90 create mode 100644 src/MCQC/read_basis.f90 create mode 100644 src/MCQC/read_geometry.f90 create mode 100644 src/MCQC/read_molecule.f90 diff --git a/src/MCQC/AOtoMO_oooooo.f90 b/src/MCQC/AOtoMO_oooooo.f90 new file mode 100644 index 0000000..10b9c0e --- /dev/null +++ b/src/MCQC/AOtoMO_oooooo.f90 @@ -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 diff --git a/src/MCQC/MP2F12.f90 b/src/MCQC/MP2F12.f90 index 5fda128..0ca0a1b 100644 --- a/src/MCQC/MP2F12.f90 +++ b/src/MCQC/MP2F12.f90 @@ -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 diff --git a/src/MCQC/NormCoeff.f90 b/src/MCQC/NormCoeff.f90 new file mode 100644 index 0000000..9e6cabf --- /dev/null +++ b/src/MCQC/NormCoeff.f90 @@ -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 diff --git a/src/MCQC/elements.f90 b/src/MCQC/elements.f90 new file mode 100644 index 0000000..22953dc --- /dev/null +++ b/src/MCQC/elements.f90 @@ -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 + diff --git a/src/MCQC/read_basis.f90 b/src/MCQC/read_basis.f90 new file mode 100644 index 0000000..2f3768c --- /dev/null +++ b/src/MCQC/read_basis.f90 @@ -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 diff --git a/src/MCQC/read_geometry.f90 b/src/MCQC/read_geometry.f90 new file mode 100644 index 0000000..60c60b8 --- /dev/null +++ b/src/MCQC/read_geometry.f90 @@ -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 diff --git a/src/MCQC/read_molecule.f90 b/src/MCQC/read_molecule.f90 new file mode 100644 index 0000000..0894e17 --- /dev/null +++ b/src/MCQC/read_molecule.f90 @@ -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