program trans c ****************** Program TRANS ***************** c This program is designed to read in the LCAO-MO c coefficient matrix, the one- and two-electron AO c integrals (in Dirac <12|12> convention), and to c transform the integrals from the AO to the MO c basis, and write these MO integrals to a file. 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 character*1 ans include "limits.h" c c Set core needed to: c 2 square arrays of order maxorb c 1 symmetric array of order maxorb c 2 arrays of dimension: c maxorb*maxorb*((maxorb*(maxorb+1))/2) c total = 2*maxorb**2 + (maxorb*(maxorb+1))/2 + c 2*maxorb*maxorb*(maxorb*(maxorb+1)/2) c dimension pcore(2*mxosq + mxo2 + 2*mxosq*mxo2) dimension ipntr(10) c c Calculate and show memory use. c itotal = 2*mxosq + mxo2 + 2*mxosq*mxo2 + 10 call modinfo('trans ') call shwmem('trans ',itotal) call goon 10 continue write(*,*) write(*,*)' Input the number of orbitals in your system: ' read(*,*)iorb if(iorb.gt.maxorb)then write(*,1010)maxorb goto 10 endif iorbsq = iorb*iorb iorb2 = iorb*((iorb+1))/2 iorb4 = iorb2*((iorb2+1))/2 ipntr(1) = 1 ipntr(2) = ipntr(1) + iorbsq ipntr(3) = ipntr(2) + iorb2 ipntr(4) = ipntr(3) + iorbsq ipntr(5) = ipntr(4) + iorbsq*iorb2 ipntr(6) = ipntr(5) + iorbsq*iorb2 c c Open the MO-AO coefficient file. c dscrpt = 'MO-AO coefficients' ids = 18 dfname = 'mocoefs.dat' idf = 11 fname = ' ' ifn = 30 iflnum = 31 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call rdmatd(fname,ifn,dscrpt,ids,pcore(ipntr(1)),iorb, & iorb,iflnum) c c Close MO coefficients file. c close(unit=31) else call rdmatu(dscrpt,ids,pcore(ipntr(1)),iorb,iorb) endif call wtmatu(dscrpt,ids,pcore(ipntr(1)),iorb,iorb) c c Open the AO integral file. c dscrpt = '1e- and 2e- AO integrals' ids = 24 dfname = 'ao_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 AO integrals available ' write(*,*)' on disk ... this can be done by running ' write(*,*)' program INTEGRAL. ' call qmexit endif 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 = 32 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) call driver(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)), & pcore(ipntr(4)),pcore(ipntr(5)),iorb,iorbsq, & iorb2,iorb4) if(yesno)then call ptintd(pcore(ipntr(2)),pcore(ipntr(5)),iorb,iorb2,iorb4) write(*,*)' The MO integrals have been written to the ', & ' file: ',FNAME c c Close 1e- and 2- MO integral file. c close(unit=32) endif write(*,*)' Do you want to see the 1e- or 2e- MO ' write(*,*)' integrals? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then call ptintu(pcore(ipntr(2)),pcore(ipntr(5)),iorb,iorb2,iorb4) endif call qmexit stop 1010 format(' Maximum number of orbitals is ',I2,'. ',/, & ' Please try again. ') end subroutine driver(c,h,p,q,g,n,nsq,n2,n4) implicit real (a-h,o-z) dimension c(n,n), h(n2), p(n,n), q(nsq*n2), g(nsq*n2) call get1e(h,n,n2) call get2e(g,n,n4) c c Close 1e- and 2- AO integral file. c close(unit=30) c c Transform the one-electron integrals. c write(*,*)' Begin 1e- integral transformation. ' call t1(c,h,p,n,nsq,n2) write(*,*)' End 1e- integral transformation. ' c c Transform the two-electron integrals. c write(*,*)' Begin 2e- integral transformation. ' c c Do the first quarter transformation. c write(*,*)' Begin 1st quarter transformation. ' call q1(c,q,g,n,nsq,n2) write(*,*)' End 1st quarter transformation. ' c c Do the second quarter transformation. c write(*,*)' Begin 2nd quarter transformation. ' call q2(c,q,g,n,nsq,n2) write(*,*)' End 2nd quarter transformation. ' c c Do the third quarter transformation. c write(*,*)' Begin 3rd quarter transformation. ' call q3(c,q,g,n,nsq,n2) write(*,*)' End 3rd quarter transformation. ' c c Do the fourth quarter transformation. c write(*,*)' Begin 4th quarter transformation. ' call q4(c,q,g,n,nsq,n2) write(*,*)' End 4th quarter transformation. ' write(*,*)' End 2e- integral transformation. ' return end subroutine q1(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) c c Do the first quarter transformation. c nl = nsq*n2 call ezero(q,nl) nsq = n*n iabsym = 0 icx = 0 do 50 ia = 1,n iam1 = ia - 1 iasq = iam1*n do 40 ib = 1,ia ibm1 = ib - 1 ibsq = ibm1*n iabsym = iabsym + 1 iabsnn = (iabsym - 1)*nsq do 30 ic = 1,ia icm1 = ic - 1 icsq = icm1*n icsym = ic*icm1/2 idend = ic if(ic.eq.ia)idend = ib do 20 id = 1,idend idm1 = id - 1 idsq = idm1*n icdsym = (icsym + id) icdsnn = (icdsym - 1)*nsq icx = icx + 1 xg = g(icx) do 10 l = 1,n issx = iabsnn + icsq + l q(issx) = q(issx) + xg*c(id,l) if(ic.ne.id)then issx = iabsnn + idsq + l q(issx) = q(issx) + xg*c(ic,l) endif if(iabsym.ne.icdsym)then issx = icdsnn + iasq + l q(issx) = q(issx) + xg*c(ib,l) if(ia.ne.ib)then issx = icdsnn + ibsq + l q(issx) = q(issx) + xg*c(ia,l) endif endif 10 continue 20 continue 30 continue 40 continue 50 continue return end subroutine q2(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) nl = nsq*n2 call ezero(g,nl) c c Do the second quarter transformation. c nsq = n*n iabsym = 0 do 150 ia = 1,n do 140 ib = 1,ia iabsym = iabsym + 1 iabsnn = (iabsym - 1)*nsq do 130 ic = 1,n icm1 = ic - 1 icsq = icm1*n do 120 l = 1,n issx = iabsnn + icsq + l xq = q(issx) do 110 k = l,n km1 = k - 1 ksq = km1*n issx = iabsnn + ksq + l g(issx) = g(issx) + xq*c(ic,k) 110 continue 120 continue 130 continue 140 continue 150 continue return end subroutine q3(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) nl = nsq*n2 call ezero(q,nl) c c Do the third quarter transformation. c nsq = n*n klsym = 0 do 250 k = 1,n km1 = k - 1 ksq = km1*n do 240 l = 1,k klsym = klsym + 1 klsnn = (klsym - 1)*nsq iabsym = 0 do 230 ia = 1,n iam1 = ia - 1 iasq = iam1*n do 220 ib = 1,ia ibm1 = ib - 1 ibsq = ibm1*n iabsym = iabsym + 1 iabsnn = (iabsym - 1)*nsq issx = iabsnn + ksq + l xg = g(issx) do 210 j = 1,n issx = klsnn + iasq + j q(issx) = q(issx) + xg*c(ib,j) if(ia.ne.ib)then issx = klsnn + ibsq + j q(issx) = q(issx) + xg*c(ia,j) endif 210 continue 220 continue 230 continue 240 continue 250 continue return end subroutine q4(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) nl = nsq*n2 call ezero(g,nl) c c Do the fourth quarter transformation. c nsq = n*n klsym = 0 do 350 k = 1,n do 340 l = 1,k klsym = klsym + 1 klsnn = (klsym - 1)*nsq do 330 ia = 1,n iam1 = ia - 1 iasq = iam1*n do 320 j = 1,n issx = klsnn + iasq + j xq = q(issx) istart = max(j,k) if(l.gt.j)istart = k + 1 do 310 i = istart,n icx = icanon(i,k,j,l) g(icx) = g(icx) + xq*c(ia,i) 310 continue 320 continue 330 continue 340 continue 350 continue return end subroutine t1(c,h,p,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n),h(n2),p(n,n) c c Transform the one-electron integrals. c nl = nsq call ezero(p,nl) iaib = 0 do 430 ia = 1,n do 420 ib = 1,ia iaib = iaib + 1 do 410 j = 1,n p(ia,j) = p(ia,j) + h(iaib)*c(ib,j) if(ia.ne.ib)then p(ib,j) = p(ib,j) + h(iaib)*c(ia,j) endif 410 continue 420 continue 430 continue nl = n2 call ezero(h,nl) do 530 ia = 1,n ij = 0 do 520 j = 1,n do 510 i = 1,j ij = ij + 1 h(ij) = h(ij) + p(ia,j)*c(ia,i) 510 continue 520 continue 530 continue return end subroutine ptintd(h,g,n,n2,n4) implicit real (a-h,o-z) dimension h(n2), g(n4) 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 = n2 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,n2,n4) implicit real (a-h,o-z) character*1 ans dimension h(n2), g(n4) 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- MO ' write(*,*)' integrals? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' The one-electron MO integrals: ' call tprntd(h,n) endif c C Write the two-electron integrals. C write(*,*)' Do you want to see the 2e- MO ' write(*,*)' integrals? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' The two-electron MO 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 1000 format(1x,f15.8) 1010 format(' i = ',i2,' j = ',i2,' k = ',i2,' l = ',i2, & ' g(i,j,k,l) = ',f15.8) end