program integral c **************** Program INTEGRAL **************** c This program is designed to calculate one- and c two-electron AO integrals and to write them out to c disk in canonical order (in Dirac <12|12> c convention). It is designed to handle only S and P c orbitals. c c Copyright (C) 1992 Jeff Nichols and Jack Simons c c ************************************************** implicit real*8 (a-h,o-z) include "limits.h" parameter (d2=2.d0,d1=1.d0,d34=0.75d0, & d12=0.5d0,d14=0.25d0,d0=0.d0) logical yesno character*30 dfname, fname, dscrpt character*1 ans dimension ixyz(0:nht,((nht+1)*(nht+2))/2,3) dimension al(nsh,nprim),co(nsh,nprim) dimension npr(nsh),nshk(nsh),xsh(3,nsh) dimension nhi(4),alp(4),x(3,4),s(ntlp1**4) dimension iptr(nsh), ssum(ntlp1**4) dimension znuc(natoms),xnuc(3,natoms) dimension ovlap(nprim,nprim) dimension cv(nprim),scr1(nprim) dimension sones(mxo2),tones(mxo2), & vones(mxo2),tvones(mxo2) dimension twos(mxo4) isym2(i,j)=(max0(i,j)*(max0(i,j)-1))/2 + min0(i,j) isym4(i,j,k,l)=(max0(isym2(i,j),isym2(k,l))* & (max0(isym2(i,j),isym2(k,l))-1))/2 + & min0(isym2(i,j),isym2(k,l)) c c calculate memory needed in reals (hence *2 for doubles below) c itotal = (nht+1*((nht+1)*(nht+2))/2*3) + &nsh*nprim*2 + nsh*nprim*2 + &nsh + nsh + 3*nsh*2 + &4 + 4*2 + 3*4*2 + ntlp1**4*2 + &nsh + ntlp1**4*2 + &natoms*2 + 3*natoms*2 + &nprim*nprim*2 + &nprim*2 + nprim*2 + &mxo2*2 + mxo2*2 + &mxo2*2 + mxo2*2 + &mxo4*2 + &natoms + natoms + & 38*2 + 3*2 + 3*3*(NHT+1)*(NHT+1)*(NHT*2+1)*2 +3*5*5*5*2 call modinfo('integral ') call shwmem('integral ',itotal) call goon c c Initialize ixyz. c do 103 i=0,nht ik = 0 do 102 ix=i,0,-1 do 101 iy=i-ix,0,-1 ik = ik + 1 ixyz(i,ik,1) = ix ixyz(i,ik,2) = iy ixyz(i,ik,3) = i - ix - iy 101 continue 102 continue 103 continue c c zero all arrays: c call dzero(al,nsh*nprim) call dzero(co,nsh*nprim) call izero(npr,nsh) call izero(nshk,nsh) call dzero(xsh,3*nsh) call izero(nhi,4) call dzero(alp,4) call dzero(x,12) call dzero(s,ntlp1**4) call izero(iptr,nsh) call dzero(ssum,ntlp1**4) call dzero(znuc,natoms) call dzero(xnuc,3*natoms) call dzero(ovlap,nprim*nprim) call dzero(cv,nprim) call dzero(scr1,nprim) call dzero(sones,mxo2) call dzero(tones,mxo2) call dzero(vones,mxo2) call dzero(tvones,mxo2) call dzero(twos,mxo4) 99 continue write(*,*)' There is a library file which contains gaussian ' write(*,*)' atomic orbital basis sets for Hydrogen - Neon. ' write(*,*)' ' write(*,*)' The basis sets available to choose from are: ' write(*,*)' ' write(*,*)' 1.) STO3G by Hehre, Stewart, and Pople ' write(*,*)' JCP, 51, 2657 (1969) ' write(*,*)' 2.) 3-21G by Brinkley, Pople, and Hehre ' write(*,*)' JACS, 102, 939 (1980) ' write(*,*)' 3.) [3s2p] by Dunning and Hay in: ' write(*,*)' Modern Theoretical Chemistry Vol 3. ' write(*,*)' Henry F. Schaefer III, Editor ' write(*,*)' 1977, Plenum Press, NY ' write(*,*)' 4.) Input your own search string. ' write(*,*)' ' write(*,*)' The alternative option is to input your own ' write(*,*)' basis. ' write(*,*) write(*,*)' Input your choice: < 1, 2, or 3 > ' write(*,*) write(*,*)' 1.) Use basis set library, ' write(*,*)' 2.) Input your own, ' write(*,*)' 3.) Quit. ' write(*,*) read(*,*)igo if(igo.eq.1)then call gbasl(al,co,npr,nshk,xsh,znuc,xnuc,nat,numsh,norbs) elseif(igo.eq.2)then call guinp(al,co,npr,nshk,xsh,znuc,xnuc,nat,numsh,norbs) elseif(igo.eq.3)then call qmexit else write(*,*)' Unrecognized option. Retry? < y/n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 99 else call qmexit endif endif pi = 3.1415926535897932385d0 pippi = (d12/pi)**d34 c c Normalize the contraction coefficients c do 9 i = 1, numsh ni = nshk(i) if(ni.eq.0)then ovx = d12 + 1.0d0 else ovx = d12 + 2.0d0 endif cvx = d12 * ovx do 2 ip = 1, npr(i) cv(ip) = co(i,ip) ai = al(i,ip) do 1 jp = 1, npr(i) aj = al(i,jp) ovlap(ip,jp) = (d2*sqrt(ai*aj)/(ai+aj))**ovx 1 continue 2 continue do 3 ip = 1, npr(i) scr1(ip) = ovlap(ip,1) * cv(1) 3 continue do 5 jp = 2, npr(i) do 4 ip = 1, npr(i) scr1(ip) = scr1(ip) + ovlap(ip,jp) * cv(jp) 4 continue 5 continue sum = 00d0 do 6 ip = 1, npr(i) sum = sum + cv(ip)*scr1(ip) 6 continue sum = 1.0d0/sqrt(sum) do 7 ip = 1, npr(i) cv(ip) = sum*cv(ip) 7 continue do 8 ip = 1, npr(i) all=4.0d00*al(i,ip) cop=cv(ip)*pippi*all**cvx co(i,ip)=cop 8 continue 9 continue iptr(1) = 0 do 10 i = 1, numsh - 1 iptr(i+1) = iptr(i) + & ((nshk(i)+1)*(nshk(i)+2))/2 10 continue c c Calculate nuclear repulsion c call crpnc(xnuc,znuc,repnuc,nat) write(*,1030)repnuc C C Open the Nuclear Repulsion file. C dscrpt = 'nuclear repulsion' ids = 17 dfname = 'NucRep' idf = 6 fname = ' ' ifn = 30 iflnum = 31 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then c c Write nuclear repulsion to disk. c write(31,1020)repnuc c c Close nuclear repulsion file. c close(unit=31) endif write(*,*)' Begin 1e- integral evaluation. ' c c Do the 1e- integrals. c do 18 ish = 1, numsh write(*,*)' Begin evaluating shell: ',ish x(1,1) = xsh(1,ish) x(2,1) = xsh(2,ish) x(3,1) = xsh(3,ish) nhi(1) = nshk(ish) nki = ((nhi(1)+1)*(nhi(1)+2))/2 do 17 jsh = 1, ish x(1,2) = xsh(1,jsh) x(2,2) = xsh(2,jsh) x(3,2) = xsh(3,jsh) nhi(2) = nshk(jsh) nkj = ((nhi(2)+1)*(nhi(2)+2))/2 lns = nki*nkj do 11 i = 1, 3*lns ssum(i) = 0d0 11 continue do 16 ip=1,npr(ish) as = al(ish,ip) xc1 = co(ish,ip) do 15 jp=1,npr(jsh) bs = al(jsh,jp) xc2 = xc1*co(jsh,jp) call mfint1(ixyz,s,lns,as,bs,x(1,1),x(1,2), & nhi(1),nhi(2),xnuc,znuc,nat) do 12 i = 1, 3*lns ssum(i) = ssum(i) + s(i)*xc2 12 continue ij = 0 do 14 i = 1, nki ir = i+iptr(ish) do 13 j = 1, nkj ij = ij + 1 jr = j+iptr(jsh) irjr = isym2(ir,jr) sones(irjr) = ssum(ij) tones(irjr) = ssum(ij + lns) vones(irjr) = ssum(ij + 2*lns) tvones(irjr) = ssum(ij + lns) + & ssum(ij + 2*lns) 13 continue 14 continue 15 continue 16 continue 17 continue 18 continue write(*,*)' End 1e- integral evaluation. ' write(*,*)' Begin 2e- integral evaluation. ' c c Do the 2e- integrals. c do 32 ish = 1,numsh write(*,*)' Begin evaluating shell: ',ish x(1,1) = xsh(1,ish) x(2,1) = xsh(2,ish) x(3,1) = xsh(3,ish) nhi(1) = nshk(ish) nki = ((nhi(1)+1)*(nhi(1)+2))/2 do 31 jsh = 1,ish x(1,2) = xsh(1,jsh) x(2,2) = xsh(2,jsh) x(3,2) = xsh(3,jsh) nhi(2) = nshk(jsh) nkj = ((nhi(2)+1)*(nhi(2)+2))/2 do 30 ksh = 1,ish x(1,3) = xsh(1,ksh) x(2,3) = xsh(2,ksh) x(3,3) = xsh(3,ksh) nhi(3) = nshk(ksh) nkk = ((nhi(3)+1)*(nhi(3)+2))/2 lshend = ksh if(ish.eq.ksh)lshend = jsh do 29 lsh = 1,lshend x(1,4) = xsh(1,lsh) x(2,4) = xsh(2,lsh) x(3,4) = xsh(3,lsh) nhi(4) = nshk(lsh) nkl = ((nhi(4)+1)*(nhi(4)+2))/2 lns = nki*nkj*nkk*nkl do 19 i=1,lns ssum(i) = 0d0 19 continue do 24 ip=1,npr(ish) alp(1) = al(ish,ip) xc1 = co(ish,ip) do 23 jp=1,npr(jsh) alp(2) = al(jsh,jp) xc2 = xc1*co(jsh,jp) do 22 kp=1,npr(ksh) alp(3) = al(ksh,kp) xc3 = xc2*co(ksh,kp) do 21 lp=1,npr(lsh) alp(4) = al(lsh,lp) xc4 = xc3*co(lsh,lp) call mfint2(ixyz,s,lns,alp,x,nhi) do 20 i = 1, lns ssum(i) = ssum(i) + s(i)*xc4 20 continue 21 continue 22 continue 23 continue 24 continue ijkl = 0 do 28 i = 1, nki ir = i + iptr(ish) do 27 j = 1, nkj jr = j + iptr(jsh) do 26 k = 1, nkk kr = k + iptr(ksh) do 25 l = 1, nkl lr = l + iptr(lsh) ijkl = ijkl + 1 twos(isym4(ir,jr,kr,lr)) = ssum(ijkl) 25 continue 26 continue 27 continue 28 continue 29 continue 30 continue 31 continue 32 continue write(*,*)' End 2e- integral evaluation. ' c c Open the integral file. c dscrpt = 'overlap integrals' ids = 17 dfname = 'overlap' idf = 7 fname = ' ' ifn = 30 iflnum = 31 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then c c Write integrals to disk. c c Write the one-electron integrals. c do 34 i = 1, norbs do 33 j = 1,norbs write(31,1020)sones(isym2(i,j)) 33 continue 34 continue c c Close overlap integral file. c close(unit=31) endif dscrpt = '1e- (T+V) and 2e- integrals' ids = 27 dfname = 'ao_integrals.dat' idf = 16 fname = ' ' ifn = 30 iflnum = 32 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call ptintd(tvones,twos,norbs) c c Close 1e- and 2e- integral file. c close(unit=32) endif call ptintu(tvones,twos,norbs) call qmexit stop 1020 format(1x,f15.8) 1030 format(/,' The nuclear repulsion is: ',f15.8,/) end subroutine crpnc(xnuc,znuc,repnuc,nat) implicit real*8 (a-h,o-z) dimension znuc(nat),xnuc(3,nat) repnuc = 0.0e0 do 4 i = 2, nat do 2 j = 1, i - 1 rab = sqrt( (xnuc(1,i) - xnuc(1,j))**2 + & (xnuc(2,i) - xnuc(2,j))**2 + & (xnuc(3,i) - xnuc(3,j))**2 ) repnuc = repnuc + znuc(i)*znuc(j)/rab 2 continue 4 continue return end subroutine gbasl(al,co,npr,nshk,xsh,znuc,xnuc,nat,numsh,norbs) implicit real*8 (a-h,o-z) include "limits.h" character*80 title(natoms) character*1 ans,itype(2) character*8 ifind(30),igoget,istrng dimension al(nsh,nprim),co(nsh,nprim) dimension npr(nsh),nshk(nsh),xsh(3,nsh) dimension znuc(natoms),xnuc(3,natoms) dimension ibas(natoms),jshll(natoms) data ifind /'H_STO3G ','H_3-21G ','H_DH32 ', & 'He_STO3G','He_3-21G','He_DH32 ', & 'Li_STO3G','Li_3-21G','Li_DH32 ', & 'Be_STO3G','Be_3-21G','Be_DH32 ', & 'B_STO3G ','B_3-21G ','B_DH32 ', & 'C_STO3G ','C_3-21G ','C_DH32 ', & 'N_STO3G ','N_3-21G ','N_DH32 ', & 'O_STO3G ','O_3-21G ','O_DH32 ', & 'F_STO3G ','F_3-21G ','F_DH32 ', & 'Ne_STO3G','Ne_3-21G','Ne_DH32 '/ data itype /'S','P'/ c c Open basis set library. c open(unit=29,file='BasisLib',access='sequential', & status='old',form = 'formatted') 97 continue nat = 0 numsh = 0 norbs = 0 write(*,*)' How many atoms are in your system? ' write(*,*) read(*,*)nat if(nat.gt.natoms)then write(*,*)' Too many atoms. See programs limits. ' call qmexit endif do 4 i = 1, nat write(*,1000)i read(*,*)znuc(i) write(*,1010)i read(*,*)xnuc(1,i),xnuc(2,i),xnuc(3,i) c xnuc(1,i) = xnuc(1,i)/0.5291771 c xnuc(2,i) = xnuc(2,i)/0.5291771 c xnuc(3,i) = xnuc(3,i)/0.5291771 98 continue write(*,1020)i read(*,*)ibas(i) c c Get bases. c if(ibas(i).lt.4)then ii = 3*(nint(znuc(i)) - 1) + ibas(i) igoget = ifind(ii) else write(*,*)' Input search string. < A8 > ' read(*,'(a8)')igoget write(*,*)' Search string is: ',igoget endif rewind(unit=29) 99 continue read(29,'(a8)',end=101)istrng if(istrng.eq.igoget)then read(29,'(a80)')title(i) read(29,*)mshlls ncntr = 0 do 3 m = 1, mshlls ncntr = ncntr + 1 numsh = numsh + 1 if(numsh.gt.nsh)then write(*,*)' Too many shells. See programs limits. ' call qmexit endif do 1 l = 1, 3 xsh(l,numsh) = xnuc(l,i) 1 continue read(29,*)npr(numsh),nshk(numsh) if(nshk(numsh).eq.0)then norbs = norbs + 1 elseif(nshk(numsh).eq.1)then norbs = norbs + 3 endif if(norbs.gt.maxorb)then write(*,*)' Too many orbitals. See programs limits. ' call qmexit endif do 2 ip = 1, npr(numsh) read(29,*)al(numsh,ip),co(numsh,ip) 2 continue 3 continue goto 102 else goto 99 endif 101 continue write(*,*)' This basis is not available for ' write(*,*)' this atom. Retry? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 98 else call qmexit endif 102 continue jshll(i) = ncntr 4 continue jc = 0 do 7 i = 1, nat write(*,1030)i,znuc(i) write(*,1040)i,(xnuc(j,i),j=1,3) write(*,*)' Would you like to view the basis set ' write(*,*)' for this atom? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,1050)i do 6 j = 1, jshll(i) jc = jc + 1 write(*,1060)j,itype(nshk(jc)+1) write(*,1070) do 5 ip = 1, npr(jc) write(*,1080)ip,al(jc,ip),co(jc,ip) 5 continue 6 continue endif 7 continue write(*,1090)norbs write(*,*)' Is this information correct? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' Begin integral evaluation. ' else write(*,*)' Retry? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 97 else write(*,*)' Exiting program. ' call qmexit endif endif c c Close basis set library file. c close(unit=29) return 1000 format(' Input the atomic number of atom ',i1) 1010 format(' Input the coordinates of atom ',i1,' in ',/, & ' atomic units ') 1020 format(' Input the basis set choice for atom ',i1,/, & ' 1.) STO3G, or ',/, & ' 2.) 3-21G, or ',/, & ' 3.) D. & H. [3s2p], ',/, & ' 4.) Search String. ',//, & ' Input choice: < 1, 2, 3, or 4 > ') 1030 format(//,' The atomic number of atom ',i1,' is: ',f4.1) 1040 format(' The coordinates of atom ',i1,' are: ',/, & 3(3x,f15.8)) 1050 format(' The basis set choice for atom ',i1,' is: ') 1060 format(/,' Shell : ',i2,' is a ',a1,' function. ') 1070 format(6x,' Primitives ',8x,' Contraction Coefficients ') 1080 format(1x,i1,2x,f15.8,8x,f15.8) 1090 format(/,' Total number of orbitals is: ',i2,/) end subroutine guinp(al,co,npr,nshk,xsh,znuc,xnuc,nat,numsh,norbs) implicit real*8 (a-h,o-z) include "limits.h" character*1 ans,itype(2) dimension al(nsh,nprim),co(nsh,nprim) dimension npr(nsh),nshk(nsh),xsh(3,nsh) dimension znuc(natoms),xnuc(3,natoms) dimension jshll(natoms) data itype /'S','P'/ 97 continue nat = 0 numsh = 0 norbs = 0 write(*,*)' How many atoms are in your system? ' write(*,*) read(*,*)nat if(nat.gt.natoms)then write(*,*)' Too many atoms. See programs limits. ' call qmexit endif do 4 i = 1, nat write(*,1000)i read(*,*)znuc(i) write(*,1010)i read(*,*)xnuc(1,i),xnuc(2,i),xnuc(3,i) 98 continue c c Get bases. c 99 continue write(*,*)' Input the TOTAL number of contracted ' write(*,*)' functions (shells) on this atom. ' write(*,*)' e.g. 3S2P = 5 ' read(*,*)mshlls ncntr = 0 do 3 m = 1, mshlls ncntr = ncntr + 1 numsh = numsh + 1 if(numsh.gt.nsh)then write(*,*)' Too many shells. See programs limits. ' call qmexit endif do 1 l = 1, 3 xsh(l,numsh) = xnuc(l,i) 1 continue 101 continue write(*,1020)m read(*,'(a)')ans if(ans.eq.'S'.or.ans.eq.'s')then nshk(numsh) = 0 elseif(ans.eq.'P'.or.ans.eq.'p')then nshk(numsh) = 1 else write(*,*)' Undefined response. ' goto 101 endif if(nshk(numsh).eq.0)then norbs = norbs + 1 elseif(nshk(numsh).eq.1)then norbs = norbs + 3 endif if(norbs.gt.maxorb)then write(*,*)' Too many orbitals. See programs limits. ' call qmexit endif write(*,1030)m read(*,*)npr(numsh) if(npr(numsh).gt.nprim)then write(*,*)' Too many primitives. See programs limits. ' call qmexit endif do 2 ip = 1, npr(numsh) 102 continue write(*,1040)ip read(*,*)al(numsh,ip),co(numsh,ip) write(*,1050)ip,al(numsh,ip),co(numsh,ip) read(*,'(a)')ans if(ans.eq.'N'.or.ans.eq.'n')goto 102 2 continue 3 continue jshll(i) = ncntr 4 continue jc = 0 do 7 i = 1, nat write(*,1060)i,znuc(i) write(*,1070)i,(xnuc(j,i),j=1,3) write(*,*)' Would you like to view the basis set ' write(*,*)' for this atom? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,1080)i do 6 j = 1, jshll(i) jc = jc + 1 write(*,1090)j,itype(nshk(jc)+1) write(*,1100) do 5 ip = 1, npr(jc) write(*,1110)ip,al(jc,ip),co(jc,ip) 5 continue 6 continue endif 7 continue write(*,1120)norbs write(*,*)' Is this information correct? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' Begin integral evaluation. ' return else write(*,*)' Retry? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 97 else write(*,*)' Exiting program. ' call qmexit endif endif return 1000 format(' Input the atomic number of atom ',i1) 1010 format(' Input the coordinates of atom ',i1,' in ',/, & ' atomic units ') 1020 format(' Is shell: ',i2,' an S or P shell < S or P > ') 1030 format(' How many primitives in shell: ',i2) 1040 format(' Input primitive exponent and contraction', & ' coefficient ',i1) 1050 format(' Primitive exponent and contraction', & ' coefficient ',i1,' is: ',/, & 1x,2(2x,f15.8),' Is this correct? < y or n > ') 1060 format(//,' The atomic number of atom ',i1,' is: ',f4.1) 1070 format(' The coordinates of atom ',i1,' are: ',/, & 3(3x,f15.8)) 1080 format(' The basis set for atom ',i1,' is: ') 1090 format(/,' Shell : ',i2,' is an ',a1,' function. ') 1100 format(6x,' Primitives ',8x,' Contraction Coefficients ') 1110 format(1x,i1,2x,f15.8,8x,f15.8) 1120 format(/,' Total number of orbitals is: ',i2,/) end C################mfint1################################################# SUBROUTINE MFINT1(IXYZ,S,LDS,A,B,CA,CB,NA,NB,X,Z,NAT) C---------------------------------------------------------------------+ C variables: | C S - FINAL MATRIX OF INTEGRALS | C 1 - (A| |B) | C 2 - (A|T|B) | C 3 - (A|V|B) | C LDS - LEADING DIMENSION FOR S | C (NA+1)*(NA+2)/2 * (NB+1)*(NB+2)/2 | C A - EXPONENT ON FIRST FUNCTION | C B - EXPONENT ON SECOND FUNCTION | C CA - ARRAY OF COORDINATES FOR A IN ATOMIC UNTIS | C CB - ARRAY OF COORDINATES FOR B IN ATOMIC UNTIS | C NA - PRINCIPLE QUAMTUM NUMBER FOR A (ie s=0, p=1 ..) | C NB - PRINCIPLE QUAMTUM NUMBER FOR B | C X - ARRAY OF NUCLEAR COORDINATES | C Z - ARRAY OF NUCLEAR CHARGES | C NAT - NUMBER OF NUCLEI | C NHT - MAXIMUM QUANTUM NUMBER | C IXYZ - CARTESIAN NUMBERS FOR X^l Y^m Z^n | C | C formulae: | C | C l1,l2 m1,m2 n1,n2 | C S(A,B) = E1 * E2 * E3 | C 0 0 0 <2.38>| C | C 1 / l1,l2 m1,m2 n1,n2 \ | C T(A,B) = - |T * E2 * E3 + ... | | C 2 \ X 0 0 / <2.54>| C | C l,p l-1,p-1 l-1,p+1 l+1,p-1 l+1,p+1 | C T = lpE1 - 2lbE1 - 2apE1 + 4abE1 | C X 0 0 0 0 <2.55>| C | C -- l1,l2 m1,m2 n1,n2 | C V(A,B) = > E1 * E2 * E3 * R(L,M,N) | C --LMN L M N | C | C mwfeyereisen 1992 | C---------------------------------------------------------------------+ IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) PARAMETER (NHT=5) DIMENSION S(LDS,3),CA(3),CB(3),X(3,NAT),Z(NAT),P(3),PA(3),PB(3), * PQ(3),E(3,0:NHT,0:NHT,0:NHT*2),R(0:4,0:4,0:4), * RSUM(0:4,0:4,0:4) DIMENSION IXYZ(0:NHT,((NHT+1)*(NHT+2))/2,3) P5 = 0.5D0 PI = 3.1415926535897932385D0 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C MAKE E +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ABI = P5/(A+B) DO 10 I=1,3 P (I) = 2D0*ABI*(A*CA(I)+B*CB(I)) PA(I) = P(I) - CA(I) 10 PB(I) = P(I) - CB(I) SAB = (2D0*PI*ABI)**1.5D0*EXP(MAX(-80D0,-2D0*A*B*ABI* * ((CA(1)-CB(1))**2+(CA(2)-CB(2))**2+(CA(3)-CB(3))**2))) CALL MAKEE(E,PA,PB,ABI,SAB,NA+1,NB+1,NHT) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C MAKE R FOR NUCLEAR REPULSION INTEGRALS++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ DO 20 I=0,4 DO 20 J=0,4 DO 20 K=0,4 R (I,J,K) = 0D0 20 RSUM(I,J,K) = 0D0 W = A+B SR = 2D0/SQRT(PI)*SQRT(W) C loop over nuclear centers & add R into RSUM DO 40 N = 1,NAT PQ(1) = P(1)-X(1,N) PQ(2) = P(2)-X(2,N) PQ(3) = P(3)-X(3,N) T = W*(PQ(1)**2+PQ(2)**2+PQ(3)**2) CALL MAKER(R,SR,T,W,PQ,NA+NB) DO 30 I=0,4 DO 30 J=0,4 DO 30 K=0,4 30 RSUM(I,J,K) = RSUM(I,J,K) + R(I,J,K)*Z(N) 40 CONTINUE C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C LOOP OVER CARTESIAN COMPONENTS AND COMPLETE INTEGRALS+++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ L = 0 DO 500 NI = 1,((NA+1)*(NA+2))/2 L1 = IXYZ(NA,NI,1) M1 = IXYZ(NA,NI,2) N1 = IXYZ(NA,NI,3) DO 500 NJ = 1,((NB+1)*(NB+2))/2 L2 = IXYZ(NB,NJ,1) M2 = IXYZ(NB,NJ,2) N2 = IXYZ(NB,NJ,3) L = L + 1 C.............................. C COMPUTE OVERLAP INTEGRALS. C.............................. S(L,1) = E(1,L1,L2,0)*E(2,M1,M2,0)*E(3,N1,N2,0) C............................... C COMPUTE KINETIC INTEGRALS. C............................... TX = 2D0*A*B*E(1,L1+1,L2+1,0) TY = 2D0*A*B*E(2,M1+1,M2+1,0) TZ = 2D0*A*B*E(3,N1+1,N2+1,0) IF (L1.GT.0) TX=TX-B*L1*E(1,L1-1,L2+1,0) IF (M1.GT.0) TY=TY-B*M1*E(2,M1-1,M2+1,0) IF (N1.GT.0) TZ=TZ-B*N1*E(3,N1-1,N2+1,0) IF (L2.GT.0) TX=TX-A*L2*E(1,L1+1,L2-1,0) IF (M2.GT.0) TY=TY-A*M2*E(2,M1+1,M2-1,0) IF (N2.GT.0) TZ=TZ-A*N2*E(3,N1+1,N2-1,0) IF (L1.GT.0.AND.L2.GT.0) TX=TX+P5*L1*L2*E(1,L1-1,L2-1,0) IF (M1.GT.0.AND.M2.GT.0) TY=TY+P5*M1*M2*E(2,M1-1,M2-1,0) IF (N1.GT.0.AND.N2.GT.0) TZ=TZ+P5*N1*N2*E(3,N1-1,N2-1,0) S(L,2) = TX *E(2,M1,M2,0)*E(3,N1,N2,0) * + E(1,L1,L2,0)*TY *E(3,N1,N2,0) * + E(1,L1,L2,0)*E(2,M1,M2,0)*TZ C......................................... C COMPUTE NUCLEAR REPULSION INTEGRALS. C......................................... S(L,3) = 0 DO 100 NT=0,N1+N2 DO 100 MT=0,M1+M2 DO 100 LT=0,L1+L2 S(L,3) = S(L,3) - E(1,L1,L2,LT)*E(2,M1,M2,MT) * * E(3,N1,N2,NT)*RSUM(LT,MT,NT) 100 CONTINUE 500 CONTINUE RETURN END C################mfint2################################################# SUBROUTINE MFINT2(IXYZ,S,LNS,ALP,X,NHI) IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) PARAMETER (NHT=5) DIMENSION S(LNS),ALP(4),X(3,4),NHI(4),P(3),Q(3),PQ(3),PA(3), * PB(3),QC(3),QD(3),R(0:4,0:4,0:4), * E1(3,0:NHT,0:NHT,0:NHT*2),E2(3,0:NHT,0:NHT,0:NHT*2) DIMENSION IXYZ(0:NHT,((NHT+1)*(NHT+2))/2,3) C---------------------------------------------------------------------+ C variables: | C S - FINAL MATRIX OF INTEGRALS | C ALP - ARRAY OF EXPONENTS | C X - ARRAY OF COORDINATES FOR A IN ATOMIC UNTIS | C NHI - ARRAY OF PRINCIPLE QUAMTUM NUMBERS (ie s=0, p=1 ..) | C NHT - MAXIMUM QUANTUM NUMBER | C IXYZ - CARTESIAN NUMBERS FOR X^l Y^m Z^n | C | C formulae: | C 1/2 5/4 __2 | C 2 PI / -abAB \ | C KAB = --------- EXP| ----- | | C a+b \ a + b / | C | C p*q 2 / 1 \1/2 | C W = --- , T = W*PQ , SR = KAB*KCD*|---| | C p+q \P+Q/ | C | C -- l1l2m1m2n1n2 l3l4m3m4n3n4 | C I = > E1 *E2 *R(L+L',M+M',N+N') | C --LMNL'M'N' LMN L'M'N' | C | C mwfeyereisen 1992 | C---------------------------------------------------------------------+ SRF = 139.94734662099890277426D0 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C MAKE E1 & E2++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ NA = NHI(1) NB = NHI(2) NC = NHI(3) ND = NHI(4) A = ALP(1) B = ALP(2) C = ALP(3) D = ALP(4) ABI = 0.5D0/(A+B) CDI = 0.5D0/(C+D) DO 10 I=1,3 P (I) = 2D0*ABI*(A*X(I,1)+B*X(I,2)) Q (I) = 2D0*CDI*(C*X(I,3)+D*X(I,4)) PA(I) = P(I) - X(I,1) PB(I) = P(I) - X(I,2) QC(I) = Q(I) - X(I,3) QD(I) = Q(I) - X(I,4) 10 PQ(I) = P(I) - Q(I) CALL MAKEE(E1,PA,PB,ABI,1D0,NA,NB,NHT) CALL MAKEE(E2,QC,QD,CDI,1D0,NC,ND,NHT) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C MAKE R++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SR = SRF*ABI*CDI/SQRT(A+B+C+D)*EXP(MAX(-80D0,-2D0*(A*B*ABI*((X(1,1 * )-X(1,2))**2+(X(2,1)-X(2,2))**2+(X(3,1)-X(3,2))**2)+C*D*CDI*( * (X(1,3)-X(1,4))**2+(X(2,3)-X(2,4))**2+(X(3,3)-X(3,4))**2)))) W = 0.5D0/(ABI+CDI) T = W*(PQ(1)**2+PQ(2)**2+PQ(3)**2) CALL MAKER(R,SR,T,W,PQ,NA+NB+NC+ND) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C LOOP OVER CARTESIAN COMPONENTS AND COMPLETE INTEGRALS+++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ L = 0 DO 30 NI = 1,((NA+1)*(NA+2))/2 L1 = IXYZ(NA,NI,1) M1 = IXYZ(NA,NI,2) N1 = IXYZ(NA,NI,3) DO 30 NJ = 1,((NB+1)*(NB+2))/2 L2 = IXYZ(NB,NJ,1) M2 = IXYZ(NB,NJ,2) N2 = IXYZ(NB,NJ,3) DO 30 NK = 1,((NC+1)*(NC+2))/2 L3 = IXYZ(NC,NK,1) M3 = IXYZ(NC,NK,2) N3 = IXYZ(NC,NK,3) DO 30 NL = 1,((ND+1)*(ND+2))/2 L4 = IXYZ(ND,NL,1) M4 = IXYZ(ND,NL,2) N4 = IXYZ(ND,NL,3) L = L + 1 S(L) = 0 DO 20 NP=0,N3+N4 DO 20 MP=0,M3+M4 DO 20 LP=0,L3+L4 E2XYZ = E2(1,L3,L4,LP)*E2(2,M3,M4,MP)*E2(3,N3,N4,NP) * * (-1D0)**(LP+MP+NP) DO 20 NT=0,N1+N2 DO 20 MT=0,M1+M2 DO 20 LT=0,L1+L2 S(L) = S(L)+R(LT+LP,MT+MP,NT+NP)*E2XYZ * * E1(1,L1,L2,LT)*E1(2,M1,M2,MT)*E1(3,N1,N2,NT) 20 CONTINUE 30 CONTINUE RETURN END C################makee################################################## SUBROUTINE MAKEE(E,PA,PB,ABI,SAB,NA,NB,NE) IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) DIMENSION E(3,0:NE,0:NE,0:NE+NE),PA(3),PB(3) C---------------------------------------------------------------------+ C | C Ex(0,0,0) = SAB | C | C E(i+1,j,N) = ABI*E(i,j,N-1)+PA*E(i,j,N)+(N+1)*E(i,j,N+1) | C | C E(i,j+1,N) = ABI*E(i,j,N-1)+PB*E(i,j,N)+(N+1)*E(i,j,N+1) | C | C---------------------------------------------------------------------+ DO 10 I=1,3 DO 10 J=0,NE DO 10 K=0,NE DO 10 L=0,2*NE 10 E(I,J,K,L) = 0D0 C.....COMPUTE E(0,0,0).................... E(1,0,0,0) = SAB E(2,0,0,0) = 1D0 E(3,0,0,0) = 1D0 DO 50 M=1,3 C.....COMPUTE E(I,0,N).................... DO 30 I=1,NA E(m,I,0,0) = PA(m)*E(m,I-1,0,0) + E(m,I-1,0,1) DO 20 N=1,I-1 20 E(m,I,0,N)=ABI*E(m,I-1,0,N-1) + PA(m)*E(m,I-1,0,N) * + (1D0+N)*E(m,I-1,0,N+1) 30 E(m,I,0,I) = ABI*E(m,I-1,0,I-1) + PA(m)*E(m,I-1,0,I) C.....COMPUTE E(I,J,N).................... DO 50 J=1,NB DO 50 I=0,NA E(m,I,J,0) = PB(m)*E(m,I,J-1,0) + E(m,I,J-1,1) DO 40 N=1,I+J-1 40 E(m,I,J,N) = ABI *E(m,I,J-1,N-1)+ PB(m)*E(m,I,J-1,N) * + (1D0+N)*E(m,I,J-1,N+1) E(m,I,J,I+J) = ABI *E(m,I,J-1,I+J-1) + PB(m)*E(m,I,J-1,I+J) 50 CONTINUE RETURN END C################maker################################################## SUBROUTINE MAKER(R,S,T,W,PQ,NX) IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) DIMENSION R(0:4,0:4,0:4),F(0:4),PQ(3) C---------------------------------------------------------------------+ C COMPUTE AUXILLARY FUNCTION | C | C j | C R(0,0,0,j) = S(-2w) Fj(T) | C | C | C R(I,K,M,j) = R(I-1,K,M,j+1)*Cx + R(I-2,K,M,j+1)*(I-1) | C = R(I,K-1,M,j+1)*Cy + R(I,K-2,M,j+1)*(K-1) | C = R(I,K,M-1,j+1)*Cz + R(I,K,M-2,j+1)*(M-1) <2.114>| C | C---------------------------------------------------------------------+ CALL GAMMAF(F,T,NX) R(0,0,0) = F(0)*S IF (NX.EQ.0) GOTO 100 PX = PQ(1) PY = PQ(2) PZ = PQ(3) R1 = F(1)*S*W*(-2D0) R(1,0,0) = R1*PX R(0,1,0) = R1*PY R(0,0,1) = R1*PZ c print*,' maker 1 :',r(1,0,0),r(0,1,0),r(0,0,1),pq IF (NX.EQ.1) GOTO 100 R2 = F(2)*S*W**2*(4D0) R(2,0,0) = R2*PX*PX+R1 R(1,1,0) = R2*PX*PY R(0,2,0) = R2*PY*PY+R1 R(1,0,1) = R2*PX*PZ R(0,1,1) = R2*PY*PZ R(0,0,2) = R2*PZ*PZ+R1 IF (NX.EQ.2) GOTO 100 R3 = F(3)*S*W**3*(-8D0) R2X = R2*PX R2Y = R2*PY R2Z = R2*PZ R3X = R3*PX**2 R3Y = R3*PY**2 R3Z = R3*PZ**2 R(3,0,0) = R3X*PX+3D0*R2X R(0,3,0) = R3Y*PY+3D0*R2Y R(0,0,3) = R3Z*PZ+3D0*R2Z R(2,1,0) = R3X*PY+R2Y R(2,0,1) = R3X*PZ+R2Z R(1,2,0) = R3Y*PX+R2X R(0,2,1) = R3Y*PZ+R2Z R(1,0,2) = R3Z*PX+R2X R(0,1,2) = R3Z*PY+R2Y R(1,1,1) = R3 *PX*PY*PZ IF (NX.EQ.3) GOTO 100 R4 = F(4)*S*W**4*(16D0) XX = PX*PX YY = PY*PY ZZ = PZ*PZ R4X = R4*XX R4Y = R4*YY R4Z = R4*ZZ R(2,1,1) =(R4X + R3)*PY*PZ R(1,2,1) =(R4Y + R3)*PX*PZ R(1,1,2) =(R4Z + R3)*PX*PY R(2,2,0) = R4X*YY + R3X+R3Y + R2 R(2,0,2) = R4X*ZZ + R3X+R3Z + R2 R(0,2,2) = R4Y*ZZ + R3Y+R3Z + R2 R(3,1,0) =(R4X + R3*3D0)*PX*PY R(3,0,1) =(R4X + R3*3D0)*PX*PZ R(1,3,0) =(R4Y + R3*3D0)*PX*PY R(0,3,1) =(R4Y + R3*3D0)*PZ*PY R(1,0,3) =(R4Z + R3*3D0)*PX*PZ R(0,1,3) =(R4Z + R3*3D0)*PY*PZ R(4,0,0) =(R4X + R3*6D0)*XX + R2*3D0 R(0,4,0) =(R4Y + R3*6D0)*YY + R2*3D0 R(0,0,4) =(R4Z + R3*6D0)*ZZ + R2*3D0 100 RETURN END C################gammaf################################################# SUBROUTINE GAMMAF(F,T,M) IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) C--------------------------------------------------+ C COMPUTE INCOMPLETE GAMMA FUNCTION | C | C /1 2m (-Tx^2) | C Fm(T) = | x e dx | C /0 | C | C--------------------------------------------------+ PARAMETER (TMAX=20, ACCY=1.D-12, SQRPI=0.88622692545275801365D0) DIMENSION F(0:M) PT5PM = 0.5D0+M IF (T.GT.TMAX) THEN F(0) = SQRPI/SQRT(T) DO 10 L=1,M 10 F(L) = F(L-1)*(1D0*L-0.5D0)/T ELSE TI2 = T*2D0 EXPTI = EXP(-T) TERM = 0.5D0/PT5PM SUMK = TERM DO 25 K=1,1000 TERM = TERM*T/(PT5PM+K) SUMK = SUMK + TERM 25 IF (TERM.LT.ACCY) GOTO 35 35 F(M) = EXPTI*SUMK DO 50 L=M-1,0,-1 50 F(L) = (F(L+1)*TI2 + EXPTI) / (2D0*L + 1D0) ENDIF RETURN END subroutine ptintd(h,g,n) implicit real*8(a-h,o-z) include "limits.h" dimension h(mxo2) dimension g(mxo4) c c Write a canonical list of 1 and 2 electron integrals c to unit 32. c rewind(unit=32) c c Write the one-electron integrals. c il = (n*(n+1))/2 do 10 i = 1,il write(32,1000)h(i) 10 continue c c Write the two-electron integrals. c ind = 0 do 50 i = 1,n do 40 k = 1,i do 30 j = 1,i lend = j if(i.eq.j)lend = k do 20 l = 1,lend ind = ind + 1 write(32,1000)g(ind) 20 continue 30 continue 40 continue 50 continue return 1000 format(1x,f15.8) end subroutine ptintu(h,g,n) implicit real*8(a-h,o-z) character*1 ans include "limits.h" dimension h(mxo2) dimension g(mxo4) c c Write a canonical list of 1 and 2 electron integrals c to the user. c c Write the one-electron integrals. c write(*,*)' do you want to see the 1e- AO ' write(*,*)' integrals? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' the one-electron AO integrals: ' ij = 0 do 10 i = 1,n do 5 j = 1,i ij = ij + 1 write(*,1005)i,j,h(ij) 5 continue 10 continue endif c c Write the two-electron integrals. c write(*,*)' do you want to see the 2e- AO ' write(*,*)' integrals? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' the two-electron AO integrals: ' write(*,*) ind = 0 do 50 i = 1,n do 40 k = 1,i do 30 j = 1,i lend = j if(i.eq.j)lend = k do 20 l = 1,lend ind = ind + 1 write(*,1010)i,j,k,l,g(ind) 20 continue 30 continue 40 continue 50 continue endif return 1005 format(' i = ',i2,' j = ',i2,' h(i,j) = ',f15.8) 1010 format(' i = ',i2,' j = ',i2,' k = ',i2,' l = ',i2, & ' g(i,j,k,l) = ',f15.8) end