program mp2 c ****************** Program MP2 ******************* c This program is designed to read in the c transformed one- and two-electron integrals and c the Fock orbital energies after which it will c compute the second order Moller Plesset c perturbation theory energy (MP2). c c Copyright (C) 1992 Jeff Nichols and Jack Simons c c ************************************************** implicit real(a-h,o-z) logical yesno character*30 dfname, fname, dscrpt include "limits.h" c c Set core needed to: c 1 vector of dimension MAXORB c 2 symmetric arrays of order MAXORB c 1 array of dimension canonical(MAXORB**4) c total = 2*(maxorb*(maxorb+1))/2 + c (maxorb*(maxorb+1))/2*(((maxorb*(maxorb+1))/2+1))/2 + c maxorb c dimension pcore(2*mxo2 + mxo4 + maxorb) dimension ipntr(10) itotal = 2*mxo2 + mxo4 + maxorb + 10 call modinfo('mp2 ') call shwmem('mp2 ',itotal) call goon 1 continue write(*,*) write(*,*)' Input the number of orbitals in your system: ' read(*,*)iorb if(iorb.gt.maxorb)then write(*,1010)maxorb goto 1 endif 5 continue write(*,*)' Please input the number of doubly occupied ' write(*,*)' orbitals. ' read(*,*)ndocc write(*,1000)ndocc if(ndocc.gt.iorb)then write(*,*)' The number of doubly occupied orbitals cannot ' write(*,*)' exceed the total number of orbitals. Please ' write(*,*)' try again. ' goto 5 endif iorb2 = iorb*(iorb+1)/2 iorb4 = iorb2*(iorb2+1)/2 ipntr(1) = 1 ipntr(2) = ipntr(1) + iorb2 ipntr(3) = ipntr(2) + iorb4 ipntr(4) = ipntr(3) + iorb ipntr(5) = ipntr(4) + iorb2 c c Open the MO integral file. c dscrpt = '1e- and 2e- MO integrals' ids = 24 dfname = 'mo_integrals.dat' idf = 16 fname = ' ' ifn = 30 iflnum = 30 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(.not.yesno)then write(*,*)' You must have the MO integrals available ' write(*,*)' on disk to construct the MP2 energy. ' call qmexit endif c c Get the 1e- AO integrals. c call get1e(pcore(ipntr(1)),iorb,iorb2) c c Get the 2e- AO integrals. c call get2e(pcore(ipntr(2)),iorb,iorb4) c c Close 1e- and 2- MO integral file. c close(unit=30) c c Open the Fock orbital energy file. c dscrpt = 'Fock orbital energies' ids = 21 dfname = 'orbital.evals' idf = 13 fname = ' ' ifn = 30 iflnum = 33 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(.not.yesno)then write(*,*)' You must have the Fock orbital energies available ' write(*,*)' on disk to construct the MP2 energy. ' call qmexit endif c c Get the orbital energies. c rewind(unit=33) do 10 i = 1,iorb read(33,*)pcore(ipntr(3)+i-1) 10 continue c c Close orbital energies file. c close(unit=33) c c Compute MP2 energy. c call domp2(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)), & pcore(ipntr(4)),ndocc,iorb,iorb2,iorb4) call qmexit stop 1000 format(' There is(are) ',i2,' doubly occupied orbital(s). ') 1010 format(' Maximum number of orbitals is ',I2,'. ',/, & ' Please try again. ') end subroutine domp2(h,g,e,ep,nc,n,n2,n4) implicit real (a-h,o-z) include "limits.h" dimension h(n2), g(n4), e(n), ep(n2) 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)) call ezero(ep,n2) c c need to do all spin combos e.g., aaaa, bbbb, abab, etc. c ib, ia, is, and ir label spin orbitals c ibspac, iaspac, isspac, and irspac label spacial orbitals c ibs, ias, iss, and irs give alpha or beta labels, 0=alpha, 1=beta c do 40 ib = 1, 2*nc ibspac = (ib+1)/2 if(ib-2*(ib/2).eq.0)then ibs=1 else ibs=0 endif do 30 ia = 1, ib-1 iaspac = (ia+1)/2 if(ia-2*(ia/2).eq.0)then ias=1 else ias=0 endif ibsias = isym2(ibspac,iaspac) do 20 is = 2*nc+1, 2*n isspac = (is+1)/2 if(is-2*(is/2).eq.0)then iss=1 else iss=0 endif do 10 ir = 2*nc+1, is-1 irspac = (ir+1)/2 if(ir-2*(ir/2).eq.0)then irs=1 else irs=0 endif iabrs = isym4(iaspac,irspac,ibspac,isspac) iabsr = isym4(iaspac,isspac,ibspac,irspac) if(ias.eq.irs)then arspin = 1.0e0 else arspin = 0.0e0 endif if(ibs.eq.iss)then bsspin = 1.0e0 else bsspin = 0.0e0 endif if(ias.eq.iss)then asspin = 1.0e0 else asspin = 0.0e0 endif if(ibs.eq.irs)then brspin = 1.0e0 else brspin = 0.0e0 endif giabrs = g(iabrs)*arspin*bsspin giabsr = g(iabsr)*asspin*brspin et = ((giabrs-giabsr)**2)/ & (e(irspac)+e(isspac)-e(iaspac)-e(ibspac)) c write(*,*)' ia,ib,ir,is,iabrs,iabsr,giabrs,giabsr,et ', c & ia,ib,ir,is,iabrs,iabsr,giabrs,giabsr,et ep(ibsias) = ep(ibsias) - et 10 continue 20 continue 30 continue 40 continue write(*,*)' Pair correlation energies: ' call tprntd(ep,nc) emp2 = 0.0e0 do 50 i = 1, n2 emp2 = emp2 + ep(i) 50 continue write(*,*)' MP2 energy = ',emp2 return end