|
Columbus
|
|
00README,
Columbus.22-jun-92.tar.Z,
Columbus.new.06-oct-94.tar.Z,
README.ccl,
aogrd.tar.Z,
autbas.tar.Z,
batra.tar.Z,
ciden.doc,
ciden.tar.Z,
cidrt.tar.Z,
cigrd.tar.Z,
ciudg.tar.Z,
ciuft.tar.Z,
columbus.new.17-mar-95.tar.Z,
gdiis.doc,
gdiis.f,
global1.3.1.tar.Z,
global2.0.tar.Z,
hermit.tar.Z,
makefile,
mcdrt.tar.Z,
mcscf.tar.Z,
mcuft.tar.Z,
newbas.tar.Z,
property.doc,
property.f,
tcgmsg.4.03.tar.Z
|
|
|
c
c pa60 - polyatom properties program for s,p,d and f-orbitals
c
c 250 basis functions allowed
c
c input modified according to the karlsruhe version of pitzer's
c scf-program by c.m.ehrhardt
c for this purpose the subroutines inpsym,getran,scfinp and comtrn
c are taken from h.f. scfaefer's gradient-program
c
c the subroutines rhftce,propa and prod were made faster (4.3.1984)
c
c the program was modified to handle f-functions by p.maier and
c c.m.ehrhardt ( 1.3.1984 )
c
c mulliken and roby population analysis included ( 12.7.1984 )
c
c modified to calculate the kinetic energy and its relativistic
c correction by s. brode in april,1984 : a subroutine opap4 was
c added and the main program and the subroutines dintfm, pint,
c propa and prnt were slightly modified
c
c boys localization included by r.ahlrichs and h.schiffer
c (1.3.1985)
c
c ft02,ft04,ft07,ft09 scratch files
c ft22 localized mos for plopr
c ft21 electrostatic potential
c ft23 file for plot
c
c input description:
c
c input files:
c ft24 integral information from the karlsruhe integral input
c program
c ft08 eigenvectors and occupation as written out by subroutine
c outgrd of scf-program
c ft05 input unit
c ft28 optional file for atomic coordinates (compatible to ints)
c
c input on unit ft05:
c
c 1. title (format 12a6)
c
c 2. isl1 , isl2 , isl4 , nsupp (format 4i5)
c isl1 = 1 full property calculation (default)
c isl1 = 2 property only in ao-basis
c isl2 = 1 no output of property matrix in ao-basis(default)
c isl2 = 2 output of these integrals
c isl2 = 3 output, and integrals are written on unit 3
c isl4 = 1 renormalisation of wavefunctions and output
c isl4 = 2 no renormalisation
c isl4 = all other values renormalisation but no output
c (default)
c nsupp= 1 output of diagonalelements of property matrix
c in mo-basis
c nsupp= 2 controll-output of integral and scf input data
c nsupp= 3 output of mo vectors in ao basis
c nsupp= 4 output of overlap population per mo
c nsupp= 0 no output (default=0)
c
c 3. m, ilnmc, icntr, (c(k),k=1,3) (format i5,6x,a4,i5,3f10.5)
c step 3 can be repeated several times
c population analysis or localization must always be the
c last option
c m specifies the property :
c potential(1), electric field(2), field gradient(3),
c dipole moment(4), quadrupole moment(5), third moment(6),
c planar density(7), linear density(8), charge density(9),
c overlap(10), second moment(11), diamagnetic shielding(12),
c relativistic correction(13), mulliken and roby population
c analysis(14), boys localization(15), plot(16)
c
c ilnmc name of the point at which property shall be
c calculated
c icntr number of atom at which property shall be
c calculated, if icntr.ne.0 you need not read in
c coordinates
c c(k) cartesian coordiantes of the point
c
c 4. special input for population analysis (m=14)
c
c nmao, (mor(im),im=1,nmao) (format i10, 20 i3)
c nmao number of maos for atom 1,2...,n
c mor numbers of selected maos
c mor(1)= 0 automatic selection
c nm, (iatnr(i),i=1,nm) (format 20i3)
c nm.ne.0 special multicenter term of order nm
c iatnr(i) atoms that form the multicenter term
c (a blank line is required if nm=0)
c
c 5. special input for localization (m=15)
c nmo (format 5x,i5)
c nmo = -1 all occupied mo's shall be localized
c nmo = 0 now the program wants to read a threshold
c for the selection of mo's, wich shall be
c localized
c nmo = nmo number of mo's, wich shall be localized
c now the program wants to read the corresponding
c mo numbers in ascending order
c either thr (5x,f10.4) for nmo = 0
c or mon(i),i=1,nmo (5x,(5i5)) for nmo > 0
c nop (5x,i5)
c nop = 0 no further input needed for localization
c nop = nop number of operators ( 1,2 or 3 )
c now the program wants to read the corresponding
c operator numbers ( x=1, y=2, z=3 ) in canonical
c order
c opn(i),i=1,nop (5x,3i5) for nop > 0
c
c
c important dimensions
c
c number of atoms: nocmx = 50
c number of primitive basis functions : nbmx = 720
c number of contracted basis functions
c and number of mos : nomx = 250
c dimension of buff : lbuff = nomx*(nomx+1)/2
c dimensions in common /big/ : maxuv = nomx**2
c max number of orbitals to be local. : mloc = 250
c (in common eig and scfinf)
c dimension of rnabcd : n*(n-3)/2+m
c with: n=nocmx*(nocmx-3)/2+nocmx-1
c m=(nocmx-2)*(nocmx-5)/2+nocmx-3
c max number of grids for 'x' direction: nmaxx = 50
c max number of grids for 'y' direction: nmaxy = 50
c max number of orbitals to be ploted : norbmx = 10
c
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375,nrabcd=748379)
parameter(nmaxx=50,nmaxy=50,norbmx=10,icomx=12)
real*8 icod(12),ilab(12),ilbl(12)
dimension eta(nbmx,25),fcntr(nbmx),
1 ftype(nbmx),frocc(nomx),c(3),dint(3,10),
2 type(20),val(10),tot(10),nstart(nomx),nstop(nomx),
3 ncontc(nomx),buff(lbuff),a(39),title(150),
4 cint(nomx,10),izahl(nbmx),itrnao(nbmx),idxxt(nbmx),
5 rnabcd(nrabcd),contrc(nomx,icomx)
dimension nocc(norbmx),cc(nmaxx*nmaxy,3),
1 dd(nmaxx*nmaxy,3),dinx(nmaxx*nmaxy),
2 cin(nmaxx*nmaxy,norbmx+1),
3 v(nmaxx*nmaxy),ddsq(nmaxx*nmaxy)
equivalence (buff(nocmx+1),rnabcd(1))
external opaa0,opaa1,opaa2,opab1,opab2,opab3,opac1,opac2, opac3,
1opad1,opab4,opaa3,opap4
external potnuc,fldnuc,grdnuc,dipnuc,qudnuc,octnuc,nonuc,magnuc,
1sldnuc
integer centr,fcntr,ftype,type
common/ioind/isl1, isl2
common/muli/vlist(nocmx,4),ianr(nomx),centr(nocmx),itnr(nomx)
common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
common/supp/nsupp
common/ncnt/ilnmcm
common/prptyp/m
common /contrc/ chomo,ncont
common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
c
data nromx/3/,nfcol/25/,ntcol/20/,iwarn/0/,ir/8/
data icod/6h1pa60 ,6h ,6h ,6hnyu pr,6hoperti,6hes pac,
1 6hkage ,6h ,6h ,6h ,6h12/66 ,6hspdf /
data title/6h 1/r,6h ex,6h fxx,6h x,6h qxx,6h oxxx,
1 6h dx,6h dxy,6h den,6h ovl,6h xx,6h fx*x*,6h qx*x*,
2 6h x*x*,6h, 6h ,6h ey,6h fyy,
3 6h y,6h qyy,6h oyyy,6h dy,6h dxz,6h ,6h ,
4 6h yy,6h fy*y*,6h qy*y*,6h y*y*,6h,
5 6h ,6h ez,6h fzz,6h z,6h qzz,6h ozzz,6h dz,
6 6h dyz,6h ,6h ,6h zz,6h fz*z*,6h qz*z*,6h z*z*,
7 6h sum , 6h ,6h ,6h fxy,6h ,
8 6h qxy,6h oxxy,6h ,6h ,6h ,6h ,6h xy,
9 6h fx*y*,6h qx*y*,6h x*y*,6h , 6h ,
1 6h ,6h fxz,6h ,6h qxz,6h oxxz,6h ,6h ,
2 6h ,6h ,6h xz,6h fx*z*,6h qx*z*,6h x*z*,6h ,
3 6h ,6h ,6h fyz,6h ,6h qyz,
4 6h oxyy,6h ,6h ,6h ,6h ,6h yz,6h fy*z*,
5 6h qy*z*,6h y*z*,6h ,
6 4*6h ,6h rsq,6h oyyz,9*6h ,
7 5*6h ,6h oxzz,9*6h ,
8 5*6h ,6h oyzz,9*6h ,
9 5*6h ,6h oxyz,9*6h /
data a/6hpotent,6helectr,6hfield ,6hdipole,6hquadru,6hthird ,
1 6hplanar,6hline d,6hcharge,6hoverla,6hsecond,6hdiamag,
& 6hrelati,
2 6hial ,6hic fie,6hgradie,6h momen,6hpole m,6hmoment,
3 6h densi,6hensity,6h densi,6hp ,6h momen,6hnetic ,
& 6hvistic,
4 6h ,6hld ,6hnt ,6ht ,6homent ,6h ,
5 6hty ,6h ,6hty ,6h ,6ht ,6hshldg ,
& 6h corr./
data type/4hs ,4hx ,4hy ,4hz ,4hxx ,4hyy ,4hzz ,
1 4hxy ,4hxz ,4hyz ,4hxxx ,4hyyy ,4hzzz ,4hxxy ,4hxxz ,
2 4hyyx ,4hyyz ,4hzzx ,4hzzy ,4hxyz /
c
c** call xuflow(0)
c file shout short output
open(unit=5,file='propin')
open(6,file='propls')
nshout=13
open(unit=nshout,file='shout',status='unknown'
&,access='sequential',form='formatted')
rewind nshout
c file vectgrd scf vectors
c ir=8
open(unit=ir,file='vectgrd',status='old',
1access='sequential',form='formatted')
c file plot: plot output
iplot=23
open(unit=iplot,file='plot',status='unknown',form='formatted')
rewind iplot
c
maxuv=nomx*nomx
lomao=22
76 do 20 i=1,12
20 ilab(i) = icod(i)
do 22 i=1,nomx
22 frocc(i) = 0.0d0
read(5,101,end=999) (ilbl(i), i=1,12)
write(6,101) (ilab(i), i=1,12)
write(6,103)
write(6,104) (ilbl(i), i=1,12)
read(5,100) isl1,isl2,isl4,nsupp
if(isl1.eq.0) isl1=1
if(isl2.eq.0) isl2=1
if(isl4.eq.0) isl4=3
if(isl2.eq.3) rewind 3
call input(eta,nfcol,fcntr,ftype,n,nbmx,vlist,centr,non,nocmx,
1 ntcol,type,ilab,ilbl,nomx,nstart,
2 nstop,contrc,ncontc,buff,izahl,nbfao,icomx,ianr,itnr)
call nrml(eta,nbmx,nfcol,nomx,contrc,nstart,nstop,fcntr,ftype,
1isl4,c,opad1)
c
c ivmax.ne.0 : transformation of mo-vectors in ao basis has already
c been done
rewind ir
read (ir,940)(alabel(i1),i1=1,19)
read(ir,940)(blabel(i1),i1=1,19)
read(ir,950)iopen,nsym,ivmax
read(ir,940)(fmt(i1),i1=1,4)
if(ivmax.eq.0) goto 170
if(ivmax.gt.maxuv) goto 1500
read(ir,fmt) (y(i),i=1,ivmax)
norb=iopen
do 180 nn=1,norb
frocc(nn)=2.0d0
180 chomo=chomo-2.0d0
write(6,970) chomo
goto 25
c
170 call getran(nbfao,frocc,norb,itrnao,idxxt)
25 read(5,102,end=999)m,ilnmcm,icntr,(c(k),k=1,3)
if(iwarn.eq.1) goto 998
if(m.eq.14.or.m.eq.15) iwarn=1
if(icntr.eq.0) go to 800
do 801 i=1,3
801 c(i)=vlist(icntr,i)
800 if(m.eq.0) go to 76
go to (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),m
1 call dintfm(opaa0, potnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
2 call dintfm(opaa1, fldnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
3 call dintfm(opaa2, grdnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
4 call dintfm(opab1, dipnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
5 call dintfm(opab2, qudnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
6 call dintfm(opab3, octnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
7 call dintfm(opac1, nonuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
8 call dintfm(opac2, nonuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
9 call dintfm(opac3, nonuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
10 call dintfm(opad1, nonuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
11 call dintfm(opab4, magnuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
12 call dintfm(opaa3,sldnuc, dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
13 call dintfm(opap4, nonuc, dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
14 call dintfm(opad1, nonuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
15 call dintfm(opab1, nonuc,dint, vlist,c,n,eta,val,nt
1,norb,frocc,tot,non,nromx,ilbl,a,title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
go to 25
16 call plotin(nmaxx,nmaxy,kpmax,cc,norbpl,mdens,nocc,norbmx,iplot,
1 iat)
if(iat.gt.0) call molcor(iplot)
call orbval(eta,nbmx,norbmx,nomx,nstart,nstop,contrc,ftype,
1frocc,cc,dd,dinx,cin,v,ddsq,kpmax,norbpl,mdens,nocc,y,iplot)
go to 25
999 stop
1500 write(6,960) ivmax,maxuv
stop
100 format(6i5)
101 format(12a6)
102 format(i5,6x,a4,i5,3f10.0)
103 format(/)
104 format(1x,12a6)
105 format(2x,10a4)
940 format(19a4)
950 format(20i4)
960 format(//' dimension of y in main program too small ! ivmax= ',
&i6,' maxuv= ',i6//)
970 format(/' charge of the molecule ',f10.5/)
c
998 write(6,106)
106 format(1x,'option 14 (population analysis) or option 15 (localizat
1ion) must be the last one'/,' ***** calculation stopped *****')
stop
end
c
c
subroutine input( eta, nfcol, fcntr, ftype, n, nbmx, vlist, centr,
1 non,nocmx,ntcol,type,ilab,ilbl,nomx,
2 nstart,nstop,contrc,ncontc,buff,izahl,nbfao,icomx,ianr,
3 itnr)
implicit real*8(a-h,o-z)
c
c parameter statement as in integral program
parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
& lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
& lpru=18000,lsu=50,lstu=8,lconu=12,
& lcu=8,lccu=182)
c
c this subroutine prepares the integral input from the karlsruhe
c input program for polyatom.
c one part is taken from routine convrt of schaefer's gradient
c program
c
real*8 ilab(12), ilbl(12)
integer fcntr, ftype, centr, type
dimension eta(nbmx,nfcol), vlist(nocmx,4),ianr(1),itnr(1)
dimension fcntr(1), ftype(1), centr(1), type(1),contrc(nomx,8)
dimension nstart(1),nstop(1),ncontc(1),buff(1),frocc(1),izahl(1)
common /contrc/ chomo,ncont
common /ioind/isl1,isl2
common/infoa/natoms,mul,nbf,nij,nelec,nalpha,nbeta
common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
1 ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
2 ccoef(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
3 z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
4 la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
5 mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu)
c
msu=lsu
mconsu=lconsu
mconu=lconu
mcu=lcu
msfu=lsfu
mgu=lgu
mnsfu=lnsfu
mpru=lpru
mccu=lccu
kaords=laords
mgcsu=lgcsu
mstu=lstu
mru=lru
mcsu=lcsu
mctu=lctu
msfru=lsfru
call inpsym(mcu,msu,mru,mctu,mcsu,icomx,mgcsu)
n=0
nbf=0
ispf=0
natoms=0
nr=0
do 500 is=1,ns
icu=nc(is)
ifu=nf(is)
do 400 ic=1,icu
isf=ispf
natoms=natoms+1
if(natoms.le.nocmx) goto 504
write(6,1200)
1200 format(1x,'too many atoms in subroutine input')
stop
504 centr(natoms)=mtype(is)
vlist(natoms,1)=x(ic,is)
vlist(natoms,2)=yy(ic,is)
vlist(natoms,3)=z(ic,is)
vlist(natoms,4)=chg(is)
chomo=chomo+chg(is)
nr=nr+1
do 400 if=1,ifu
isf=isf+1
icons=mcons(isf)
ncns=ncon(icons)
jtyp=lmnp1(icons)
if(jtyp.le.4) goto 11
write(6,1000) jtyp
1000 format(1x,' jtyp = ',i3,4x,' only s,p,d and f orbitals are allowed
1 ',/,'****** calculation stopped ******')
stop
11 continue
jtypm=jtyp*(jtyp+1)/2
jtyp = (jtypm*(jtyp-1))/3 +1
do 400 i=1,jtypm
nbf=nbf+1
if(nbf.lt.nomx) goto 503
write(6,1100) nbf,nomx
1100 format(1x,' dimension too small in input, nbf =',i5,3x,
1' nomx =',i5)
stop
503 ncontc(nbf)=ncns
itnr(nbf)=jtyp-1+i
ianr(nbf)=nr
do 400 icon=1,ncns
n=n+1
if(n.le.nbmx) goto 505
write(6,1300)
1300 format(1x,'too many primitive functions in subroutine input')
stop
505 izahl(n)=nr
fcntr(n)=mtype(is)
ftype(n)=type(jtyp-1+i)
eta(n,4)=zet(icon,icons)
contrc(nbf,icon)=ccoef(icon,icons)
400 continue
ispf=ispf+ifu
500 continue
non=natoms
ncont=nbf
nbfao=nbf
call rdetac(eta,nfcol,fcntr,ftype,n,nbmx,vlist,centr,non,
1 nocmx,ntcol,type,nstart,nstop,contrc,ncontc,nomx,
2 izahl)
return
end
c
c
subroutine rdetac(eta,nfcol,cntr,typ,n,nfmx,vlist,centre,non,
1ncmx,ntcol,type,nstart,nstop,contrc,ncontc,nomx,izahl)
c
c polyatom inputroutine for basis functions
c
implicit real*8(a-h,o-z)
common /contrc/ chomo,ncont
dimension ncontc(1),nstart(1),nstop(1),contrc(nomx,8)
dimension eta(nfmx,nfcol),vlist(ncmx,4)
dimension cntr(1),typ(1),centre(1),type(1),izahl(1)
integer cntr,typ,centre,blank,type
write(6,105)
c.....center specifications
write(6,104)
write(6,106)
do 1 i=1,non
write(6,107) centre(i), (vlist(i,j), j=1,4)
1 continue
if(ncont.ne.n) go to 202
do 201 i=1,n
201 ncontc(i) = 1
202 nstart(1)=1
nstop(1)=ncontc(1)
do 31 i=2,ncont
nstart(i)=nstart(i-1)+ncontc(i-1)
31 nstop(i)=nstop(i-1)+ncontc(i)
c.....set up coordinates of function in eta
icon=1
jcon=0
do 3 i=1,n
if(i-nstop(icon)) 60,60,61
60 jcon=jcon+1
go to 62
61 icon=icon+1
jcon=1
62 j=izahl(i)
do 6 m=1,3
6 eta(i,m)=vlist(j,m)
c.....set up unnormalized hybridization coeffs of function in eta
do 7 j=1,ntcol
if(noteqv(type(j),typ(i))) 7, 8, 7
7 continue
write(6,2000) type(j),typ(i)
2000 format('0***** types not equal ',2a4)
stop
8 continue
do 9 m=1,ntcol
9 eta(i,m+5)=0.0d0
do 2 m=1,ntcol
if(typ(i).ne.type(m)) goto 2
eta(i,5)=m
eta(i,m+5)=1.0
goto 3
2 continue
3 continue
return
101 format(2i3)
1020 format(a4,6x,4f12.0)
1030 format(36i2)
103 format(a4,6x,a4,i3,3x,2f12.0)
104 format(/ ,1x,15hnuclear centres)
105 format(///,1x,32hgaussian function specifications)
106 format(/ ,1x,6hcentre,18x,12hco-ordinates,17x,6hcharge,//)
107 format(1x,a6,6x,3f12.8,6x,f4.1)
110 format(///,1x,8hfunction,3x,11hcontraction,3x,6hcentre,3x,4htype,
111x,8hexponent,3x,11hcoefficient)
111 format(1x,i5,8x,i5,8x,i3,a3,3x,a6,3x,2f13.7)
114 format(//,1x,18hadditional centres,//)
115 format(//,1x,6hcentre,18x,12hco-ordinates//)
end
c
c
subroutine nrml(eta,nfmx,nfcol,nomx,c,nstart,nstop,nc,nt,isl4,d,
1ovl)
implicit real*8 (a-h,o-z)
dimension a(20),eta(nfmx,nfcol),c(nomx,8),nstart(1),nstop(1),
1nc(1),nt(1),d(3)
data cons/15.7496099d0/,c1/1.0d0/,zero/0.d0/,c3/3.0d0/,c4/4.0d0/,
1 delt/1.0d-6/,c15/15.0d0/
common /contrc/chomo,ncont
c normalize primitive functions.
do 5 i=1,ncont
is=nstart(i)
if=nstop(i)
do 4 ii=is,if
z=eta(ii,4)*c4
do 1 m=1,20
1 a(m)=eta(ii,m+5)
s=cons*(a(1)+(a(2)+a(3)+a(4)+(a(8)+a(9)+a(10)+c3*(a(5)+a(6)+a(7))+
1 (a(20)+c3*(a(14)+a(15)+a(16)+a(17)+a(18)+a(19))+c15*(a(11)+a(12)
2 +a(13))) /z)/z)/z)/(z*sqrt(z))
s=c1/sqrt(s)
do 3 m=6,nfcol
3 eta(ii,m)=eta(ii,m)*s
4 continue
5 continue
if(isl4.eq.2) return
c renormalize contraction coefficients
d(1)=zero
d(2)=zero
d(3)=zero
if(isl4.ne.1) go to 2
write(6,101)
write(6,110)
2 do 9 i=1,ncont
is=nstart(i)
if=nstop(i)
ip=0
s=zero
do 7 ii=is,if
ip=ip+1
jp=0
do 7 jj=is,if
call propa(ovl,ii,jj,eta,a,d,1,nfmx)
jp=jp+1
s=s+c(i,ip)*c(i,jp)*a(1)
7 continue
s=c1/sqrt(s)
ip=0
do 8 ii=is,if
ip=ip+1
c(i,ip)=c(i,ip)*s
if(abs(c(i,ip)+c1).le.delt) c(i,ip)=c1
if(isl4.ne.1) go to 8
write(6,111)ii,i,nc(ii),nt(ii),eta(ii,4),c(i,ip)
8 continue
9 continue
110 format(//,1x,8hfunction,3x,11hcontraction,3x,6hcentre,3x,4htype,
111x,8hexponent,3x,11hcoefficient)
111 format(1x,i5,8x,i5,8x,a6,3x,a6,3x,2f13.7)
101 format(///,10x,15hrenormalization)
return
end
c
c
subroutine dintfm(opa, nuc, dint, vlist, c, n, eta,
1 val, nt, norb, frocc, tot, noc, nro, ilbl, a, title,
2 nstart,nstop,contrc,cint,buff,rnabcd)
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
common /contrc/ chomo,ncont
common /prptyp/ m
common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
dimension dint(nro,*), vlist(nocmx,4), c(3), eta(nbmx,16), val(1),
1 frocc(1), tot(1)
dimension nstart(1),nstop(1),contrc(nomx,8),cint(nomx,1),
1 buff(lbuff),rnabcd(1)
real*8 ilbl
common /ioind/isl1,isl2
dimension ilbl(1), a(1), title(1), dinx(310)
dimension lin(84),min(84),nin(84)
data pi/3.14159265359d0/, cl/137.03604d0/
data lin/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,1,0,1,0,1,4,0,0,3,3,1,0,1,0
1 ,2,2,0,2,1,1,5,0,0,3,3,2,2,0,0,4,4,1,0,0,1,1,3,1,2,2,1,6,0,0,3,
2 3,0,5,5,1,0,0,1,4,4,2,0,2,0,3,3,1,2,2,1,4,1,1,2/,min/0,0,1,0,0,
3 2,0,1,0,1,0,3,0,1,0,2,2,0,1,1,0,4,0,1,0,3,3,0,1,2,0,2,1,2,1,0,5,
4 0,2,0,3,0,3,2,1,0,4,4,1,0,1,1,3,2,1,2,0,6,0,3,0,3,1,0,0,1,5,5,2,0
5,0,2,4,4,2,1,3,1,3,2,1,4,1,2/,nin/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,0,
6 1,2,2,1,0,0,4,0,1,0,1,3,3,0,2,2,1,1,2,0,0,5,0,2,0,3,2,3,0,1,0,1,
7 4,4,3,1,1,1,2,2,0,0,6,0,3,3,0,1,5,5,1,0,0,2,4,4,0,2,1,2,2,3,1,3,
8 1,1,4,2/
c
c
call dumdum(nuc,noc,vlist,c,dint,nro,nocmx,nt)
if(m.eq.15) nt=3
if(m.eq.13) nt=3
call pint(opa, n, eta, c, val, nt, nro, nbmx, a, title,
1 ilbl,nstart,nstop,contrc,nomx,buff,lbuff)
if (m.ne.13) go to 7
c
c calculation of the darwin-term
c
valval=0.0d0
c + loop over the centers
do 80 ic=1,noc
c(1)=vlist(ic,1)
c(2)=vlist(ic,2)
c(3)=vlist(ic,3)
ipp=0
c + + loop over the c-gto's
do 40 iz=1,ncont
dinx(iz)=0.d0
istart=nstart(iz)
istop=nstop(iz)
iprim=0
c + + + loop over the gto's
do 40 ii=istart,istop
ipp=ipp+1
iprim=iprim+1
ieta=eta(ii,5)+0.00001
li=lin(ieta)
mi=min(ieta)
ni=nin(ieta)
xc=c(1)-eta(ii,1)
yc=c(2)-eta(ii,2)
zc=c(3)-eta(ii,3)
expo=-eta(ii,4)*(xc*xc+yc*yc+zc*zc)
din=1.d0
if(li.ne.0) din=xc**li
if(mi.ne.0) din=din*yc**mi
if(ni.ne.0) din=din*zc**ni
din=din*exp(expo)
ieta=ieta+5
dinx(iz)=dinx(iz)+din*contrc(iz,iprim)*eta(ipp,ieta)
40 continue
vval=0.d0
ii=0
c + + loop over the orbitals
do 70 mi=1,norb
dens=0.d0
c + + + loop over the c-gtos
do 75 mj=1,ncont
ii=ii+1
dens=dens+dinx(mj)*y(ii)
75 continue
vval=vval+frocc(mi)*dens*dens
70 continue
valval=valval+vlist(ic,4)*vval
80 continue
dint(1,3)=valval*0.5d0*pi/cl**2
c
c
7 if (m.ne.14) go to 6
c
c mulliken and roby population analysis
c
call cntrcm(cint,norb,noc,buff,frocc,rnabcd)
return
c
6 go to (4,5), isl1
4 call output(dint,norb,n,nt,nro,frocc,tot,c,ilbl,a,title,
1 cint,buff)
5 return
end
c
c
subroutine pint(opa,n,eta,c,val,nt,nro,nbmx,a,title,
1 ilbl,nstart,nstop,contrc,nomx,buff,lbuff)
c
c calculation of integrals in ao-basis. data are written on
c file ft002
c
implicit real*8(a-h,o-z)
common /contrc/ chomo,ncont
dimension nstart(1),nstop(1),contrc(nomx,8)
dimension a(13,3), title(14,10)
dimension c(1),eta(nbmx,16),val(1),dinx(10),buff(lbuff)
real*8 ilbl(1)
common/ioind/isl1,isl2
common/prptyp/mm
ipr=ncont*(ncont+1)/2
if(ipr.gt.lbuff.and.mm.eq.14) goto 99
nu = 3
rewind 2
ix=0
do 50 i=1,ncont
do 50 j=1,i
nu = nu + 1
istart = nstart(i)
istop = nstop(i)
jstart = nstart(j)
jstop = nstop(j)
iprim = 0
do 2 m=1,nt
2 dinx(m)=0.0d0
do 40 ii=istart,istop
iprim = iprim + 1
jprim = 0
do 40 jj=jstart,jstop
jprim = jprim + 1
call propa(opa,ii,jj,eta,val,c,nt,nbmx)
do 30 m=1,nt
30 dinx(m) = contrc(i,iprim)*contrc(j,jprim)*val(m) + dinx(m)
40 continue
if(mm.ne.14) goto 55
ix=ix+1
buff(ix)=dinx(1)
goto 50
55 do 60 m=1,nt
ix=ix+1
buff(ix)=dinx(m)
if(ix.lt.lbuff)go to 60
ix=0
write(2)buff
60 continue
50 continue
if(mm.eq.14) return
write(2)buff
if(isl2.eq.1)return
if(isl2.eq.3) go to 1
write(6,19) (ilbl(i), i=1,12)
write(6,20) (a(mm,i), i=1,3)
write(6,21)
rewind 2
ix=lbuff
do 6 i=1,ncont
do 7 j=1,i
write(6,21)
do 70 m=1,nt
ix=ix+1
if(ix.le.lbuff)go to 70
ix=1
read(2)buff
70 write(6,10)i,j,buff(ix)
7 continue
6 continue
return
1 rewind 2
ix=lbuff
do 66 i=1,ncont
do 77 j=1,i
do 770 m=1,nt
ix=ix+1
if(ix.le.lbuff) go to 770
ix=1
read(2) buff
write(3) buff
770 continue
77 continue
66 continue
return
c
99 write(6,7000) ipr,lbuff
7000 format(1x,'dimension too small in subroutine pint ipr =',i8,
1 ' lbuff = ',i8,/,' population analysis is not possible')
stop
c
19 format(1h0,/ ,1x,12a6)
20 format( //,1x,34hintegrals over basis functions for,5x,3a6//)
21 format(/)
10 format(1x,2i5,d20.5)
end
c
c
subroutine cntrac(cint,norb,n,nt,nro,buff)
c
c transformation of integrals in mo-basis. only diagonalelements
c are calculated
c
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
common /contrc/ chomo,ncont
common/ioind/isl1, isl2
common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
dimension cint(nomx,1),buff(lbuff)
rewind 2
do 10 m=1,nt
do 10 i=1,norb
10 cint(i,m)=0.0d0
nend=ncont*norb
ix=lbuff
do 50 k=1,ncont
do 50 l=1,k
do 40 m=1,nt
ix=ix+1
if(ix.le.lbuff)go to 20
ix=1
read(2)buff
20 do 30 i=1,norb
kk=ncont*(i-1)+k
ll=ncont*(i-1)+l
tmp=y(kk)*y(ll)
if(k.ne.l) tmp=tmp+tmp
30 cint(i,m)=buff(ix)*tmp+cint(i,m)
40 continue
50 continue
return
end
c
c
subroutine cntrcm(cint,norb,non,buff,frocc,rnabcd)
c
c mulliken and roby population analysis
c
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
integer centr
common/scfinf/dmh(nomx),roc(nomx)
common /big/ y(nomx*nomx),su(nomx*nomx),uv(nomx*nomx)
common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
common/muli/vlist(nocmx,4),ianr(nomx),centr(nocmx),itnr(nomx)
common /contrc/ chomo,ncont
common/ioind/isl1, isl2
common/supp/nsupp
data ifile/4/,jfile/7/,kfile/9/,tol/0.d0/
dimension cint(2),buff(lbuff),rnmy(20),frocc(2),mor(20)
c number of aos : 250
c number of atoms : 16
c dimension of ovpop : 16*(16+1)/2
c dimension of nat and jatnr : 16+1
dimension ovpop(136),mos(250),iatnr(16),nat(17),rmopo(250,250)
dimension rnab(15,16),jatnr(17),rnabc(14,15,16),rnamo(250,16)
dimension rnabcd(1),katnr(16)
c rnabcd(1) equivalenced to buff(16+1) in main program
equivalence (uv(1),rmopo(1,1)), (rnabc(1,1,1),rnamo(1,1))
tri=1.d0/3.d0
four=1.d0/4.d0
ipr=ncont*(ncont+1)/2
ncmax=ncont*ncont
maxuv=nomx*nomx
idmax=maxuv/2
nocm1=nocmx-1
nocm2=nocmx-2
nocp1=nocmx+1
c
c mulliken population analysis
c
c overlap * mos, written on kfile (for roby analysis)
call hm(buff,y,su,norb,ncont)
momax=norb*ncont
rewind kfile
write(kfile) (su(ir),ir=1,momax)
rewind kfile
c
c calculation of density , written on jfile
do 205 i=1,ipr
205 su(i)=0.d0
do 210 k=1,norb
ij=0
ist=(k-1)*ncont
do 220 i=1,ncont
iy=ist+i
do 220 j=1,i
ij=ij+1
jy=ist+j
220 su(ij)=su(ij)+y(iy)*y(jy)*frocc(k)
210 continue
rewind jfile
write(jfile) (su(ir),ir=1,ipr)
rewind jfile
c
c
do 10 i=1,ncont
10 cint(i)=0.0d0
ix=0
iend=non*(non+1)/2
ij=0
do 15 i=1,iend
ovpop(i)=0.d0
do 15 j=1,norb
15 rmopo(i,j)=0.d0
c calculation of overlap charges
do 50 k=1,ncont
do 50 l=1,k
ij=ij+1
ix=ix+1
tmp=0.d0
do 30 i=1,norb
kk=ncont*(i-1)+k
ll=ncont*(i-1)+l
ka=ianr(k)
la=ianr(l)
kla=(ka*(ka-1))/2+la
tmpb=y(kk)*y(ll)*frocc(i)*buff(ix)
if(ka.ne.la) tmpb=tmpb*2.0d0
rmopo(kla,i)=rmopo(kla,i)+tmpb
if(k.eq.l) goto 30
if(ka.eq.la) rmopo(kla,i)=rmopo(kla,i)+tmpb
30 continue
tmpb=su(ij)*buff(ix)
cint(k)=cint(k)+tmpb
if(k.ne.l) cint(l)=cint(l)+tmpb
50 continue
ren=0.d0
do 150 ii=1,ncont
150 ren=ren+cint(ii)
write(6,1000)
c calculation of occupation numbers
do 70 j=1,non
do 71 k=1,20
71 rnmy(k)=0.d0
rnmy2=0.d0
rnmy3=0.d0
rnmy4=0.d0
rnm=0.d0
do 60 i=1,ncont
if(ianr(i).ne.j) goto 60
itp=itnr(i)
rnmy(itp)=rnmy(itp)+cint(i)
60 continue
rnmy2=rnmy2+rnmy(2)+rnmy(3)+rnmy(4)
do 80 k=5,10
80 rnmy3=rnmy3+rnmy(k)
do 90 k=11,20
90 rnmy4=rnmy4+rnmy(k)
rnm=rnmy(1)+rnmy2+rnmy3+rnmy4
chg=vlist(j,4)-rnm
write(6,2000)
write(6,3000) j,centr(j),rnmy(1),rnmy2,rnmy3,rnmy4,rnm,chg
if(rnmy2.ne.0.d0) write(6,4000) (rnmy(ii),ii=2,4)
if(rnmy3.ne.0.d0) write(6,4001) (rnmy(ii),ii=5,10)
if(rnmy4.ne.0.d0) write(6,4002) (rnmy(ii),ii=11,20)
70 continue
write(6,5000) ren
write(6,6000)
fakt=1.0d0
do 115 i=1,iend
do 115 j=1,norb
115 ovpop(i)=ovpop(i)+rmopo(i,j)
call outpak(ovpop,non,1,1)
if(nsupp.ne.5) goto 218
ii=0
write(6,6400)
do 110 i=1,non
do 110 j=1,i
write(6,6500) i,centr(i),j,centr(j)
ii=ii+1
write(6,1001) (rmopo(ii,k),k=1,norb)
110 continue
c
c roby population analysis
c
218 continue
write(6,8000)
nnc=ncont
c
c calculate blocks of atoms
itmp=1
ij=0
do 400 i=2,ncont
if(ianr(i).ne.ianr(i-1)) goto 405
itmp=itmp+1
goto 400
405 ij=ij+1
iatnr(ij)=itmp
itmp=1
400 continue
iatnr(ij+1)=itmp
c
ndim=0
do 733 i=1,non
ia=iatnr(i)
733 ndim=ndim+(ia*ia+ia)/2
c
c atom blocked overlap matrix on su
call sort(buff,su,iatnr,non,0,ncmax)
c
c calculation of s**1/2 and s**-1/2 (diagonal)
ist=1
jst=1
ij=0
ll=0
do 530 i=1,non
ia=iatnr(i)
do 540 j=1,ia
do 540 k=1,ia
ij=ij+1
uv(ij)=0.d0
if(j.eq.k) uv(ij)=1.d0
540 continue
call erduw(su(jst),uv(ist),ia,1.d-12)
ii=jst-1
do 560 l=1,ia
ii=ii+l
ll=ll+1
cint(ll)=sqrt(su(ii))
560 dmh(ll)=1.d0/cint(ll)
ist=ist+ia*ia
jst=jst+(ia*ia+ia)/2
530 continue
itest=ist-1
if(itest.le.idmax) goto 94
write(6,7100) itest,idmax
7100 format(1x,'itest = ',i8,' greater than idmax = ',i8,
1 ' in subroutine cntrcm',/,' enlarge nomx in main program')
stop
94 continue
c
c calculation of s**1/2 (atom blocked)
call uduk(uv,cint,su,iatnr,non,ndim)
rewind jfile
read(jfile) (y(ii),ii=1,ipr)
rewind jfile
c
c atom blocked density matrix on uv
call sort(y,uv(idmax+1),iatnr,non,0,ncmax)
c s**1/2 * density * s**1/2
call h1h2h1(su,uv(idmax+1),cint,y,iatnr,non,ipr)
c
c calculation of s**-1/2 (atom blocked)
call uduk(uv,dmh,su,iatnr,non,ndim)
c
c calculation of maos (written on uv)
ij=0
ist=1
jst=1
do 570 i=1,non
ia=iatnr(i)
do 580 j=1,ia
iuv1=jst-1+(j-1)*ia
iuv2=jst-1+j
do 580 k=1,j
ij=ij+1
iuv=iuv1+k
uv(iuv)=su(ij)
iuv=iuv2+(k-1)*ia
580 uv(iuv)=su(ij)
call erduw(y(ist),uv(jst),ia,1.d-12)
ist=ist+(ia*ia+ia)/2
jst=jst+ia*ia
570 continue
c full mao matrix on y
call sort(uv,y,iatnr,non,1,ncmax)
c calculation of s * mao
call hm(buff,y,uv,ncont,ncont)
c read in density and calculate (mao*s)# * d * (mao*s)
rewind jfile
read(jfile) (buff(ii),ii=1,ipr)
rewind jfile
call ukhu(uv,buff,cint,su,ncont,ncont)
rewind ifile
write(ifile) (su(ir),ir=1,ncmax)
rewind ifile
c
ii=ncont+1
ij=1
do 585 i=1,ncont
roc(i)=su(ij)
585 ij=ij+ii
c mao * (s*mao), written on jfile
call mkm(y,uv,su,ncont,ncont,ncont)
rewind jfile
write(jfile) (su(ir),ir=1,ncmax)
rewind jfile
c
do 649 i=1,ncont
649 mos(i)=0
c
c roby occupation numbers
nnc=0
ii=0
ist=0
do 600 i=1,non
jatnr(i)=iatnr(i)
ia=iatnr(i)
rat=0.d0
do 610 j=1,ia
ii=ii+1
610 rat=rat+roc(ii)
write(6,6700) i,centr(i),rat
write(6,6800) (roc(ir),ir=ist+1,ii)
c
c selection of maos by user
rat=0.d0
read(5,1500) nmao,(mor(im),im=1,nmao)
if(mor(1).eq.0) goto 851
do 854 im=1,nmao
kk=mor(im)
mos(kk)=kk
cint(im)=roc(kk)
854 rat=rat+roc(kk)
iatnr(i)=nmao
nnc=nnc+nmao
goto 855
851 iatnr(i)=nmao
c automatic selection of maos according to occupation numbers
if(nmao.eq.0) write(6,6900) nmao,i
if(nmao.eq.0) goto 600
nnc=nnc+nmao
iij=0
do 850 k=1,nmao
rpr=0.d0
do 860 l=1,ia
if(roc(ist+l).le.rpr) goto 860
kk=ist+l
rpr=roc(kk)
860 continue
iij=iij+1
mos(kk)=kk
rat=rat+rpr
cint(iij)=rpr
roc(kk)=-rpr
850 continue
855 rnmy(i)=rat
write(6,6900) nmao,i,rat
write(6,7500) (cint(ir),ir=1,nmao)
600 ist=ist+ia
c
do 871 i=1,non
871 ianr(i)=jatnr(i)-iatnr(i)
min=1000
do 872 i=1,non
if(ianr(i).lt.min) min=ianr(i)
872 continue
c hypervalence contributions
call hypval(jatnr,mos,tol,ipr,ifile,non,ncmax,ren,
1 jfile,rtot)
c
c
nnc2=nnc*nnc
mm=0
do 870 m=1,ncont
if(mos(m).eq.0) goto 870
mm=mm+1
mos(mm)=mos(m)
roc(mm)=abs(roc(m))
870 continue
c
c calculate trace (mao*b*mao * s**-1)
c calculation of unassigned charge
rewind ifile
read(ifile) (uv(iw),iw=1,ncmax)
rewind ifile
rewind jfile
read(jfile) (su(iw),iw=1,ncmax)
rewind jfile
c order according to selected maos
call order(uv,su,mos,nnc,ncont)
c ordered overlap(mao basis) written on jfile
rewind jfile
write(jfile) (su(iw),iw=1,nnc2)
rewind jfile
c
if(nsupp.ne.4) goto 655
write(6,9100)
call prtblk(su,nnc,nnc,'mao ','mao ',mos)
write(6,9200)
655 call osinv1(su,nnc,d,tol,itnr,ianr)
c
rna=0.d0
do 910 i=1,nnc
a=0.d0
i1=i
i20=(i-1)*nnc
do 920 j=1,nnc
a=a+su(i1)*uv(i20+j)
920 i1=i1+nnc
rna=rna+a
910 continue
uach=ren-rna
write(6,7700) uach
rtot=(uach-rtot)/non
do 500 i=1,non
500 dmh(i)=dmh(i)+rtot
write(6,7800) rtot
c
c selection and reordering of mao vector y
imo=0
ie=0
ij=0
nnc=0
nst=1
nat(1)=1
do 650 l=1,non
nmao=iatnr(l)
nat(l+1)=nat(l)+nmao
if(nmao.eq.0) goto 650
do 660 i=1,nmao
imo=imo+1
imos=mos(imo)-1
ist=imos*ncont+1
ie=ist+ncont-1
do 680 k=ist,ie
ij=ij+1
680 y(ij)=y(k)
660 continue
write(6,7600) l,centr(l)
call prtblk(y(nst),ncont,nmao,' ao ','mao ',mos(nnc+1))
nst=nst+nmao*ncont
nnc=nnc+nmao
650 continue
c write out maos for plot
ndimw=nnc*ncont
write(lomao,8900) (alabel(i1),i1=1,19)
write(lomao,8900) (blabel(i1),i1=1,19)
write(lomao,9000) nnc,ncont,ndimw
write(lomao,8900) (fmt(i1),i1=1,4)
write(lomao,fmt) (y(i),i=1,ndimw)
c
c calculation of occupation numbers (mo-contributions)
rewind kfile
read(kfile) (uv(ir),ir=1,momax)
rewind kfile
c (maos)# * (sc) = su
call mkm(y,uv,su,nnc,ncont,norb)
if(nsupp.eq.5) write(6,8650)
do 960 m=1,non
rat=0.d0
mm=0
do 970 j=1,norb
ind=(j-1)*nnc
do 980 k=nat(m),nat(m+1)-1
mm=mm+1
980 uv(mm)=su(ind+k)
970 continue
ndim=iatnr(m)
i0=-ndim
do 950 i=1,norb
v=0.d0
i0=i0+ndim
do 955 j=1,ndim
955 v=v+uv(i0+j)*uv(i0+j)
rnamo(i,m)=v*frocc(i)
rat=rat+rnamo(i,m)
950 continue
if(nsupp.eq.5) write(6,8655) m,centr(m),(rnamo(ir,m),ir=1,norb)
960 continue
c
c two center shared electron numbers
write(6,8350)
ncen=1
ishare=2
rewind jfile
read(jfile) (uv(ir),ir=1,nnc2)
rewind jfile
do 715 i=1,non
715 buff(i)=0.d0
nat(non+1)=nnc+1
c
do 700 i=1,non-1
ij=0
if(iatnr(i).eq.0) goto 700
do 710 ii=nat(i),nat(i+1)-1
ij=ij+1
710 mos(ij)=ii
do 720 j=i+1,non
if(iatnr(j).eq.0) goto 720
ik=ij
do 730 jj=nat(j),nat(j+1)-1
ik=ik+1
730 mos(ik)=jj
itest=ik*ik
if(itest.gt.idmax) goto 95
c
call sen(su,uv,y(1),y(idmax),mos,cint,rna,ik,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc)
rnab(i,j)=rna
write(6,7900) i,centr(i),j,centr(j)
rna=0.d0
do 810 im=1,norb
roc(im)=rnamo(im,i)+rnamo(im,j)-roc(im)
810 rna=rna+roc(im)
write(6,8100) (roc(ir),ir=1,norb)
buff(i)=buff(i)+rna*0.5d0
buff(j)=buff(j)+rna*0.5d0
write(6,8200) rna
720 continue
700 continue
c
c three center shared electron numbers
c
if(non.le.3) goto 515
ncen=-1
ishare=3
write(6,8400)
do 740 i=1,non-2
ij=0
if(iatnr(i).eq.0) goto 740
do 745 ii=nat(i),nat(i+1)-1
ij=ij+1
745 mos(ij)=ii
do 750 j=i+1,non-1
if(iatnr(j).eq.0) goto 750
ik=ij
do 755 jj=nat(j),nat(j+1)-1
ik=ik+1
755 mos(ik)=jj
do 760 k=j+1,non
if(iatnr(k).eq.0) goto 760
il=ik
do 765 kk=nat(k),nat(k+1)-1
il=il+1
765 mos(il)=kk
itest=il*il
if(itest.gt.idmax) goto 95
c
call sen(su,uv,y(1),y(idmax),mos,cint,rna,il,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc)
rnabc(i,j,k)=rna
rna=rnmy(i)+rnmy(j)+rnmy(k)-rnab(i,j)-rnab(i,k)-rnab(j,k)+rna
buff(i)=buff(i)-tri*rna
buff(j)=buff(j)-tri*rna
buff(k)=buff(k)-tri*rna
if(abs(rna).ge.0.01d0) write(6,8500) i,j,k,rna
760 continue
750 continue
740 continue
c
c four center shared electron numbers
c
if(non.le.4) goto 515
ncen=1
ishare=4
write(6,8450)
do 770 i=1,non-3
ij=0
if(iatnr(i).eq.0) goto 770
do 775 ii=nat(i),nat(i+1)-1
ij=ij+1
775 mos(ij)=ii
do 780 j=i+1,non-2
if(iatnr(j).eq.0) goto 780
ik=ij
do 785 jj=nat(j),nat(j+1)-1
ik=ik+1
785 mos(ik)=jj
do 790 k=j+1,non-1
if(iatnr(k).eq.0) goto 790
il=ik
do 795 kk=nat(k),nat(k+1)-1
il=il+1
795 mos(il)=kk
do 800 l=k+1,non
if(iatnr(l).eq.0) goto 800
im=il
do 805 ll=nat(l),nat(l+1)-1
im=im+1
805 mos(im)=ll
itest=im*im
if(itest.gt.idmax) goto 95
c
call sen(su,uv,y(1),y(idmax),mos,cint,rna,im,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc)
index1=j*(j-3)/2+i
index2=l*(l-3)/2+k
index=index2*(index2-3)/2+index1
rnabcd(index)=rna
rna=rnmy(i)+rnmy(j)+rnmy(k)+rnmy(l)-rnab(i,j)-rnab(i,k)-rnab(i,l)
1 -rnab(j,k)-rnab(j,l)-rnab(k,l)+rnabc(i,j,k)+rnabc(i,j,l)
2 +rnabc(i,k,l)+rnabc(j,k,l)-rna
buff(i)=buff(i)+four*rna
buff(j)=buff(j)+four*rna
buff(k)=buff(k)+four*rna
buff(l)=buff(l)+four*rna
if(abs(rna).ge.0.002d0) write(6,8550) i,j,k,l,rna
800 continue
790 continue
780 continue
770 continue
515 read(5,6100) nm,(iatnr(i),i=1,nm)
if(nm.eq.0) goto 933
call spemu(su,uv,y(1),y(idmax),mos,cint,rna,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc,iatnr,rnmy,rnab,rnabc,rnabcd,nm,nat,idmax,
2 jatnr,katnr,nocm1,nocm2,nocp1)
c
933 continue
c
c calculation of atomic charges
c chomo: charge of molecule, cmult: multicenter contribution
sum=0.d0
do 786 i=1,non
rnmy(i)=vlist(i,4)-rnmy(i)+buff(i)-dmh(i)
sum=sum+rnmy(i)
786 continue
sum=sum-chomo
cmult=sum*ncen
if(non.eq.3) write(6,8840) cmult
if(non.eq.4) write(6,8850) cmult
if(non.ge.5) write(6,8800) cmult
write(6,8600)
sum=sum/non
do 787 i=1,non
q=rnmy(i)-sum
write(6,8700) i,centr(i),q
787 continue
c
c
1000 format(///,' mulliken population analysis',//)
1001 format(1x,12f10.6)
1500 format(i10,20i3)
2000 format(//,' atom nr. symbol s-occupation p-occupation ',
1 'd-occupation f-occupation total charge ')
3000 format(1x,i5,4x,a4,2x,6(f12.5,2x))
4000 format(/,1x,'n(px) =',f9.5,2x,'n(py) =',f9.5,2x,'n(pz) =',f9.5)
4001 format(1x,'n(dxx) =',f9.5,2x,'n(dyy) =',f9.5,2x,'n(dzz) =',f9.5,
1 2x,'n(dxy) =',f9.5,2x,'n(dxz) =',f9.5,2x,'n(dyz) =',f9.5)
4002 format(1x,'n(fxxx)=',f9.5,2x,'n(fyyy)=',f9.5,2x,'n(fzzz)=',f9.5,
1 2x,'n(fxxy)=',f9.5,2x,'n(fxxz)=',f9.5,2x,'n(fyyx)=',f9.5,/,
2 1x,'n(fyyz)=',f9.5,2x,'n(fzzx)=',f9.5,2x,'n(fzzy)=',f9.5,
3 2x,'n(fxyz)=',f9.5)
5000 format(//,' total number of electrons = ',f15.8)
6000 format(///,' mulliken overlap population ',/)
6100 format(20i3)
6400 format(//,' overlap population per mo',/)
6500 format(/,' atom ',i5,a4,' / ','atom ',i5,a4)
6700 format(//,' atom nr ',i5,2x,a4,' occupation',f10.6)
6800 format(1x,'occupation of maos',/,(10f10.6))
6900 format(1x,i5,' maos selected for atom ',i4,' occupation ',f10.6)
7500 format(1x,'occupation of maos ',/,(10f10.6))
7600 format(///,' selected maos for atom ',i3,a4)
7700 format(/,10x,'unassigned charge = ',f10.4,/)
7800 format(1x,'rest charge per atom = ',f10.6)
7900 format(//,' atom ',i3,a4,' / atom ',i3,a4)
8000 format(///,' roby population analysis with modified atomic',
1 ' orbitals ')
8100 format(1x,'contributions of mos to shared electron number',
1 /,(10f10.5))
8200 format(1x,' shared electron number = ',f10.4)
8300 format(1x,'n( 1 2 3 ) = ',f10.5)
8350 format(///,' two center shared electron numbers '/)
8400 format(///,' three center shared electron numbers ge 0.01'/)
8450 format(///,' four center shared electron numbers ge 0.002'/)
8500 format(1x,'n(',3i3,') = ',f10.4,f10.6)
8550 format(1x,'n(',4i3,') = ',f10.4)
8600 format(//,' atomic charges with multicenter corrections',
1 /,' atom charge')
8650 format(//,' contributions of mos to occupation numbers')
8655 format(/,' atom nr ',i5,2x,a4,/,(12f10.6))
8700 format(1x,i5,a4,f10.4)
8800 format(//,' multicenter rest contribution = ',f10.6,
1 ' (correct sign for five center terms)')
8840 format(//,' three center shared electron number = ',f10.6)
8850 format(//,' four center shared electron number = ',f10.6)
8900 format(19a4)
9000 format(20i4)
9100 format(///' s matrix in mao basis'/)
9200 format(//)
9500 format(a1,i4,14i5)
9550 format(1x,2i5)
return
c
95 write(6,7200) itest,idmax,ishare
7200 format(1x,'itest = ',i8,' greater than idmax = ',i8,
1 ' in subroutine cntrcm for',i3,'-center sens',/,
2 ' enlarge all arrays of dimension nomx')
stop
end
c
c
subroutine spemu(su,uv,ssu,suv,mos,cint,rna,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc,iatnr,rnmy,rnab,rnabc,rnabcd,nm,nat,idmax,
2 jatnr,katnr,nocm1,nocm2,nocp1)
implicit real*8(a-h,o-z)
dimension su(1),uv(1),ssu(1),suv(1),mos(1),cint(1),itnr(1),ianr(1)
dimension roc(1),frocc(1),iatnr(1),rnmy(1),rnab(nocm1,1),nat(1)
dimension rnabc(nocm2,nocm1,1),rnabcd(1),jatnr(nocp1),katnr(1)
c
iend=nm
do 310 j=1,nm-1
max=iatnr(1)
ind=1
do 320 i=2,iend
if(iatnr(i).lt.max) goto 320
max=iatnr(i)
ind=i
320 continue
jj=nm-j+1
iend=jj-1
iatnr(ind)=iatnr(jj)
310 iatnr(jj)=max
write(6,2000) (iatnr(i),i=1,nm)
c
ij=0
do 50 i=1,nm
i1=iatnr(i)
do 60 ii=nat(i1),nat(i1+1)-1
ij=ij+1
60 mos(ij)=ii
50 continue
itest=ij*ij
if(itest.gt.idmax) goto 99
call sen(su,uv,ssu,suv,mos,cint,rna,ij,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc)
tmult=0.d0
do 100 i=1,nm
i1=iatnr(i)
tmult=tmult+rnmy(i1)
100 continue
tmult=tmult-rna
write(6,3500) tmult
if(nm.le.4) return
ipha=(-1)**(nm-1)
sum=ipha*rna
c
lmax=int(4.e+09)
do 10 i=1,nm
i1=iatnr(i)
sum=sum+rnmy(i1)
if(i.eq.1) goto 10
do 20 j=1,i-1
j1=iatnr(j)
sum=sum-rnab(j1,i1)
if(j.eq.1) goto 20
do 30 k=1,j-1
k1=iatnr(k)
sum=sum+rnabc(k1,j1,i1)
if(k.eq.1) goto 30
do 40 l=1,k-1
l1=iatnr(l)
index1=k1*(k1-3)/2+l1
index2=i1*(i1-3)/2+j1
index=index2*(index2-3)/2+index1
40 sum=sum-rnabcd(index)
30 continue
20 continue
10 continue
if(nm.eq.5) goto 410
c
do 110 i=5,nm-1
do 115 ii=1,i
katnr(ii)=ii
115 jatnr(ii)=iatnr(ii)
ij=0
do 81 k=1,i
i1=jatnr(k)
do 91 ii=nat(i1),nat(i1+1)-1
ij=ij+1
91 mos(ij)=ii
81 continue
itest=ij*ij
if(itest.gt.idmax) goto 99
call sen(su,uv,ssu,suv,mos,cint,rna,ij,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc)
ipha=(-1)**(i-1)
sum=sum+ipha*rna
c
jatnr(i+1)=1000
jst=i
do 140 l=1,lmax
do 120 j=jst,i
do 130 k=jst+1,nm
if(iatnr(k).le.jatnr(j)) goto 130
if(iatnr(k).ge.jatnr(j+1)) goto 130
jatnr(j)=iatnr(k)
katnr(j)=k
ij=0
do 80 ll=1,i
i1=jatnr(ll)
do 90 ii=nat(i1),nat(i1+1)-1
ij=ij+1
90 mos(ij)=ii
80 continue
itest=ij*ij
if(itest.gt.idmax) goto 99
call sen(su,uv,ssu,suv,mos,cint,rna,ij,d,tol,itnr,ianr,
1 norb,nnc,roc,frocc)
ipha=(-1)**(i-1)
sum=sum+ipha*rna
130 continue
120 continue
ki=0
do 145 kk=1,i
if(jatnr(kk).ne.iatnr(nm-i+kk)) goto 145
ki=ki+1
if(ki.eq.i) goto 110
if(ki.eq.1) jst=kk-1
kkk=katnr(kk-1)+1
if(ki.eq.1) kkk=katnr(kk-1)+2
jatnr(kk)=iatnr(kkk)
katnr(kk)=kkk
145 continue
140 continue
110 continue
c
410 write(6,3000) sum
return
c
1000 format(a1,i2,20i3)
2000 format(///,' special multi center contribution for the atoms',/,
1 (20i5))
3000 format(1x,'multi center shared electron number =',f15.6,/)
3500 format(1x,'total contribution =',f15.6)
c
99 write(6,4000) itest,idmax
4000 format(1x,'itest = ',i8,' greater than idmax = ',i8,
1 ' in subroutine spemu for special multi center contributions',/,
2 ' enlarge all arrays of dimension nomx')
return
end
c
c
subroutine sen(su,uv,ssu,suv,mos,cint,rna,idim,d,tol,itnr,
1 ianr,norb,nnc,roc,frocc)
c
c calculation of shared electron numbers
c
implicit real*8 (a-h,o-z)
dimension su(1),uv(1),ssu(1),suv(1),mos(1),cint(1),itnr(1)
dimension ianr(1),roc(1),frocc(1)
c
c sort (mao)#*s*c on ssu
ij=0
do 10 i=1,norb
ind0=(i-1)*nnc
do 20 j=1,idim
ij=ij+1
ind=ind0+mos(j)
20 ssu(ij)=su(ind)
10 continue
c
c sort s(mao basis) on suv
ij=0
do 30 j=1,idim
ind0=(mos(j)-1)*nnc
do 40 i=1,idim
ind=ind0+mos(i)
ij=ij+1
40 suv(ij)=uv(ind)
30 continue
c
call osinv1(suv,idim,d,tol,itnr,ianr)
ij=0
do 60 i=1,idim
ind=(i-1)*idim
do 60 j=1,i
ij=ij+1
suv(ij)=suv(ind+j)
60 continue
call ukhud(ssu,suv,cint,roc,idim,norb)
c
rna=0.d0
do 50 i=1,norb
roc(i)=roc(i)*frocc(i)
50 rna=rna+roc(i)
return
end
c
c
subroutine ukhud(a,b,c,d,idim1,idim2)
c
c u# * triangle matrix * u , diagonalelements only
implicit real*8(a-h,o-z)
dimension a(1),b(1),c(1),d(1)
kl=0
ij=0
i10=-idim1
do 40 j=1,idim2
i10=i10+idim1
ij=0
do 50 k=1,idim1
cc=0.d0
ij=ij+1
i20=(k*k-k)/2
do 60 l=1,k
60 cc=cc+a(i10+l)*b(i20+l)
if(k.eq.idim1) goto 55
i20=k*(k+3)/2
do 70 l=k+1,idim1
cc=cc+a(i10+l)*b(i20)
70 i20=i20+l
55 c(ij)=cc
50 continue
m1=0
m2=(j-1)*idim1
kl=kl+1
dd=0.0d0
do 90 n=1,idim1
90 dd=dd+c(m1+n)*a(m2+n)
d(kl)=dd
40 continue
return
end
c
subroutine hypval(jatnr,mos,tol,ipr,ifile,non,
1 ncmax,ren,jfile,rtot)
c
c calculation of hypervalence populations
c
implicit real*8 (a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
integer centr
dimension jatnr(1),mos(1)
common/scfinf/dmh(nomx),roc(nomx)
common /big/ y(nomx*nomx),su(nomx*nomx),uv(nomx*nomx)
common/muli/vlist(nocmx,4),ianr(nomx),centr(nocmx),itnr(nomx)
common /contrc/ chomo,ncont
iend=0
rtot=0.d0
write(6,1000)
do 100 ia=1,non
ist=iend+1
iend=iend+jatnr(ia)
nnc=ist-1
if(nnc.eq.0) goto 111
do 50 i=1,nnc
50 itnr(i)=i
111 continue
do 110 ii=ist,iend
if(mos(ii).eq.0) goto 110
nnc=nnc+1
itnr(nnc)=ii
110 continue
if(ia.eq.non) goto 121
do 120 ii=iend+1,ncont
nnc=nnc+1
itnr(nnc)=ii
120 continue
c
121 rewind ifile
read(ifile) (uv(iw),iw=1,ncmax)
rewind ifile
rewind jfile
read(jfile) (su(iw),iw=1,ncmax)
rewind jfile
call order(uv,su,itnr,nnc,ncont)
call osinv1(su,nnc,d,tol,itnr,ianr)
c
c calculate trace(mao*s*d*s*mao * s**-1)
rna=0.d0
do 910 i=1,nnc
a=0.d0
i1=i
i20=(i-1)*nnc
do 920 j=1,nnc
a=a+su(i1)*uv(i20+j)
920 i1=i1+nnc
rna=rna+a
910 continue
rna=ren-rna
dmh(ia)=rna
rtot=rtot+rna
write(6,2000) ia,centr(ia),rna
100 continue
c
1000 format(/' hypervalence contributions'//' unshared atomic',
1 ' populations')
2000 format(1x,'atom nr.',i3,1x,a4,5x,f10.6)
return
end
c
subroutine order(a,b,iord,idims,idimb)
implicit real*8(a-h,o-z)
dimension a(1),b(1),iord(1)
c sort matrix a an b according to vector iord
ij=0
do 10 j=1,idims
ind2=(iord(j)-1)*idimb
do 20 i=1,idims
ind=ind2+iord(i)
ij=ij+1
a(ij)=a(ind)
20 b(ij)=b(ind)
10 continue
return
end
c
c
subroutine hm(a,b,c,norb,ncont)
c triangle matrix(a) * full matrix(b)
implicit real*8(a-h,o-z)
dimension a(1),b(1),c(1)
ij=0
do 20 i=1,norb
ibst=(i-1)*ncont
do 30 j=1,ncont
ct=0.d0
ij=ij+1
iast=j*(j-1)/2
do 40 k=1,j
40 ct=ct+a(iast+k)*b(ibst+k)
if(j.eq.ncont) goto 29
iast=j*(j+3)/2
do 50 k=j+1,ncont
ct=ct+a(iast)*b(ibst+k)
50 iast=iast+k
29 c(ij)=ct
30 continue
20 continue
return
end
c
subroutine mkm(a,b,c,idima,idimab,idimb)
implicit real*8(a-h,o-z)
dimension a(1),b(1),c(1)
c (full matrix)# * full matrix
ij=0
do 10 i=1,idimb
iast=0
ibst=(i-1)*idimab
do 20 j=1,idima
ij=ij+1
v=0.d0
do 30 k=1,idimab
30 v=v+a(iast+k)*b(ibst+k)
c(ij)=v
20 iast=iast+idimab
10 continue
return
end
c
c
c
c
subroutine sort(a,b,nr,non,iswi,ncmax)
c sorting of matrices in atom blocked matrices
implicit real*8(a-h,o-z)
dimension a(1),b(1),nr(1)
c iswi=0: full matrix in atom blocked matrix
ij=0
if(iswi.eq.1) goto 30
jend=0
do 10 i=1,non
jst=jend+1
jend=jend+nr(i)
do 10 j=jst,jend
kkst=(j*j-j)/2
do 10 k=jst,j
ij=ij+1
kl=kkst+k
b(ij)=a(kl)
10 continue
return
c iswi=1: atom blocked matrix in full matrix
30 iofs=0
ndim=0
do 35 i=1,non
35 ndim=ndim+nr(i)
do 38 i=1,ncmax
38 b(i)=0.d0
iaa=0
kl=0
do 40 i=1,non
ia=nr(i)
iaa=iaa+ia
do 50 j=1,ia
do 60 k=1,ia
ij=ij+1
kl=kl+1
60 b(kl)=a(ij)
kl=kl+ndim-iaa+iofs
50 continue
kl=kl+ia
iofs=iofs+ia
40 continue
return
end
c
c
c
c
c
c
subroutine uduk(a,b,c,nr,non,ndim)
c u# * vector * u
implicit real*8(a-h,o-z)
dimension a(1),b(1),c(1),nr(1)
do 10 i=1,ndim
10 c(i)=0.d0
ii=0
ijst=0
iist=0
do 20 l=1,non
ia=nr(l)
do 30 i=1,ia
ii=ii+1
ij=ijst
ist=iist+(i-1)*ia
do 40 j=1,ia
iu=ist+j
do 40 k=1,j
ij=ij+1
ju=ist+k
40 c(ij)=c(ij)+a(iu)*b(ii)*a(ju)
30 continue
ijst=ij
iist=iu
20 continue
return
end
c
c
subroutine ukhu(a,b,c,d,ncont,nc)
c u# * triangle matrix * u
implicit real*8(a-h,o-z)
dimension a(1),b(1),c(1),d(1)
i10=-ncont
do 40 j=1,nc
i10=i10+ncont
ij=0
i11=i10+1
do 50 k=1,ncont
ij=ij+1
i20=(k*k-k)/2
i21=i20+1
cij = 0.0d0
eij = 0.0d0
if (k. eq. 1) go to 110
do 60 l = 1, k-1, 2
cij = cij + a(i10+l) * b(i20+l)
60 eij = eij + a(i11+l) * b(i21+l)
if (mod(k,2) .eq. 0) go to 120
110 continue
cij = cij + a(i10+k) * b(i20+k)
120 continue
if(k.eq.ncont) goto 49
loff = (k+3) * k / 2
if (k+1 .eq .ncont) go to 130
do 70 l=k+1,ncont-1,2
cij=cij+a(i10+l) * b(loff)
eij=eij+a(i11+l) * b(loff+l)
70 loff = loff + l + l + 1
if (mod(ncont-k,2) .eq. 0) go to 49
130 cij=cij+a(i10+ncont) * b(loff)
49 c(ij) = cij+eij
50 continue
m2=0
do 80 m=1,nc
kl=(m-1)*nc+j
dkl=0.d0
fkl=0.0d0
m21 = m2 + 1
if (ncont .eq. 1) go to 140
do 90 n=1,ncont-1,2
dkl=dkl+c(n)*a(m2+n)
90 fkl=fkl+c(n+1)*a(m21+n)
if(mod(ncont,2) .eq. 0) go to 150
140 dkl = dkl + c(ncont) * a(m2+ncont)
150 continue
d(kl)=dkl+fkl
m2=m2+ncont
80 continue
40 continue
return
end
c
c
subroutine osinv1(a,n,d,tol,l,m)
c inversion of a full matrix
implicit real*8 (a-h,o-z)
dimension a(1),m(1),l(1)
if(n.ne.1) goto 5
a(1)=1.0/a(1)
return
5 d=1.d0
nk=-n
do 180 k=1,n
nk=nk+n
l(k)=k
m(k)=k
kk=nk+k
biga=a(kk)
do 20 j=k,n
iz=n*(j-1)
do 20 i=k,n
ij=iz+i
if(abs(biga)-abs(a(ij))) 10,20,20
10 biga=a(ij)
l(k)=i
m(k)=j
20 continue
j=l(k)
if(j-k) 50,50,30
30 ki=k-n
do 40 i=1,n
ki=ki+n
holo=-a(ki)
ji=ki-k+j
a(ki)=a(ji)
40 a(ji)=holo
50 i=m(k)
if(i-k) 80,80,60
60 jp=n*(i-1)
do 70 j=1,n
jk=nk+j
ji=jp+j
holo=-a(jk)
a(jk)=a(ji)
70 a(ji)=holo
80 if(abs(biga)-tol) 90,100,100
90 d=0.d0
return
100 continue
biginv = -1.0d0 / biga
do 120 i=1,n
a(nk+i)=a(nk+i)*biginv
120 continue
a(nk+k) = a(nk+k) * (-biga)
joff = 0
do 155 j=1,n
if(j. eq. k) go to 154
ajk = a(joff+k)
if (k .eq. 1) go to 152
do 151 i=1,k-1
151 a(joff+i) = a(joff+i) + a(nk+i) * ajk
152 continue
if (k .eq. n) go to 154
do 153 i=k+1,n
153 a(joff+i) = a(joff+i) + a(nk+i) * ajk
154 continue
joff = joff + n
155 continue
kj=k-n
biginv=-biginv
do 170 j=1,n
kj=kj+n
a(kj)=a(kj)*biginv
170 continue
kadd = (k-1)*n + k
a(kadd) = a(kadd) * biga
d=d*biga
a(kk)=biginv
180 continue
k=n
190 k=k-1
if(k) 260,260,200
200 i=l(k)
if(i-k) 230,230,210
210 jq=n*(k-1)
jr=n*(i-1)
do 220 j=1,n
jk=jq+j
holo=a(jk)
ji=jr+j
a(jk)=-a(ji)
220 a(ji)=holo
230 j=m(k)
if(j-k) 190,190,240
240 ki=k-n
do 250 i=1,n
ki=ki+n
holo=a(ki)
ji=ki+j-k
a(ki)=-a(ji)
250 a(ji)=holo
goto 190
260 return
end
c
c
subroutine h1h2h1(a,b,c,d,nr,non,max)
c matrix * matrix * matrix (three triangle matrices)
implicit real*8(a-h,o-z)
dimension a(1),b(1),c(1),d(1),nr(1)
do 10 i=1,max
10 d(i)=0.d0
iofs=0
kl=0
do 20 i=1,non
ia=nr(i)
do 30 j=1,ia
ij=0
i10=iofs+j*(j-1)/2
jj=iofs+j
do 40 k=1,ia
i2=iofs+k*(k-1)/2
ij=ij+1
c(ij)=0.d0
do 50 l=1,k
i1=i10+l
if(l.gt.j) i1=l*(l-1)/2+jj
i2=i2+1
50 c(ij)=c(ij)+a(i1)*b(i2)
if(k.eq.ia) goto 40
i20=iofs+k
do 55 l=k+1,ia
i1=i10+l
if(l.gt.j) i1=l*(l-1)/2+jj
i2=i20+l*(l-1)/2
55 c(ij)=c(ij)+a(i1)*b(i2)
40 continue
do 60 m=1,j
mm=iofs+m
kl=kl+1
m20=m*(m-1)/2+iofs
do 70 n=1,ia
m2=m20+n
if(n.gt.m) m2=n*(n-1)/2+mm
70 d(kl)=d(kl)+c(n)*a(m2)
60 continue
30 continue
iofs=iofs+ia*(ia+1)/2
20 continue
return
end
c
c
function noteqv(i,j)
implicit real*8(a-h,o-z)
noteqv = i-j
return
end
c
c
real*8 function ggntde(i,j,eta,nfmx,nfcol)
implicit real*8(a-h,o-z)
dimension eta(nfmx,nfcol)
dimension aa(3),bb(3),a(10),b(10),e(3)
data c2/2.0d0/,c3/3.0d0/,const/15.7496099d0/
do 1 m=1,3
aa(m) = eta(i,m)
1 bb(m) = eta(j,m)
iff1=eta(i,5)+0.0001
iff2=eta(j,5)+0.0001
do 2 m=1,10
a(m) = eta(i,m+5)
2 b(m) = eta(j,m+5)
call divpt(aa,eta(i,4),bb,eta(j,4),e,eexp,efact)
call rhftce(a,aa,e,iff1)
call rhftce(b,bb,e,iff2)
fexp=c2*eexp
ggntde = efact*const*(a(1)*b(1)+(a(2)*b(2
1)+a(3)*b(3)+a(4)*b(4)+a(1)*(b(5)+b(6)+b(7))+b(1)*(a(5)+a(6)+a(7))+
2 (c3*(a(5)*b(5)+a(6)*b(6)+a(7)*b(7))+a(8)*b(8)+a(9)*b(9)+a(10)*b(1
30)+a(5)*(b(6)+b(7))+a(6)*(b(5)+b(7))+a(7)*(b(5)+b(6)))/fexp)/fexp)
4/(fexp*sqrt(fexp))
return
end
c
c
subroutine rhftce(cfs,a,e,iff)
implicit real*8(a-h,o-z)
dimension cfs(20) ,a(3),e(3)
data c2/2.0d0/,c3/3.0d0/
aex = e(1)-a(1)
aey = e(2)-a(2)
aez = e(3)-a(3)
goto (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20),iff
1 return
2 cfs(1)=aex*cfs(2)
return
3 cfs(1)=aey*cfs(3)
return
4 cfs(1)=aez*cfs(4)
return
5 cfs(1)=aex*aex*cfs(5)
cfs(2)=c2*aex*cfs(5)
return
6 cfs(1)=aey*aey*cfs(6)
cfs(3)=c2*aey*cfs(6)
return
7 cfs(1)=aez*aez*cfs(7)
cfs(4)=c2*aez*cfs(7)
return
8 cfs(1)=aex*aey*cfs(8)
cfs(2)=aey*cfs(8)
cfs(3)=aex*cfs(8)
return
9 cfs(1)=aex*aez*cfs(9)
cfs(2)=aez*cfs(9)
cfs(4)=aex*cfs(9)
return
10 cfs(1)=aey*aez*cfs(10)
cfs(3)=aez*cfs(10)
cfs(4)=aey*cfs(10)
return
11 cfs(1)=aex*aex*aex*cfs(11)
cfs(2)=c3*aex*aex*cfs(11)
cfs(5)=c3*aex*cfs(11)
return
12 cfs(1)=aey*aey*aey*cfs(12)
cfs(3)=c3*aey*aey*cfs(12)
cfs(6)=c3*aey*cfs(12)
return
13 cfs(1)=aez*aez*aez*cfs(13)
cfs(4)=c3*aez*aez*cfs(13)
cfs(7)=c3*aez*cfs(13)
return
14 cfs(1)=aex*aex*aey*cfs(14)
cfs(2)=c2*aex*aey*cfs(14)
cfs(3)=aex*aex*cfs(14)
cfs(5)=aey*cfs(14)
cfs(8)=c2*aex*cfs(14)
return
15 cfs(1)=aex*aex*aez*cfs(15)
cfs(2)=c2*aex*aez*cfs(15)
cfs(4)=aex*aex*cfs(15)
cfs(5)=aez*cfs(15)
cfs(9)=c2*aex*cfs(15)
return
16 cfs(1)=aey*aey*aex*cfs(16)
cfs(2)=aey*aey*cfs(16)
cfs(3)=c2*aey*aex*cfs(16)
cfs(6)=aex*cfs(16)
cfs(8)=c2*aey*cfs(16)
return
17 cfs(1)=aey*aey*aez*cfs(17)
cfs(3)=c2*aey*aez*cfs(17)
cfs(4)=aey*aey*cfs(17)
cfs(6)=aez*cfs(17)
cfs(10)=c2*aey*cfs(17)
return
18 cfs(1)=aez*aez*aex*cfs(18)
cfs(2)=aez*aez*cfs(18)
cfs(4)=c2*aez*aex*cfs(18)
cfs(7)=aex*cfs(18)
cfs(9)=c2*aez*cfs(18)
return
19 cfs(1)=aez*aez*aey*cfs(19)
cfs(3)=aez*aez*cfs(19)
cfs(4)=c2*aez*aey*cfs(19)
cfs(7)=aey*cfs(19)
cfs(10)=c2*aez*cfs(19)
return
20 cfs(1)=aex*aey*aez*cfs(20)
cfs(2)=aez*aey*cfs(20)
cfs(3)=aex*aez*cfs(20)
cfs(4)=aex*aey*cfs(20)
cfs(8)=aez*cfs(20)
cfs(9)=aey*cfs(20)
cfs(10)=aex*cfs(20)
return
end
c
c
subroutine dumdum(nuc,noc,vlist,c,dint,nro,nocmx,nt)
implicit real*8(a-h,o-z)
dimension vlist(nocmx,4), c(1), dint(nro,1)
call nuc(noc,vlist,c,dint,nro,nocmx,nt)
return
end
c
c
subroutine nonuc(noc,vlist,c,dint,nro,nocmx,nt)
c for properties without nuclear part
implicit real*8(a-h,o-z)
dimension dint(nro,1)
dint(1,1)=0.0d0
nt = 1
return
end
c
c
subroutine fldnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of electric field
implicit real*8(a-h,o-z)
dimension dint(nro,1),vlist(nocmx ,4),c(3) ,fnuc(3)
do20 i=1,3
20 fnuc(i)=0.0d0
do21 i=1,noc
dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
if(dst.eq.0.0d0) go to 21
do 23 j=1,3
fnuc(j)=fnuc(j)+vlist(i,4)*(vlist(i,j)-c(j))/(sqrt(dst))**3
23 continue
21 continue
do22 l=1,3
22 dint(1,l) = -fnuc(l)
nt = 3
return
end
c
c
subroutine potnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of potential
implicit real*8(a-h,o-z)
dimension dint(nro,1),vlist(nocmx ,4),c(3)
pnuc=0.0d0
do 12 i=1,noc
dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
if(dst.eq.0.0d0)go to 12
pnuc=pnuc+vlist(i,4)/sqrt(dst)
12 continue
dint(1,1)=pnuc
nt = 1
return
end
c
c
subroutine grdnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of field gradient
implicit real*8(a-h,o-z)
dimension dint(nro,*),vlist(nocmx ,4),c(3) ,gnuc(6)
do 30 i=1,6
30 gnuc(i)=0.0d0
do 31 i=1,noc
dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
if(dst.eq.0.0d0) go to 31
mu=0
do 32 j=1,3
do 33 k=1,j
mu=mu+1
gnuc(mu)=gnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
1/(sqrt(dst))**5
33 continue
32 continue
31 continue
dst=gnuc(1)+gnuc(3)+gnuc(6)
dint(1,1)=3.0d0*gnuc(1)-dst
dint(1,2)=3.0d0*gnuc(3)-dst
dint(1,3)=3.0d0*gnuc(6)-dst
dint(1,4)=gnuc(2)*3.0d0
dint(1,5)=gnuc(4)*3.0d0
dint(1,6)=gnuc(5)*3.0d0
do 40 i=1,6
40 dint(1,i) = - dint(1,i)
nt = 6
return
end
c
c
subroutine magnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of second moment
implicit real*8(a-h,o-z)
dimension dint(nro,*),vlist(nocmx ,4),c(3) ,qnuc(6)
do50 i=1,6
50 qnuc(i)=0.0d0
do 51 i=1,noc
mu=0
do52 j=1,3
do53 k=1,j
mu=mu+1
qnuc(mu)=qnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
53 continue
52 continue
51 continue
dint(1,1)=qnuc(1)
dint(1,2)=qnuc(3)
dint(1,3)=qnuc(6)
dint(1,4)=qnuc(2)
dint(1,5)=qnuc(4)
dint(1,6)=qnuc(5)
nt = 6
return
end
c
c
subroutine qudnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of quadropole moment
implicit real*8(a-h,o-z)
dimension dint(nro,*),vlist(nocmx ,*),c(3) ,qnuc(6)
data c3/3.0d0/
do50 i=1,6
50 qnuc(i)=0.0d0
do 51 i=1,noc
mu=0
do52 j=1,3
do53 k=1,j
mu=mu+1
qnuc(mu)=qnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
53 continue
52 continue
51 continue
dint(1,1)=qnuc(1)
dint(1,2)=qnuc(3)
dint(1,3)=qnuc(6)
dint(1,4)=qnuc(2)*1.5d0
dint(1,5)=qnuc(4)*1.5d0
dint(1,6)=qnuc(5)*1.5d0
dint(1,7)=dint(1,1)+dint(1,2)+dint(1,3)
dint(1,1)=(c3*dint(1,1)-dint(1,7))/2.d0
dint(1,2)=(c3*dint(1,2)-dint(1,7))/2.d0
dint(1,3) =(c3*dint(1,3)-dint(1,7))/2.d0
nt=7
return
end
c
c
subroutine dipnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of dipole moment
implicit real*8(a-h,o-z)
dimension dint(nro,*),vlist(nocmx ,4),c(3) ,dnuc(3)
do40 i=1,3
40 dnuc(i)=0.0d0
do41 i=1,noc
do41 j=1,3
41 dnuc(j)=dnuc(j)+vlist(i,4)*(vlist(i,j)-c(j))
do42 k=1,3
42 dint(1,k)=dnuc(k)
nt = 3
return
end
c
c
subroutine octnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of third moment
implicit real*8(a-h,o-z)
dimension dint(nro,*),vlist(nocmx,4),c(*),onuc(10)
do60 i=1,10
60 onuc(i)=0.0d0
do61 i=1,noc
mu=0
do62 j=1,3
do63 k=1,j
do64 m=1,k
mu=mu+1
onuc(mu)=onuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))*(
1vlist(i,m)-c(m))
64 continue
63 continue
62 continue
61 continue
dint(1,1)=onuc(1)
dint(1,2)=onuc(4)
dint(1,3)=onuc(10)
dint(1,4)=onuc(2)
dint(1,5)=onuc(5)
dint(1,6)=onuc(3)
dint(1,7)=onuc(7)
dint(1,8)=onuc(8)
dint(1,9)=onuc(9)
dint(1,10)=onuc(6)
nt = 10
return
end
c
c
subroutine sldnuc(noc,vlist,c,dint,nro,nocmx,nt)
c nuclear part of diamagnetic shielding
implicit real*8(a-h,o-z)
dimension dint(nro,*),vlist(nocmx,4),c(3),pnuc(6)
do 30 i=1,6
30 pnuc(i)=0.d0
do 31 i=1,noc
dst=(vlist(i,1)-c(1))**2+(vlist(i,2)-c(2))**2+(vlist(i,3)-c(3))**2
if(dst.lt.1.d-10) go to 31
mu=0
do 32 j=1,3
do 33 k=1,j
mu=mu+1
pnuc(mu)=pnuc(mu)+vlist(i,4)*(vlist(i,j)-c(j))*(vlist(i,k)-c(k))
1/(sqrt(dst))**3
33 continue
32 continue
31 continue
dint(1,1)=pnuc(1)
dint(1,2)=pnuc(3)
dint(1,3)=pnuc(6)
dint(1,4)=pnuc(2)
dint(1,5)=pnuc(4)
dint(1,6)=pnuc(5)
dint(1,7)=dint(1,1)+dint(1,2)+dint(1,3)
nt=7
return
end
c
c
subroutine propa(name,ij,kl,eta,va,c,nt,nbmx)
implicit real*8(a-h,o-z)
common /abfunc/ ra(3),rb(3),ga,gb,ia,ib
c this common-block is used in subroutine opap4 only
dimension eta(nbmx,15)
dimension va(1),c(1)
dimension a(3),b(3),e(3),d(3),aa(20),bb(20),dd(84),v(10),val(10)
dimension lin(84),min(84),nin(84)
data lin/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,1,0,1,0,1,4,0,0,3,3,1,0,1,0
1 ,2,2,0,2,1,1,5,0,0,3,3,2,2,0,0,4,4,1,0,0,1,1,3,1,2,2,1,6,0,0,3,
2 3,0,5,5,1,0,0,1,4,4,2,0,2,0,3,3,1,2,2,1,4,1,1,2/,min/0,0,1,0,0,
3 2,0,1,0,1,0,3,0,1,0,2,2,0,1,1,0,4,0,1,0,3,3,0,1,2,0,2,1,2,1,0,5,
4 0,2,0,3,0,3,2,1,0,4,4,1,0,1,1,3,2,1,2,0,6,0,3,0,3,1,0,0,1,5,5,2,0
5,0,2,4,4,2,1,3,1,3,2,1,4,1,2/,nin/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,0,
6 1,2,2,1,0,0,4,0,1,0,1,3,3,0,2,2,1,1,2,0,0,5,0,2,0,3,2,3,0,1,0,1,
7 4,4,3,1,1,1,2,2,0,0,6,0,3,3,0,1,5,5,1,0,0,2,4,4,0,2,1,2,2,3,1,3,
8 1,1,4,2/
do 10 i=1,3
a(i)=eta(ij,i)
b(i)=eta(kl,i)
ra(i)=a(i)
rb(i)=b(i)
10 continue
ga=eta(ij,4)
gb=eta(kl,4)
do 20 i=1,20
aa(i)=eta(ij,i+5)
bb(i)=eta(kl,i+5)
20 continue
iff1=eta(ij,5)+0.0001
iff2=eta(kl,5)+0.0001
ia=iff1
ib=iff2
call divpt(a,eta(ij,4),b,eta(kl,4),e,gama,efact)
call rhftce(aa,a,e,iff1)
call rhftce(bb,b,e,iff2)
call prod(aa,bb,dd,iff1,iff2)
do 30 i=1,nt
val(i)=0.d0
30 continue
do 25 i=1,3
d(i)=e(i)-c(i)
25 continue
c name represents an external function
if(iff1.gt.10.or.iff2.gt.10) goto 110
if(iff1.gt.4.or.iff2.gt.4) goto 120
if(iff1.gt.1.or.iff2.gt.1) goto 130
c s - s
call name(lin(1),min(1),nin(1),gama,v,d)
do 50 j=1,nt
val(j)=dd(1)*v(j)+val(j)
50 continue
do 70 i=1,nt
va(i)=efact*val(i)
70 continue
return
130 continue
if(iff1.ne.1.and.iff2.ne.1) goto 131
c s - p
do 61 i=1,4
if(abs(dd(i))-1.d-8) 61,61,41
41 call name(lin(i),min(i),nin(i),gama,v,d)
do 51 j=1,nt
val(j)=dd(i)*v(j)+val(j)
51 continue
61 continue
do 71 i=1,nt
va(i)=efact*val(i)
71 continue
return
131 continue
c p - p
do 62 i=1,10
if(abs(dd(i))-1.d-8) 62,62,42
42 call name(lin(i),min(i),nin(i),gama,v,d)
do 52 j=1,nt
val(j)=dd(i)*v(j)+val(j)
52 continue
62 continue
do 72 i=1,nt
va(i)=efact*val(i)
72 continue
return
120 continue
if(iff1.ne.1.and.iff2.ne.1) goto 121
c s - d
do 63 i=1,10
if(abs(dd(i))-1.d-8) 63,63,43
43 call name(lin(i),min(i),nin(i),gama,v,d)
do 53 j=1,nt
val(j)=dd(i)*v(j)+val(j)
53 continue
63 continue
do 73 i=1,nt
va(i)=efact*val(i)
73 continue
return
121 continue
if(iff1.gt.4.and.iff2.gt.4) goto 122
c p - d
do 64 i=1,20
if(abs(dd(i))-1.d-8) 64,64,44
44 call name(lin(i),min(i),nin(i),gama,v,d)
do 54 j=1,nt
val(j)=dd(i)*v(j)+val(j)
54 continue
64 continue
do 74 i=1,nt
va(i)=efact*val(i)
74 continue
return
122 continue
c d - d
do 65 i=1,35
if(abs(dd(i))-1.d-8) 65,65,45
45 call name(lin(i),min(i),nin(i),gama,v,d)
do 55 j=1,nt
val(j)=dd(i)*v(j)+val(j)
55 continue
65 continue
do 75 i=1,nt
va(i)=efact*val(i)
75 continue
return
110 continue
if(iff1.ne.1.and.iff2.ne.1) goto 111
c s - f
do 66 i=1,20
if(abs(dd(i))-1.d-8) 66,66,46
46 call name(lin(i),min(i),nin(i),gama,v,d)
do 56 j=1,nt
val(j)=dd(i)*v(j)+val(j)
56 continue
66 continue
do 76 i=1,nt
va(i)=efact*val(i)
76 continue
return
111 continue
if(iff1.gt.4.and.iff2.gt.4) goto 112
c p - f
do 67 i=1,35
if(abs(dd(i))-1.d-8) 67,67,47
47 call name(lin(i),min(i),nin(i),gama,v,d)
do 57 j=1,nt
val(j)=dd(i)*v(j)+val(j)
57 continue
67 continue
do 77 i=1,nt
va(i)=efact*val(i)
77 continue
return
112 continue
if(iff1.gt.10.and.iff2.gt.10) goto 113
c d - f
do 68 i=1,56
if(abs(dd(i))-1.d-8) 68,68,48
48 call name(lin(i),min(i),nin(i),gama,v,d)
do 58 j=1,nt
val(j)=dd(i)*v(j)+val(j)
58 continue
68 continue
do 78 i=1,nt
va(i)=efact*val(i)
78 continue
return
113 continue
c f - f
do 69 i=1,84
if(abs(dd(i))-1.d-8) 69,69,49
49 call name(lin(i),min(i),nin(i),gama,v,d)
do 59 j=1,nt
val(j)=dd(i)*v(j)+val(j)
59 continue
69 continue
do 79 i=1,nt
va(i)=efact*val(i)
79 continue
return
end
c
c
real*8 function olap(l,m,n,gama)
implicit real*8(a-h,o-z)
dimension dftr(7)
data pi /3.14159265359d0/
data dftr /1.d0,1.d0,3.d0,15.d0,105.d0,945.d0,10395.d0/
lh=l/2
mh=m/2
nh=n/2
if(2*lh-l) 50,1,50
1 if(2*mh-m) 50,2,50
2 if(2*nh-n) 50,3,50
3 hf=.5d0
olap=(sqrt(pi/gama))**3*(hf/gama)**(lh+mh+nh)*dftr(lh+1)*dftr(mh+
11)*dftr(nh+1)
return
50 olap=0.d0
return
end
c
c
subroutine prod(c,d,s,iff1,iff2)
implicit real*8(a-h,o-z)
dimension c(20),d(20),s(84)
do 10 i=1,84
10 s(i)=0.d0
if(iff1.gt.10.or.iff2.gt.10) goto 30
if(iff1.gt.4.or.iff2.gt.4) goto 20
s( 1)=c( 1)*d( 1)
c end of s - s
if(iff1.eq.1.and.iff2.eq.1) return
s( 2)=c( 1)*d( 2)+c( 2)*d( 1)
s( 3)=c( 1)*d( 3)+c( 3)*d( 1)
s( 4)=c( 1)*d( 4)+c( 4)*d( 1)
c end of s - p
if(iff1.eq.1.or.iff2.eq.1) return
s( 5)=c( 2)*d( 2)
s( 6)=c( 3)*d( 3)
s( 7)=c( 4)*d( 4)
s( 8)=c( 2)*d( 3)+c( 3)*d( 2)
s( 9)=c( 2)*d( 4)+c( 4)*d( 2)
s(10)=c( 3)*d( 4)+c( 4)*d( 3)
c end of p - p
return
20 continue
s( 1)=c( 1)*d( 1)
s( 2)=c( 1)*d( 2)+c( 2)*d( 1)
s( 3)=c( 1)*d( 3)+c( 3)*d( 1)
s( 4)=c( 1)*d( 4)+c( 4)*d( 1)
s( 5)=c( 1)*d( 5)+c( 5)*d( 1)+c( 2)*d( 2)
s( 6)=c( 1)*d( 6)+c( 6)*d( 1)+c( 3)*d( 3)
s( 7)=c( 1)*d( 7)+c( 7)*d( 1)+c( 4)*d( 4)
s( 8)=c( 1)*d( 8)+c( 8)*d( 1)+c( 2)*d( 3)+c( 3)*d( 2)
s( 9)=c( 1)*d( 9)+c( 9)*d( 1)+c( 2)*d( 4)+c( 4)*d( 2)
s(10)=c( 1)*d(10)+c(10)*d( 1)+c( 3)*d( 4)+c( 4)*d( 3)
c end of s - d
if(iff1.eq.1.or.iff2.eq.1) return
s(11)=c( 2)*d( 5)+c( 5)*d( 2)
s(12)=c( 3)*d( 6)+c( 6)*d( 3)
s(13)=c( 4)*d( 7)+c( 7)*d( 4)
s(14)=c( 2)*d( 8)+c( 8)*d( 2)+c( 3)*d( 5)+c( 5)*d( 3)
s(15)=c( 2)*d( 9)+c( 9)*d( 2)+c( 4)*d( 5)+c( 5)*d( 4)
s(16)=c( 2)*d( 6)+c( 6)*d( 2)+c( 3)*d( 8)+c( 8)*d( 3)
s(17)=c( 3)*d(10)+c(10)*d( 3)+c( 4)*d( 6)+c( 6)*d( 4)
s(18)=c( 2)*d( 7)+c( 7)*d( 2)+c( 4)*d( 9)+c( 9)*d( 4)
s(19)=c( 3)*d( 7)+c( 7)*d( 3)+c( 4)*d(10)+c(10)*d( 4)
s(20)=c( 2)*d(10)+c(10)*d( 2)+c( 3)*d( 9)
1 +c( 9)*d( 3)+c( 4)*d( 8)+c( 8)*d( 4)
c end of p - d
if(iff1.lt.5.or.iff2.lt.5) return
s(21)=c( 5)*d( 5)
s(22)=c( 6)*d( 6)
s(23)=c( 7)*d( 7)
s(24)=c( 5)*d( 8)+c( 8)*d( 5)
s(25)=c( 5)*d( 9)+c( 9)*d( 5)
s(26)=c( 6)*d( 8)+c( 8)*d( 6)
s(27)=c( 6)*d(10)+c(10)*d( 6)
s(28)=c( 7)*d( 9)+c( 9)*d( 7)
s(29)=c( 7)*d(10)+c(10)*d( 7)
s(30)=c( 5)*d( 6)+c( 6)*d( 5)+c( 8)*d( 8)
s(31)=c( 5)*d( 7)+c( 7)*d( 5)+c( 9)*d( 9)
s(32)=c( 6)*d( 7)+c( 7)*d( 6)+c(10)*d(10)
s(33)=c( 5)*d(10)+c(10)*d( 5)+c( 8)*d( 9)+c( 9)*d( 8)
s(34)=c( 6)*d( 9)+c( 9)*d( 6)+c( 8)*d(10)+c(10)*d( 8)
s(35)=c( 7)*d( 8)+c( 8)*d( 7)+c( 9)*d(10)+c(10)*d( 9)
c end of d - d
return
30 continue
s( 1)=c( 1)*d( 1)
s( 2)=c( 1)*d( 2)+c( 2)*d( 1)
s( 3)=c( 1)*d( 3)+c( 3)*d( 1)
s( 4)=c( 1)*d( 4)+c( 4)*d( 1)
s( 5)=c( 1)*d( 5)+c( 5)*d( 1)+c( 2)*d( 2)
s( 6)=c( 1)*d( 6)+c( 6)*d( 1)+c( 3)*d( 3)
s( 7)=c( 1)*d( 7)+c( 7)*d( 1)+c( 4)*d( 4)
s( 8)=c( 1)*d( 8)+c( 8)*d( 1)+c( 2)*d( 3)+c( 3)*d( 2)
s( 9)=c( 1)*d( 9)+c( 9)*d( 1)+c( 2)*d( 4)+c( 4)*d( 2)
s(10)=c( 1)*d(10)+c(10)*d( 1)+c( 3)*d( 4)+c( 4)*d( 3)
s(11)=c( 2)*d( 5)+c( 5)*d( 2)+c(1)*d(11)+c(11)*d(1)
s(12)=c( 3)*d( 6)+c( 6)*d( 3)+c(1)*d(12)+c(12)*d(1)
s(13)=c( 4)*d( 7)+c( 7)*d( 4)+c(1)*d(13)+c(13)*d(1)
s(14)=c( 2)*d( 8)+c( 8)*d( 2)+c( 3)*d( 5)+c( 5)*d( 3)+c(1)*d(14)+
1 c(14)*d(1)
s(15)=c( 2)*d( 9)+c( 9)*d( 2)+c( 4)*d( 5)+c( 5)*d( 4)+c(1)*d(15)+
1 c(15)*d(1)
s(16)=c( 2)*d( 6)+c( 6)*d( 2)+c( 3)*d( 8)+c( 8)*d( 3)+c(1)*d(16)+
1 c(16)*d(1)
s(17)=c( 3)*d(10)+c(10)*d( 3)+c( 4)*d( 6)+c( 6)*d( 4)+c(1)*d(17)+
1 c(17)*d(1)
s(18)=c( 2)*d( 7)+c( 7)*d( 2)+c( 4)*d( 9)+c( 9)*d( 4)+c(1)*d(18)+
1 c(18)*d(1)
s(19)=c( 3)*d( 7)+c( 7)*d( 3)+c( 4)*d(10)+c(10)*d( 4)+c(1)*d(19)+
1 c(19)*d(1)
s(20)=c( 2)*d(10)+c(10)*d( 2)+c( 3)*d( 9)+c(9)*d(3)+c(4)*d(8)+
1 c(8)*d(4)+c(1)*d(20)+c(20)*d(1)
c end of s - f
if(iff1.eq.1.or.iff2.eq.1) return
s(21)=c( 5)*d( 5)+c(2)*d(11)+c(11)*d(2)
s(22)=c( 6)*d( 6)+c(3)*d(12)+c(12)*d(3)
s(23)=c( 7)*d( 7)+c(4)*d(13)+c(13)*d(4)
s(24)=c( 5)*d( 8)+c( 8)*d( 5)+c(3)*d(11)+c(11)*d(3)+c(2)*d(14)+
1 c(14)*d(2)
s(25)=c( 5)*d( 9)+c( 9)*d( 5)+c(2)*d(15)+c(15)*d(2)+c(4)*d(11)+
1 c(11)*d(4)
s(26)=c( 6)*d( 8)+c( 8)*d( 6)+c(2)*d(12)+c(12)*d(2)+c(3)*d(16)+
1 c(16)*d(3)
s(27)=c( 6)*d(10)+c(10)*d( 6)+c(3)*d(17)+c(17)*d(3)+c(4)*d(12)+
1 c(12)*d(4)
s(28)=c( 7)*d( 9)+c( 9)*d( 7)+c(2)*d(13)+c(13)*d(2)+c(4)*d(18)+
1 c(18)*d(4)
s(29)=c( 7)*d(10)+c(10)*d( 7)+c(3)*d(13)+c(13)*d(3)+c(4)*d(19)+
1 c(19)*d(4)
s(30)=c( 5)*d( 6)+c( 6)*d( 5)+c( 8)*d( 8)+c(2)*d(16)+c(16)*d(2)+
1 c(3)*d(14)+c(14)*d(3)
s(31)=c( 5)*d( 7)+c( 7)*d( 5)+c( 9)*d( 9)+c(2)*d(18)+c(18)*d(2)+
1 c(4)*d(15)+c(15)*d(4)
s(32)=c( 6)*d( 7)+c( 7)*d( 6)+c(10)*d(10)+c(3)*d(19)+c(19)*d(3)+
1 c(4)*d(17)+c(17)*d(4)
s(33)=c( 5)*d(10)+c(10)*d( 5)+c( 8)*d( 9)+c( 9)*d( 8)+c(3)*d(15)+
1 c(15)*d(3)+c(4)*d(14)+c(14)*d(4)+c(2)*d(20)+c(20)*d(2)
s(34)=c( 6)*d( 9)+c( 9)*d( 6)+c( 8)*d(10)+c(10)*d( 8)+c(2)*d(17)+
1 d(2)*c(17)+c(3)*d(20)+c(20)*d(3)+c(4)*d(16)+c(16)*d(4)
s(35)=c( 7)*d( 8)+c( 8)*d( 7)+c( 9)*d(10)+c(10)*d( 9)+c(2)*d(19)+
1 c(19)*d(2)+c(3)*d(18)+c(18)*d(3)+c(4)*d(20)+c(20)*d(4)
c end of p - f
if(iff1.eq.2.or.iff2.eq.2) return
s(36)=c(5)*d(11)+c(11)*d(5)
s(37)=c(6)*d(12)+c(12)*d(6)
s(38)=c(7)*d(13)+c(13)*d(7)
s(39)=c(6)*d(11)+c(11)*d(6)+c(5)*d(16)+c(16)*d(5)+c(8)*d(14)+
1 c(14)*d(8)
s(40)=c(7)*d(11)+c(11)*d(7)+c(5)*d(18)+c(18)*d(5)+c(9)*d(15)+
1 c(15)*d(9)
s(41)=c(5)*d(12)+c(12)*d(5)+c(6)*d(14)+c(14)*d(6)+c(8)*d(16)+
1 c(16)*d(8)
s(42)=c(5)*d(13)+c(13)*d(5)+c(7)*d(15)+c(15)*d(7)+c(9)*d(18)+
1 c(18)*d(9)
s(43)=c(7)*d(12)+c(12)*d(7)+c(6)*d(19)+c(19)*d(6)+c(10)*d(17)+
1 c(17)*d(10)
s(44)=c(6)*d(13)+c(13)*d(6)+c(7)*d(17)+c(17)*d(7)+c(10)*d(19)+
1 c(19)*d(10)
s(45)=c(8)*d(11)+c(11)*d(8)+c(5)*d(14)+c(14)*d(5)
s(46)=c(9)*d(11)+c(11)*d(9)+c(5)*d(15)+c(15)*d(5)
s(47)=c(8)*d(12)+c(12)*d(8)+c(6)*d(16)+c(16)*d(6)
s(48)=c(10)*d(12)+c(12)*d(10)+c(6)*d(17)+c(17)*d(6)
s(49)=c(10)*d(13)+c(13)*d(10)+c(7)*d(19)+c(19)*d(7)
s(50)=c(9)*d(13)+c(13)*d(9)+c(7)*d(18)+c(18)*d(7)
s(51)=c(8)*d(13)+c(13)*d(8)+c(7)*d(20)+c(20)*d(7)+c(9)*d(19)+
1 c(19)*d(9)+c(10)*d(18)+c(18)*d(10)
s(52)=c(10)*d(11)+c(11)*d(10)+c(5)*d(20)+c(20)*d(5)+c(9)*d(14)+
1 c(14)*d(9)+c(8)*d(15)+c(15)*d(8)
s(53)=c(9)*d(12)+c(12)*d(9)+c(6)*d(20)+c(20)*d(6)+c(10)*d(16)+
1 c(16)*d(10)+c(8)*d(17)+c(17)*d(8)
s(54)=c(5)*d(17)+c(17)*d(5)+c(6)*d(15)+c(15)*d(6)+c(14)*d(10)+
1 d(14)*c(10)+c(9)*d(16)+c(16)*d(9)+c(8)*d(20)+c(20)*d(8)
s(55)=c(5)*d(19)+c(19)*d(5)+c(7)*d(14)+c(14)*d(7)+c(10)*d(15)+
1 c(15)*d(10)+c(9)*d(20)+c(20)*d(9)+c(8)*d(18)+c(18)*d(8)
s(56)=c(6)*d(18)+c(18)*d(6)+c(7)*d(16)+c(16)*d(7)+c(10)*d(20)+
1 c(20)*d(10)+c(9)*d(17)+c(17)*d(9)+c(8)*d(19)+c(19)*d(8)
c end of d - f
if(iff1.eq.3.or.iff2.eq.3) return
s(57)=c(11)*d(11)
s(58)=c(12)*d(12)
s(59)=c(13)*d(13)
s(60)=c(11)*d(12)+c(12)*d(11)+c(14)*d(16)+c(16)*d(14)
s(61)=c(11)*d(13)+c(13)*d(11)+c(15)*d(18)+c(18)*d(15)
s(62)=c(12)*d(13)+c(13)*d(12)+c(17)*d(19)+c(19)*d(17)
s(63)=c(11)*d(14)+c(14)*d(11)
s(64)=c(11)*d(15)+c(15)*d(11)
s(65)=c(13)*d(18)+c(18)*d(13)
s(66)=c(13)*d(19)+c(19)*d(13)
s(67)=c(12)*d(17)+d(12)*c(17)
s(68)=c(12)*d(16)+c(16)*d(12)
s(69)=c(11)*d(16)+c(16)*d(11)+c(14)*d(14)
s(70)=c(11)*d(18)+c(18)*d(11)+c(15)*d(15)
s(71)=c(13)*d(15)+c(15)*d(13)+c(18)*d(18)
s(72)=c(13)*d(17)+c(17)*d(13)+c(19)*d(19)
s(73)=c(12)*d(14)+c(14)*d(12)+c(16)*d(16)
s(74)=c(12)*d(19)+c(19)*d(12)+c(17)*d(17)
s(75)=c(11)*d(17)+c(17)*d(11)+c(14)*d(20)+c(20)*d(14)+
1 c(15)*d(16)+c(16)*d(15)
s(76)=c(11)*d(19)+c(19)*d(11)+c(20)*d(15)+c(15)*d(20)+
1 c(14)*d(18)+c(18)*d(14)
s(77)=c(12)*d(18)+c(18)*d(12)+c(16)*d(19)+c(19)*d(16)+
1 c(17)*d(20)+c(20)*d(17)
s(78)=c(13)*d(14)+c(14)*d(13)+c(15)*d(19)+c(19)*d(15)+
1 c(18)*d(20)+c(20)*d(18)
s(79)=c(12)*d(15)+c(15)*d(12)+c(14)*d(17)+c(17)*d(14)+
1 c(16)*d(20)+c(20)*d(16)
s(80)=c(13)*d(16)+c(16)*d(13)+c(17)*d(18)+c(18)*d(17)+
1 c(19)*d(20)+c(20)*d(19)
s(81)=c(11)*d(20)+c(20)*d(11)+c(14)*d(15)+c(15)*d(14)
s(82)=c(12)*d(20)+c(20)*d(12)+c(16)*d(17)+c(17)*d(16)
s(83)=c(13)*d(20)+c(20)*d(13)+c(18)*d(19)+c(19)*d(18)
s(84)=c(14)*d(19)+c(19)*d(14)+c(15)*d(17)+c(17)*d(15)+
1 c(16)*d(18)+c(18)*d(16)+c(20)*d(20)
c end of f - f
return
end
c
c
subroutine fmc(mvar,xvar,expmx,fmch)
implicit real*8(a-h,o-z)
m=mvar
x=xvar
expmx=exp(-x)
if(x -20.d0) 10,10,20
c x less than or equal to 20.0.
10 a=m
a=a+0.5d0
term=1.d0/a
ptlsum=term
do 11 i=2,60
a=a+1.d0
term=term*x/a
ptlsum=ptlsum+term
if (term/ptlsum- 1.d-8 ) 12, 11, 11
11 continue
print 999,m,x
12 fmch=.5d0*ptlsum*expmx
return
c x greater than 20.0.
20 a=m
xd=1.d0/x
b=a+0.5d0
a=a-0.5d0
approx=.886226925428d0*(sqrt(xd)*xd**m)
if (m) 21, 23, 21
21 do 22 i=1,m
b=b-1.d0
22 approx=approx*b
23 fimult=.5d0*expmx*xd
fiprop=fimult/approx
term=1.d0
ptlsum=term
notrms=x
notrms=notrms+m
do 24 i=2,notrms
term=term*a*xd
ptlsum=ptlsum+term
if (abs(term*fiprop/ptlsum) - 1.d-8 ) 25,25,24
24 a=a-1.d0
print 999,m,x
25 fmch=approx-fimult*ptlsum
return
999 format (24h no convergence for fmch, i6, d16.9)
end
c
c
subroutine opaa0(l,m,n,ga,v,d)
c electronic part of potential
implicit real*8(a-h,o-z)
dimension v(1),d(3),fnu(7),fn(7),fd(7),dp(4)
data pi /3.14159265359d0/, fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,
1 720.d0/
data fd/1.d0,1.d0,0.5d0,0.1666666666667d0,0.416666666667d-1,
1 .833333333333d-2,.1388888888889d-2/
itot=l+m+n
itoth=itot/2
pre=(2.d0*pi/ga)*fn(l+1)*fn(m+1)*fn(n+1)
if(2*itoth-itot) 1,2,1
1 pre=-pre
2 del=.25d0/ga
x=ga*(d(1)**2+d(2)**2+d(3)**2)
xx=2.d0*x
call fmc(6,x,expmx,fmch)
fnu(7)=fmch
fnu(6)=(expmx+xx*fnu(7))/11.d0
fnu(5)=(expmx+xx*fnu(6))/9.d0
fnu(4)=(expmx+xx*fnu(5))/7.d0
fnu(3)=(expmx+xx*fnu(4))/5.d0
fnu(2)=(expmx+xx*fnu(3))/3.d0
fnu(1)=expmx+xx*fnu(2)
dp(1)=1.d0
dp(2)=del
dp(3)=del**2
v(1)=pre*aainer(0,0,0,l,m,n,d,dp,fnu,fn,fd)*(-1.0d0)
return
end
c
c
subroutine opaa1(l,m,n,ga,v,d)
c electronic part of electric field
implicit real*8(a-h,o-z)
dimension v(3),d(3),fnu(9),fd(9),fn(9),dp(4)
data pi /3.14159265359d0/, fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,
1 720.d0,5040.d0,40320.d0/
data fd /1.d0,1.d0,.5d0,.1666666666667d0,.416666666667d-1,
1 .833333333333d-2,.1388888888889d-2,.1984126984d-3,0.248015873d-4/
itot=l+m+n
itoth=itot/2
pre=4.d0*pi*fn(l+1)*fn(m+1)*fn(n+1)
if(2*itoth-itot) 1,2,1
1 pre=-pre
2 del=.25d0/ga
x=ga*(d(1)**2+d(2)**2+d(3)**2)
xx=2.d0*x
call fmc(8,x,expmx,fmch)
fnu(9)=fmch
fnu(8)=(expmx+xx*fnu(9))/15.d0
fnu(7)=(expmx+xx*fnu(8))/13.d0
fnu(6)=(expmx+xx*fnu(7))/11.d0
fnu(5)=(expmx+xx*fnu(6))/9.d0
fnu(4)=(expmx+xx*fnu(5))/7.d0
fnu(3)=(expmx+xx*fnu(4))/5.d0
fnu(2)=(expmx+xx*fnu(3))/3.d0
fnu(1)=expmx+xx*fnu(2)
dp(1)=1.d0
dp(2)=del
dp(3)=del**2
dp(4)=del*dp(3)
v(1)=pre*aainer(1,0,0,l,m,n,d,dp,fnu,fn,fd)
v(2)=pre*aainer(0,1,0,l,m,n,d,dp,fnu,fn,fd)
v(3)=pre*aainer(0,0,1,l,m,n,d,dp,fnu,fn,fd)
return
end
c
c
subroutine opaa2(l,m,n,ga,v,d)
c electronic part of field gradient
implicit real*8(a-h,o-z)
dimension v(6),d(3),fnu(9),fd(9),fn(9),dp(5)
data pi /3.14159265359d0/, fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,
1 720.d0,5040.d0,40320.d0/
data fd /1.d0,1.d0,.5d0,.1666666666667d0,.416666666667d-1,
1 .833333333333d-2,.1388888888889d-2,.1984126984d-3,0.248015873d-4/
itot=l+m+n
itoth=itot/2
pre=(8.d0*pi*ga)*fn(l+1)*fn(m+1)*fn(n+1)
if(2*itoth-itot)1,2,1
1 pre=-pre
2 del=.25d0/ga
x=ga*(d(1)**2+d(2)**2+d(3)**2)
xx=2.d0*x
call fmc(8,x,expmx,fmch)
fnu(9)=fmch
fnu(8)=(expmx+xx*fnu(9))/15.d0
fnu(7)=(expmx+xx*fnu(8))/13.d0
fnu(6)=(expmx+xx*fnu(7))/11.d0
fnu(5)=(expmx+xx*fnu(6))/9.d0
fnu(4)=(expmx+xx*fnu(5))/7.d0
fnu(3)=(expmx+xx*fnu(4))/5.d0
fnu(2)=(expmx+xx*fnu(3))/3.d0
fnu(1)=expmx+xx*fnu(2)
dp(1)=1.d0
dp(2)=del
dp(3)=del**2
dp(4)=del*dp(3)
dp(5)=del*dp(4)
v(1)=pre*aainer(2,0,0,l,m,n,d,dp,fnu,fn,fd)
v(2)=pre*aainer(0,2,0,l,m,n,d,dp,fnu,fn,fd)
v(3)=pre*aainer(0,0,2,l,m,n,d,dp,fnu,fn,fd)
v(4)=pre*aainer(1,1,0,l,m,n,d,dp,fnu,fn,fd)
v(5)=pre*aainer(1,0,1,l,m,n,d,dp,fnu,fn,fd)
v(6)=pre*aainer(0,1,1,l,m,n,d,dp,fnu,fn,fd)
call opac3(l,m,n,ga,dp(1),d)
dp(1)=1.333333333333d0*pi*dp(1)
v(1)=v(1)+dp(1)
v(2)=v(2)+dp(1)
v(3)=v(3)+dp(1)
return
end
c
c
subroutine opaa3(l,m,n,ga,v,d)
c electronic part of diamagnetic shielding
implicit real*8(a-h,o-z)
dimension vx(3),d(3),v(7)
call opaa1(l,m,n,ga,vx,d)
v(1)=vx(1)*d(1)
v(2)=vx(2)*d(2)
v(3)=vx(3)*d(3)
v(4)=vx(2)*d(1)
v(5)=vx(3)*d(1)
v(6)=vx(3)*d(2)
call opaa1(l+1,m,n,ga,vx,d)
v(1)=v(1)+vx(1)
v(4)=v(4)+vx(2)
v(5)=v(5)+vx(3)
call opaa1(l,m+1,n,ga,vx,d)
v(2)=v(2)+vx(2)
v(6)=v(6)+vx(3)
call opaa1(l,m,n+1,ga,vx,d)
v(3)=v(3)+vx(3)
v(7)=v(1)+v(2)+v(3)
do 300 i=1,7
v(i)=-v(i)
300 continue
return
end
c
c
subroutine opab1(l,m,n,ga,v,d)
c electronic part of dipole moment
implicit real*8(a-h,o-z)
dimension v(3),d(3),g(4)
g(1)=olap(l,m,n,ga)
g(2)=olap(l+1,m,n,ga)
g(3)=olap(l,m+1,n,ga)
g(4)=olap(l,m,n+1,ga)
v(1)=g(2)+d(1)*g(1)
v(2)=g(3)+d(2)*g(1)
v(3)=g(4)+d(3)*g(1)
do 300 i=1,3
300 v(i)=-v(i)
return
end
c
c
subroutine opab2(l,m,n,ga,v,d)
c electronic part of quadrupole moment
implicit real*8(a-h,o-z)
dimension v(7),d(3),g(10)
g(1)=olap(l,m,n,ga)
g(2)=olap(l+1,m,n,ga)
g(3)=olap(l,m+1,n,ga)
g(4)=olap(l,m,n+1,ga)
g(5)=olap(l+2,m,n,ga)
g(6)=olap(l,m+2,n,ga)
g(7)=olap(l,m,n+2,ga)
g(8)=olap(l+1,m+1,n,ga)
g(9)=olap(l+1,m,n+1,ga)
g(10)=olap(l,m+1,n+1,ga)
v(1)=g(5)+2.d0*d(1)*g(2)+d(1)**2*g(1)
v(2)=g(6)+2.d0*d(2)*g(3)+d(2)**2*g(1)
v(3)=g(7)+2.d0*d(3)*g(4)+d(3)**2*g(1)
v(4)=(g(8)+d(1)*g(3)+d(2)*g(2)+d(1)*d(2)*g(1))/2.d0
v(5)=(g(9)+d(1)*g(4)+ d(3)*g(2)+d(1)*d(3)*g(1))/2.d0
v(6)=(g(10)+d(2)*g(4)+d(3)*g(3)+d(2)*d(3)*g(1))/2.d0
v(7)=v(1)+v(2)+v(3)
v(1)=(3.d0*v(1)-v(7))/2.d0
v(2)=(3.d0*v(2)-v(7))/2.d0
v(3)=(3.d0*v(3)-v(7))/2.d0
v(4)=3.d0*v(4)
v(5)=3.d0*v(5)
v(6)=3.d0*v(6)
do 300 i=1,7
300 v(i)=-v(i)
return
end
c
c
subroutine opab3(l,m,n,ga,v,dd)
c electronic part of third moment
implicit real*8(a-h,o-z)
dimension v(10),dd(3),g(20),d(20)
g(1)=olap( l,m,n,ga)
g(2)=olap(l+1,m,n,ga)
g(3)=olap(l,m+1,n,ga)
g(4)=olap(l,m,n+1,ga)
g(5)=olap(l+2,m,n,ga)
g(6)=olap(l,m+2,n,ga)
g(7)=olap(l,m,n+2,ga)
g(8)=olap(l+1,m+1,n,ga)
g(9)=olap(l+1,m,n+1,ga)
g(10)=olap(l,m+1,n+1,ga)
g(11)=olap(l+3,m,n,ga)
g(12)=olap(l,m+3,n,ga)
g(13)=olap(l,m,n+3,ga)
g(14)=olap(l+2,m+1,n,ga)
g(15)=olap(l+2,m,n+1,ga)
g(16)=olap(l+1,m+2,n,ga)
g(17)=olap(l,m+2,n+1,ga)
g(18)=olap(l+1,m,n+2,ga)
g(19)=olap(l,m+1,n+2,ga)
g(20)=olap(l+1,m+1,n+1,ga)
d(2)=dd(1)
d(3)=dd(2)
d(4)=dd(3)
d(5)=dd(1)**2
d(6)=dd(2)**2
d(7)=dd(3)**2
d(8)=d(2)*d(3)
d(9)=d(2)*d(4)
d(10)=d(3)*d(4)
d(11)=d(5)*d(2)
d(12)=d(6)*d(3)
d(13)=d(7)*d(4)
d(14)=d(5)*d(3)
d(15)=d(5)*d(4)
d(16)=d(6)*d(2)
d(17)=d(6)*d(4)
d(18)=d(7)*d(2)
d(19)=d(7)*d(3)
d(20)=d(8)*d(4)
v(1)=g(11)+3.d0*d(2)*g(5)+3.d0*d(5)*g(2)+d(11)*g(1)
v(2)=g(12)+3.d0*d(3)*g(6)+3.d0*d(6)*g(3)+d(12)*g(1)
v(3)=g(13)+3.d0*d(4)*g(7)+3.d0*d(7)*g(4)+d(13)*g(1)
v(4)=g(14)+2.d0*d(2)*g(8)+d(3)*g(5)+d(5)*g(3)+2.d0*d(8)*g(2)
1 +d(14)*g(1)
v(5)=g(15)+2.d0*d(2)*g(9)+d(4)*g(5)+d(5)*g(4)+2.d0*d(9)*g(2)
1 +d(15)*g(1)
v(6)=g(16)+2.d0*d(3)*g(8)+d(2)*g(6)+d(6)*g(2)+2.d0*d(8)*g(3)
1 +d(16)*g(1)
v(7)=g(17)+2.d0*d(3)*g(10)+d(4)*g(6)+d(6)*g(4)+2.d0*d(10)*g(3)
1 +d(17)*g(1)
v(8)=g(18)+2.d0*d(4)*g(9)+d(2)*g(7)+d(7)*g(2)+2.d0*d(9)*g(4)
1 +d(18)*g(1)
v(9)=g(19)+2.d0*d(4)*g(10)+d(3)*g(7)+d(7)*g(3)+2.d0*d(10)*g(4)
1 +d(19)*g(1)
v(10)=g(20)+d(2)*g(10)+d(3)*g(9)+d(4)*g(8)+d(10)*g(2)+d(9)*g(3)+
1d(8)*g(4)+d(20)*g(1)
do 300 i=1,10
300 v(i)=-v(i)
return
end
c
c
subroutine opab4(l,m,n,ga,v,d)
c electronic part of second moment
implicit real*8(a-h,o-z)
dimension v(6),d(3),g(10)
g(1)=olap(l,m,n,ga)
g(2)=olap(l+1,m,n,ga)
g(3)=olap(l,m+1,n,ga)
g(4)=olap(l,m,n+1,ga)
g(5)=olap(l+2,m,n,ga)
g(6)=olap(l,m+2,n,ga)
g(7)=olap(l,m,n+2,ga)
g(8)=olap(l+1,m+1,n,ga)
g(9)=olap(l+1,m,n+1,ga)
g(10)=olap(l,m+1,n+1,ga)
v(1)=g(5)+2.d0*d(1)*g(2)+d(1)**2*g(1)
v(2)=g(6)+2.d0*d(2)*g(3)+d(2)**2*g(1)
v(3)=g(7)+2.d0*d(3)*g(4)+d(3)**2*g(1)
v(4)=g(8)+d(1)*g(3)+d(2)*g(2)+d(1)*d(2)*g(1)
v(5)=g(9)+d(1)*g(4)+d(3)*g(2)+d(1)*d(3)*g(1)
v(6)=g(10)+d(2)*g(4)+d(3)*g(3)+d(2)*d(3)*g(1)
do 300 i=1,6
300 v(i)=-v(i)
return
end
c
c
subroutine opac1(l,m,n,ga,v,d)
c electronic part of planar density
implicit real*8(a-h,o-z)
dimension v(3),d(3),dd(3),dftr(6)
data pi /3.14159265359d0/, dftr /1.d0,1.d0,3.d0,5.d0,15.d0,105.d0/
do 1 i=1,3
dd(i)=-d(i)
v(i)=0.d0
1 continue
g=.5d0/ga
lh=l/2
mh=m/2
nh=n/2
pre=2.d0*g*pi
if(2*nh-n) 9,3,9
3 if(2*mh-m) 7,4,7
4 v(1)=1.d0
if(l) 39,41,39
39 v(1)=dd(1)**l
41 v(1)=pre*v(1)* exp (-ga*dd(1)**2)*g**(mh+nh)*dftr(mh+1)*dftr(nh
1+1)
if(2*lh-l) 50,5,50
5 v(2)=1.d0
if(m) 49,51,49
49 v(2)=dd(2)**m
51 v(2)=pre*v(2)* exp (-ga*dd(2)**2)*g**(lh+nh)*dftr(lh+1)*dftr(nh
1+1)
6 v(3)=1.d0
if(n) 59,61,59
59 v(3)=dd(3)**n
61 v(3)=pre*v(3)* exp (-ga*dd(3)**2)*g**(lh+mh)*dftr(lh+1)*dftr(mh
1+1)
return
7 if(2*lh-l) 50,8,50
8 v(2)=1.d0
if(m) 79,81,79
79 v(2)=dd(2)**m
81 v(2)=pre*v(2)* exp (-ga*dd(2)**2)*g**(lh+nh)*dftr(lh+1)*dftr(nh
1+1)
return
9 if(2*mh-m) 50,10,50
10 if(2*lh-l) 50,6,50
50 return
end
c
c
subroutine opac2(l,m,n,ga,v,d)
c electronic part of linear density
implicit real*8(a-h,o-z)
dimension v(3),d(3),dd(3),ds(3),dftr(5)
data sqrpi /1.772453850906d0/, dftr /1.d0,1.d0,3.d0,15.d0,105.d0/
do 1 i=1,3
dd(i)=-d(i)
ds(i)=d(i)**2
v(i)=0.d0
1 continue
g=.5d0/ga
pre=sqrpi*sqrt(2.d0*g)
lh=l/2
mh=m/2
nh=n/2
if(2*nh-n) 4,3,4
3 v(1)=1.d0
if(l.le.0.and.dd(1).eq.0.d0) go to 29
28 v(1)=dd(1)**l
29 if(m.le.0.and.dd(2).eq.0.d0) go to 31
30 v(1)=v(1)*dd(2)**m
31 v(1)=pre*g**nh*v(1)*exp(-ga*(ds(1)+ds(2)))*dftr(nh+1)
4 if(2*mh-m) 6,5,6
5 v(2)=1.d0
if(l.le.0.and.dd(1).eq.0.d0) go to 49
48 v(2)=dd(1)**l
49 if(n.le.0.and.dd(1).eq.0.d0) go to 51
50 v(2)=v(2)*dd(1)**n
51 v(2)=pre*g**mh*v(2)*exp(-ga*(ds(1)+ds(3)))*dftr(mh+1)
6 if(2*lh-l) 8,7,8
7 v(3)=1.d0
if(m.le.0.and.dd(2).eq.0.d0) go to 69
68 v(3)=dd(2)**m
69 if(n.le.0.and.dd(3).eq.0.d0) go to 71
70 v(3)=v(3)*dd(3)**n
71 v(3)=pre*g**lh*v(3)*exp(-ga*(ds(2)+ds(3)))*dftr(lh+1 )
8 return
end
c
c
subroutine opac3(l,m,n,ga,v,d)
c electronic part of charge density
implicit real*8(a-h,o-z)
dimension v(1),d(1),dd(3)
ddsq=0.d0
do 1 i=1,3
dd(i)=-d(i)
ddsq=ddsq+d(i)**2
1 continue
v(1) =1.d0
if(l) 2,3,2
2 v(1)=dd(1)**l
3 if(m) 4,5,4
4 v(1)=v(1)*dd(2)**m
5 if(n) 6,7,6
6 v(1)=v(1)*dd(3)**n
7 v(1)=v(1)*exp(-ga*ddsq)
return
end
c
c
subroutine opad1(l,m,n,ga,v,d)
c electronic part of overlap
implicit real*8(a-h,o-z)
dimension v(1),d(1)
v(1)=olap(l,m,n,ga)
return
end
c
c
subroutine opap4(l,m,n,gc,v,rc)
c
c this subroutine calculates the expectation-values of
c -p**2/2 (kinetic energy) and p**4/(8*cl**2) (relativistic
c correction of the kinetic energy)
c
c the program works correct only for s and p and the irreducible
c linear combinations of the primitive d,f,... functions
c
c written by s. brode in april,1984
c
c
implicit real*8 (a-h,o-z)
logical lodd,modd,nodd,leven,meven,neven
dimension v(3),rc(3),rca(3),rcb(3),srcab(3),
&lmni(35),o(25),dftr(7)
data a0/0.0d0/,a1/1.0d0/,a2/2.0d0/,a4/4.0d0/,a8/8.0d0/,
&lmni/0,3*1,6*2,10*3,15*4/,cl/137.03604d0/,pi/3.1415926559d0/,
&dftr/1.0d0,1.0d0,3.0d0,15.0d0,105.0d0,945.0d0,103935.0d0/
common /abfunc/ ra(3),rb(3),ga,gb,ia,ib
c
c define the overlap-integral
c
ola(ll,mm,nn)=gch**(ll/2+mm/2+nn/2)*dftr(ll/2+1)*dftr(mm/2+1)*
&dftr(nn/2+1)
c
c some previous calculations
c
do 10 i=1,25
10 o(i)=a0
gch=a1/(a2*gc)
ropigc=sqrt(pi/gc)*pi/gc
lodd=mod(l,2).eq.1
leven=.not.lodd
modd=mod(m,2).eq.1
meven=.not.modd
nodd=mod(n,2).eq.1
neven=.not.nodd
ga2=ga*a2
gb2=gb*a2
gab4=ga2*gb2
prcaa=a0
prcbb=a0
do 20 i=1,3
rca(i)=rc(i)-ra(i)
rcb(i)=rc(i)-rb(i)
srcab(i)=rca(i)+rcb(i)
prcaa=prcaa+rca(i)*rca(i)
prcbb=prcbb+rcb(i)*rcb(i)
20 continue
zwa=ga2*prcaa-dfloat(2*lmni(ia)+3)
zwb=gb2*prcbb-dfloat(2*lmni(ib)+3)
c
c calculation of the overlap-integrals
c
if(lodd.or.modd.or.nodd) go to 110
o( 1)=ola(l ,m ,n )
o( 5)=ola(l+2,m ,n )
o( 6)=ola(l ,m+2,n )
o( 7)=ola(l ,m ,n+2)
o(11)=ola(l+4,m ,n )
o(12)=ola(l ,m+4,n )
o(13)=ola(l ,m ,n+4)
o(20)=ola(l+2,m+2,n )
o(21)=ola(l+2,m ,n+2)
o(22)=ola(l ,m+2,n+2)
110 continue
c
if(leven.or.modd.or.nodd) go to 120
o( 2)=ola(l+1,m ,n )
o( 8)=ola(l+3,m ,n )
o(14)=ola(l+1,m+2,n )
o(15)=ola(l+1,m ,n+2)
120 continue
c
if(lodd.or.meven.or.nodd) go to 130
o( 3)=ola(l ,m+1,n )
o( 9)=ola(l ,m+3,n )
o(16)=ola(l+2,m+1,n )
o(17)=ola(l ,m+1,n+2)
130 continue
c
if(lodd.or.modd.or.neven) go to 140
o( 4)=ola(l ,m ,n+1)
o(10)=ola(l ,m ,n+3)
o(18)=ola(l+2,m ,n+1)
o(19)=ola(l ,m+2,n+1)
140 continue
c
if(leven.or.meven.or.nodd) go to 150
o(23)=ola(l+1,m+1,n )
150 continue
if(leven.or.modd.or.neven) go to 160
o(24)=ola(l+1,m ,n+1)
160 continue
if(lodd.or.meven.or.neven) go to 170
o(25)=ola(l ,m+1,n+1)
170 continue
c
c calculation of kinetic energy
c
p2=-ga*(zwa*o(1)
& +ga2*(a2*(rca(1)*o(2)+rca(2)*o(3)+rca(3)*o(4))
& +o(5)+o(6)+o(7)))
c
c calculation of relativistic correction
c
p4=(gab4*(o(11)+o(12)+o(13)+a2*(o(20)+o(21)+o(22))
& +a2*(srcab(1)*(o( 8)+o(16)+o(18))
& +srcab(2)*(o(14)+o( 9)+o(19))
& +srcab(3)*(o(15)+o(17)+o(10)))
& +a4*(rca(1)*rcb(1)*o(5)
& +rca(2)*rcb(2)*o(6)
& +rca(3)*rcb(3)*o(7)
& +a2*((rca(1)*rcb(2)+rca(2)+rcb(1))*o(23)
& +(rca(1)*rcb(3)+rca(3)+rcb(1))*o(24)
& +(rca(2)*rcb(3)+rca(3)+rcb(2))*o(25))))
& +(gb2*zwa+ga2*zwb)*(o(5)+o(6)+o(7))
& +a2*(gb2*zwa*(rcb(1)*o(2)+rcb(2)*o(3)+rcb(3)*o(4))
& +ga2*zwb*(rca(1)*o(2)+rca(2)*o(3)+rca(3)*o(4)))
& +zwa*zwb*o(1))
& *(-gab4)/(a8*cl**2)
c
c store p2 and p4 to the correct storage-locations
c
v(1)=p2*ropigc
v(2)=p4*ropigc
v(3)=v(1)+v(2)
return
end
c
c
real*8 function aainer(r,s,t,l,m,n,d,dp,fnu,fn,fd)
implicit real*8(a-h,o-z)
dimension d(3),dp(1),fnu(1),fn(1),fd(1)
integer r,s,t,u,v,w,uvwt,rstt,uvwth
rstt=r+s+t
lmnt=l+m+n
lmnrst=rstt+lmnt
lh=l/2
mh=m/2
nh=n/2
aainer=0.0d0
last1 = lh+1
do 6 ii=1,last1
i = ii-1
il=l-2*i+1
ilr=r+il
ilrh=(ilr-1)/2
fi=fn(ilr)*fd(il)*fd(i+1)
last2 = mh+1
do 5 jj=1,last2
j = jj-1
jm=m-2*j+1
jms=s+jm
jmsh=(jms-1)/2
fij=fn(jms)*fd(jm)*fd(j+1)*fi
last3 = nh+1
do 4 kk=1,last3
k = kk-1
kn=n-2*k+1
knt=t+kn
knth=(knt-1)/2
ijkt=i+j+k
fijk=fn(knt)*fd(kn)*fd(k+1)*dp(ijkt+1)*fij
lmrsij=lmnrst-2*ijkt
last4 = ilrh+1
do 3 iu=1,last4
u = iu-1
ilru=ilr-2*u
fu=fd(u+1)*fd(ilru)
if(abs(d(1))-0.0000001d0)10,10,11
10 if(ilru-1) 12,12, 3
11 fu=fu*d(1)**(ilru-1)
12 last5 = jmsh+1
do 2 iv=1,last5
v = iv-1
jmsv=jms-2*v
fuv=fu*fd(v+1)*fd(jmsv)
if(abs(d(2))-0.0000001d0)20,20,21
20 if(jmsv-1) 22,22,2
21 fuv=fuv*d(2)**(jmsv-1)
22 last6 = knth+1
do 1 iw=1,last6
w = iw-1
kntw=knt-2*w
fuvw=fuv*fd(w+1) *fd(kntw)
if(abs(d(3))-0.0000001d0)30,30,31
30 if(kntw-1) 32,32,1
31 fuvw=fuvw*d(3)**(kntw-1)
32 uvwt=u+v+w
uvwth=uvwt/2
if(2*uvwth-uvwt) 33,34,33
33 fuvw=-fuvw
34 nuindx=lmrsij-uvwt
fuvw=fuvw*fnu(nuindx +1)*dp(uvwt+1)
aainer=fijk*fuvw+aainer
1 continue
2 continue
3 continue
4 continue
5 continue
6 continue
return
end
c
c
subroutine divpt(a,alpha,b,beta,ecoord,eexp,efact)
implicit real*8(a-h,o-z)
c this computes the centre, exponent, and multiplying factor of
c a single gaussian which can replace the product of two gaussia
c centres a and b, and exponents alpha and beta.
dimension a(3), b(3), ecoord(3)
do 1 i=1,3
1 ecoord(i)=(alpha*a(i)+beta*b(i))/(alpha+beta)
eexp=alpha+beta
absqd=0.0d0
do 2 i=1,3
2 absqd=absqd+(b(i)-a(i))*(b(i)-a(i))
efact=exp(-alpha*beta*absqd/(alpha+beta))
return
end
c
c
subroutine elcmom(nt,frocc,dint,nro,norb,tot,cint,nomx)
c total property in mo basis * occupation numbers
implicit real*8(a-h,o-z)
dimension frocc(*),dint(nro,*),tot(*),cint(nomx,*)
do 11 nu=1,nt
tot(nu)=0.0d0
do 12 n=1,norb
tot(nu)=tot(nu) + frocc(n)*cint(n,nu)
dint(2,nu)=tot(nu)
12 continue
11 continue
return
end
c
c
subroutine total(nt,dint,nro)
c nuclear + electronic part of property
implicit real*8(a-h,o-z)
dimension dint(nro,*)
do 20 k=1,nt
20 dint(3,k)=dint(1,k)+dint(2,k)
return
end
c
c
subroutine output(dint,norb,n,nt,nro,frocc,tot,c,ilbl,a,
& title,cint,buff)
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
real*8 c(1), ilbl(1), a(1), title(1)
parameter (mloc=250)
dimension dint(nro,*),frocc(*),tot(nt),cint(nomx,*),buff(lbuff)
c
c dimensions for localization have to be:
c help(ncont,norb,nt), xlv(ncont,norb), op(nnorb,nt), d(norb,norb)
c eig(norb), mon(norb)
c
dimension mon(mloc),xlv(2),repv(2)
logical ok
common /big/ y(nomx*nomx), help(nomx*mloc*3),
& d(mloc*(mloc+1)*3/2),op(mloc*mloc)
common /eig/ eig(mloc)
common /prptyp/ m
common /contrc/ chomo,ncont
equivalence (help(1),xlv(1),repv(1))
c
if(m.eq.15) go to 17
call cntrac(cint,norb,n,nt,nro,buff)
call elcmom(nt,frocc,dint,nro,norb,tot,cint,nomx)
call total(nt,dint,nro)
call prnt(dint,nro,nt,norb,c,ilbl,a,title,cint,nomx)
return
c
17 write(6,6000)
do 10 i=1,norb
10 mon(i)=i
c
read(5,1000) nmo
if(nmo) 300,200,100
100 read(5,1000) (mon(i),i=1,nmo)
c
ok=.true.
ialt=mon(1)
do 220 ip=2,nmo
if(mon(ip).gt.ialt) go to 210
ok=.false.
write(6,2222) ip,mon(ip),ialt
210 ialt=mon(ip)
220 continue
if (ok) go to 250
stop
c
200 read(5,2000) thr
write(6,6002) thr
ii=0
do 20 i=1,norb
if(eig(i).lt.thr) go to 20
ii=ii+1
mon(ii)=i
20 continue
nmo=ii
c
250 ij=0
do 30 i=1,nmo
eig(i)=eig(mon(i))
ist=(mon(i)-1)*ncont+1
ie=ist+ncont-1
do 40 k=ist,ie
ij=ij+1
40 y(ij)=y(k)
30 continue
norb=nmo
c
300 write(6,6001) norb,(mon(i),i=1,norb)
nnorb=norb*(norb+1)/2
call local(nt,norb,ncont,op,nnorb,buff,lbuff,y,help,d,xlv,repv)
return
c
1000 format(5x,5i5)
2000 format(5x,f10.4)
2222 format(' ====> danger : ordering of mos illegal ',i5,5x,2i5)
6000 format(//' boys localization '//)
6001 format(1x,i5,' occupied mo s selected :',5i5/(31x,5i5))
6002 format(6x,' threshold for selection of mo s :',f10.4/)
end
c
c
subroutine local(nt,norb,ncont,op,nnorb,buff,lbuff,y,help,d,xlv,
1 repv)
implicit real*8(a-h,o-z)
integer opn(3)
logical ok
parameter (mloc=250)
dimension y(ncont,norb),help(ncont,norb,nt),buff(lbuff),
& op(nnorb,nt),d(norb,norb),xlv(ncont,norb),repv(2)
common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
common/supp/nsupp
common/eig/eig(mloc)
c
c help(i,j), xlv(i,j) and repv(i) are equivalenced in subroutine
c output
c
rewind 2
do 10 m=1,nt
do 10 i=1,norb
do 10 k=1,ncont
10 help(k,i,m)=0.0d0
do 20 j=1,nt
20 opn(j)=j
c
ix=lbuff
do 30 k=1,ncont
do 30 l=1,k
fac=1.0d0
if(k.eq.l) fac=0.5d0
do 50 m=1,nt
ix=ix+1
if(ix.le.lbuff)go to 40
ix=1
read(2) buff
40 continue
temp=-buff(ix)*fac
do 60 i=1,norb
help(k,i,m)=help(k,i,m)+temp*y(l,i)
60 help(l,i,m)=help(l,i,m)+temp*y(k,i)
50 continue
30 continue
c
do 70 m=1,nt
ij=0
do 80 i=1,norb
do 80 j=1,i
ij=ij+1
top=0.0d0
do 90 k=1,ncont
90 top=top+y(k,i)*help(k,j,m)
op(ij,m)=top
80 continue
70 continue
c
c
read(5,5000) nop
if(nop.eq.0) go to 100
read(5,5000) (opn(nn),nn=1,nop)
c
if(nop.eq.1) go to 130
ok=.true.
ialt=opn(1)
do 110 ip=2,nop
if(opn(ip).gt.ialt) go to 120
ok=.false.
write(6,6000) ip,opn(ip),ialt
120 ialt=opn(ip)
110 continue
if (ok) go to 130
stop
c
130 if(nsupp.lt.2) go to 140
write(6,6010)
do 150 i=1,nnorb
150 write(6,6020) (op(i,kk),kk=1,nt)
c
140 ii=0
do 160 n=1,nop
nopn=opn(n)
do 160 l=1,nnorb
160 op(l,n)=op(l,nopn)
nt=nop
c
100 write(6,6030) nt,(opn(nn),nn=1,nt)
thi=1.0d-13
call boys(norb,nt,nnorb,op,thi,d)
c
ij=0
do 300 i=1,norb
do 300 j=1,i
ij=ij+1
repv(ij)=0.d0
do 300 k=1,norb
300 repv(ij)=repv(ij)+d(k,i)*eig(k)*d(k,j)
write(6,6140)
call outpak(repv,norb,1,2)
c
do 170 i=1,norb
do 180 k=1,ncont
xlv(k,i)=y(k,1)*d(1,i)
if(norb.eq.1) go to 180
do 190 j=2,norb
190 xlv(k,i)=xlv(k,i)+y(k,j)*d(j,i)
180 continue
170 continue
nxlv=norb*ncont
c
if(nsupp.lt.2) go to 200
write(6,6040)
do 210 i=1,nnorb
210 write(6,6020) (op(i,kk),kk=1,nt)
c
200 write(6,6050)
write(6,6060) (opn(i),i=1,nt)
ii=0
do 220 i=1,norb
ii=ii+i
220 write(6,6070) i,(op(ii,kk),kk=1,nt)
c
write(6,6080)
do 230 i=1,norb
230 write(6,6090) i,(d(j,i),j=1,norb)
c
if(nsupp.lt.2) go to 240
write(6,6100)
do 250 i=1,norb
250 write(6,6110) i,(xlv(ll,i),ll=1,ncont)
c
240 write(lomao,6120) (alabel(i1),i1=1,19)
write(lomao,6120) (blabel(i1),i1=1,19)
write(lomao,6130) norb,ncont,nxlv
write(lomao,6120) (fmt(i1),i1=1,4)
write(lomao,fmt) ((xlv(ll,i),ll=1,ncont),i=1,norb)
c
return
5000 format(5x,5i5)
6000 format(' ====> danger : ordering of operators illegal ',i5,5x,2i5)
6010 format(/' operators in mo basis '/)
6020 format(4f20.14)
6030 format(/1x,i5,' operators selected :',10i5)
6040 format(/' operators in lmo basis '/)
6050 format(//' diagonalelements of the operators in lmo basis :'/)
6060 format(8x,'lmo',12x,3(i2,12x)/)
6070 format(5x,i5,5x,3f14.8)
6080 format(/' lmo s as linearcombination of mo s ')
6090 format(/' lmo nb.',i5/(1x,5f14.8))
6100 format(/' lmo vectors '/)
6110 format(' lmo nb.',i5/(1x,4f14.8))
6120 format(19a4)
6130 format(20i4)
6140 format(/' fockmatrix in lmo basis',/)
end
c
c
subroutine boys(n,no,nn1,op,thi,d)
c
c subroutine for boys localization r.ahlrichs 2/85
c to determine vectors |i>, i=1,n
c which maximize
c f = sum(i,j=1,n)sum(k=1,no) (-)**2
c equivalent to maximization of
c g = sum(i=1,n)sum(k=1,no) **2 + **2
c equivalent to minimization of
c h = sum(i**2
c stationarity condition
c sum(k) op(ij,k)*(op(ii,k)-op(jj,k)) = 0 all i.ne.j
c
implicit real*8(a-h,o-z)
dimension d(n,n),op(nn1,no)
c
c input
c n = dimension of lin. space
c no = # of operators in f or g
c op(ij,k) matrix representation of k'th operator
c symmetric and stored triangular,nn1=n(n+1)/2
c thi = rotation threshold (1.-6 to 1.-12 sngl or dbl prc)
c result
c d contains the solution vectors as columns
c op the transformed operators
c
a0=0.d0
a1=1.d0
a2=2.d0
a10=10.d0
aqu=0.25d0
accr=1.d-15
th=dmax1(thi,accr)
c
c set d to unit matrix
c
do 20 i=1,n
do 10 j=1,n
10 d(j,i)=a0
20 d(i,i)=a1
if(n.eq.1) return
c
c opn = operator norm needed for cut off threshold
c
opn=a0
do 22 k=1,no
do 21 i=1,nn1
21 opn=opn+op(i,k)**2
22 continue
if(opn.eq.a0) return
opn=sqrt(opn/dfloat(nn1*no))*accr
ths=aqu
it=0
do 200 isw=1,9999
c
c loop over sweeps
c
thsw=a0
ii=1
ij=2
do 110 i=2,n
ii=ii+i
jj=0
im1=i-1
ip1=i+1
li=ii-i
do 100 j=1,im1
c
c loop over i,j rotations
c maximize sum(k) **2 + **2
c like in jacobi
c
jj=jj+j
c
c get ax , b determining the tan of the rotation angle
c
ax=a0
b=a0
do 30 k=1,no
aux1=op(jj,k)-op(ii,k)
aux2=op(ij,k)
ax=ax+aux1*aux2
30 b=b+aqu*aux1**2-aux2**2
c
c rotation angle x0
c
if (abs(ax)+abs(b).lt.opn) go to 100
x0=aqu*datan2(ax,b)
thsw=dmax1(thsw,abs(x0))
if (abs(x0).lt.ths) go to 100
c
c rotation
c
it=it+1
c c=dcos(x0)
s=dsin(x0)
s2=s*s
c2=a1-s2
c=sqrt(c2)
do 40 l=1,n
aux1=d(l,j)
d(l,j)=c*d(l,j)+s*d(l,i)
40 d(l,i)=c*d(l,i)-s*aux1
jm1=j-1
c c2=c*c
c s2=s*s
cs=c*s
c2ms2=c2-s2
lj=jj-j
jp1=j+1
do 95 k=1,no
aux1=op(ii,k)
aux2=op(jj,k)
op(ii,k)=aux1*c2+aux2*s2-a2*op(ij,k)*cs
op(jj,k)=aux1*s2+aux2*c2+a2*op(ij,k)*cs
op(ij,k)=cs*(aux1-aux2)+op(ij,k)*c2ms2
if(jm1.le.0) go to 60
do 50 l=1,jm1
aux1=c*op(li+l,k)-s*op(lj+l,k)
op(lj+l,k)=s*op(li+l,k)+c*op(lj+l,k)
50 op(li+l,k)=aux1
60 if (im1.eq.j) go to 80
jl=jj+j
do 70 l=jp1,im1
aux1=c*op(li+l,k)-s*op(jl,k)
op(jl,k)=c*op(jl,k)+s*op(li+l,k)
op(li+l,k)=aux1
70 jl=jl+l
80 if (i.eq.n) go to 95
jl=ii+j
il=ii+i
do 90 l=ip1,n
aux1=c*op(jl,k)+s*op(il,k)
op(il,k)=c*op(il,k)-s*op(jl,k)
op(jl,k)=aux1
il=il+l
90 jl=jl+l
95 continue
100 ij=ij+1
110 ij=ij+1
c
c check convergence
c
if (thsw.le.th) go to 210
ths=dmin1(thsw**2,ths*aqu)
ths=dmax1(th,ths)
c
c write (6,199) isw,it,thsw,ths
c199 format (' isw,it,thsw,,ths',2i5,2d16.8)
c
200 continue
write (6,300) isw,thsw
return
210 write (6,310) it,isw,thsw
return
300 format (///' boys not converged in',i6,' sweeps thsw =',
*d16.8/)
310 format (/' boys converged after',i6,' pivots, in ',i6
*,' sweeps, thsw = ',d16.8/)
end
c
c
subroutine mmult(a,b,bs,n)
implicit real*8(a-h,o-z)
dimension a(n,1), b(n,1), bs(n,1)
do 6 i=1,n
do 6 l=1,n
bs(i,l) = 0.0d0
do 3 j=1,n
do 3 k=1,n
bs(i,l) = bs(i,l) + a(j,i)*b(j,k)*a(k,l)
3 continue
6 continue
do 9 i=1,n
do 9 j=1,n
9 b(i,j) = bs(i,j)
return
end
c
c
subroutine prnt(dint,nro,nt,norb,c,ilbl,a,title,cint,nomx)
c
c this subroutine transforms into the principal axis system and
c produces the final output
c
implicit real*8(a-h,o-z)
dimension dint(nro,*),c(3),a(13,3),title(15,10),cint(nomx,*)
real*8 ilbl(12)
dimension pr1(3,3), pr2(3,3), u(3,3)
dimension name(3)
data name/2hx*,2hy*,2hz*/
data audeb/2.541770d0/
common/ioind/isl1, isl2
common/prptyp/mm
common/ncnt/ilnmcm
common/supp/nsupp
c file shout short output, nshout=13
nshout=13
mmm=0
if(mm.eq.13) mmm=2
c the titels for the relat. correction are stored at title(15,j)
c not at title(13,j) !
write(6,108) (ilbl(i), i=1,12)
if(nsupp.lt.1) goto 50
write(6,104) (a(mm,i), i=1,3)
write(6,201)
write(6,300) (title(mm+mmm,i), i=1,nt)
write(6,301)
do 1 i=1,norb
write(6,106)i,i,(cint(i,k),k=1,nt)
1 continue
50 write(6,108) (ilbl(i), i=1,12)
write(6,200) (a(mm,i), i=1,3), ilnmcm, (c(j), j=1,3)
write(nshout,1200) (a(mm,i), i=1,3), ilnmcm, (c(j), j=1,3)
if(mm.eq.13) write(6,109)
write(6,201)
write(nshout,1201)
write(6,300) (title(mm+mmm,i), i=1,nt)
write(nshout,1300) (title(mm+mmm,i), i=1,nt)
write(6,301)
if(mm.ne.13) go to 6000
write(6,6101) (dint(1,l), l=1,nt)
write(6,6102) (dint(2,m), m=1,nt)
write(6,6103) (dint(3,n), n=1,nt)
go to 6001
6000 write(6,101) (dint(1,l), l=1,nt)
write(6,102) (dint(2,m), m=1,nt)
write(6,103) (dint(3,n), n=1,nt)
write(nshout,1101) (dint(1,l), l=1,nt)
write(nshout,1102) (dint(2,m), m=1,nt)
write(nshout,1103) (dint(3,n), n=1,nt)
c output of electrostatic potential on file 21
if(mm.eq.1) write(21,20000) c(1),c(2),c(3),dint(3,1)
20000 format(4f10.6)
6001 continue
if(mm.ne.4) goto 61
sum=0.d0
do 60 i=1,3
60 sum=sum+dint(3,i)**2
sum=sqrt(sum)*audeb
write(6,403) sum
write(nshout,1403) sum
c only the following properties are transformed into the principal
c axis system: field gradient, quadrupole moment, second moment,
c diamagnetic shielding
61 if(mm.eq.3) go to 12
if(mm.eq.5) go to 12
if(mm.eq.11) go to 12
if(mm.eq.12)go to 12
return
12 flag =abs(dint(3,4)) +abs(dint(3,5)) +abs(dint(3,6))
if(flag.lt.3.0d-05) return
if(mm.eq.3) ll=12
if(mm.eq.5) ll=13
if(mm.eq.11) ll=14
if(mm.eq.12) ll=12
if(nsupp.lt.1) goto 51
write(6,108) (ilbl(i), i=1,12)
write(6,104) (a(mm,i), i=1,3)
write(6,400)
write(6,300) (title(ll,i), i=1,nt)
write(6,1400)
write(6,1300) (title(ll,i), i=1,nt)
write(6,301)
51 pr1(1,1) = dint(3,1)
pr1(1,2) = dint(3,4)
pr1(1,3) = dint(3,5)
pr1(2,2) = dint(3,2)
pr1(2,3) = dint(3,6)
pr1(3,3) = dint(3,3)
call hdiag(pr1,3,0,u,nr)
c transform electronic, total, and nuclear parts.
do 14 i=1,3
pr2(1,1)=dint(i,1)
pr2(1,2)=dint(i,4)
pr2(1,3)=dint(i,5)
pr2(2,2)=dint(i,2)
pr2(2,3)=dint(i,6)
pr2(3,3)=dint(i,3)
do 18 ii=1,3
do 18 ij=ii,3
18 pr2(ij,ii)=pr2(ii,ij)
call mmult(u,pr2,pr1,3)
dint(i,1)=pr2(1,1)
dint(i,4)=pr2(1,2)
dint(i,5)=pr2(1,3)
dint(i,2)=pr2(2,2)
dint(i,6)=pr2(2,3)
dint(i,3)=pr2(3,3)
14 continue
c transform mo integrals.
do 22 i=1,norb
pr2(1,1)=cint(i,1)
pr2(1,2)=cint(i,4)
pr2(1,3)=cint(i,5)
pr2(2,2)=cint(i,2)
pr2(2,3)=cint(i,6)
pr2(3,3)=cint(i,3)
do 15 ii=1,3
do 15 ij=ii,3
15 pr2(ij,ii) = pr2(ii,ij)
call mmult(u,pr2,pr1,3)
cint(i,1)=pr2(1,1)
cint(i,4)=pr2(1,2)
cint(i,5)=pr2(1,3)
cint(i,2)=pr2(2,2)
cint(i,6)=pr2(2,3)
cint(i,3)=pr2(3,3)
if(nsupp.lt.1) goto 22
write(6,106)i,i,(cint(i,k),k=1,nt)
22 continue
write(6,108) (ilbl(i), i=1,12)
write(6,200) (a(mm,i), i=1,3), ilnmcm, (c(k), k=1,3)
write(6,400)
write(6,304) (title(ll,i), i=1,nt)
write(6,101) (dint(1,l), l=1,6)
write(6,102) (dint(2,l), l=1,6)
write(6,103) (dint(3,l), l=1,6)
37 write(6,303)
write(6,401)
do 24 i=1,3
24 write(6,202) name(i), (u(j,i), j=1,3)
return
200 format( //,1x, 3a6,/,1x,38hcalculated with reference to the poin
1t,5x,a6,/,1x,11hcoordinates,/, 6x,3hx =,f12.6,5x,3hy =,f12.6,5x,
23hz =,f12.6)
1200 format(1x, 3a6,/,1x,38hcalculated with reference to the point,
15x,a6,/,1x,11hcoordinates,/, 6x,3hx =,f12.6,5x,3hy =,f12.6,5x,
23hz =,f12.6)
201 format(/,1x,42hin atomic units and for global coordinates)
1201 format(1x,42hin atomic units and for global coordinates)
202 format(1x,a2,1x,3d15.6)
300 format(//,10x,7(9x,a6),/,10x,3(9x,a6))
1300 format(10x,7(9x,a6),/,10x,3(9x,a6))
301 format(/)
303 format(///)
304 format(//,12x,7(9x,a6),/,12x,3(9x,a))
400 format(/,1x,49hin atomic units and for the principle axis system)
1400 format(1x,49hin atomic units and for the principle axis system)
401 format(1x,67hthe new coordinate system, (x*,y*,z*), is related to
1the old one by,//,16x,1hx,14x,1hy,14x,1hz,//)
403 format(1x,'norm ',8x,d15.7,' debeye')
1403 format(1x,'norm ',8x,d15.7,' debeye')
101 format(/,1x,7hnuclear,5x,7(1x,d15.7),/,13x,3(1x,d15.7))
102 format(1x,10helectronic,2x,7(1x,d15.7),/,13x,3(1x,d15.7))
103 format(/,1x,5htotal,7x,7(1x,d15.7),/,13x,3(1x,d15.7))
1101 format(1x,7hnuclear,5x,7(1x,d15.7),/,13x,3(1x,d15.7))
1102 format(1x,10helectronic,2x,7(1x,d15.7),/,13x,3(1x,d15.7))
1103 format(1x,5htotal,7x,7(1x,d15.7),/,13x,3(1x,d15.7))
104 format( //,1x,23hintegrals over mo s for,5x,3a6)
106 format(1x,2i5,7(1x,d15.7),/,11x,3(1x,d15.7))
108 format(1h0,12a6)
109 format(/,1x,'the sum of the nuclear part means the darwin-term.',
&' the relativisticly corrected '/' kinetic energy is printed as ',
&'the total sum.')
6101 format(/,1x,7hdarwin ,5x,3(1x,f25.6))
6102 format(1x,10helectronic,2x,3(1x,f25.6))
6103 format(/,1x,5htotal,7x,3(1x,f25.6))
end
c
c
subroutine hdiag(h,n,iegen,u,nr)
c diagonalization of a real symmetric matrix by the jacobi method
c where h is the array to be diagonalized.
c n is the order of the matrix, h.
c returns eigval in x
c iegen must be set unequal to zero if only eigenvalues are
c to be computed.
c iegen must be set equal to zero if eigenvalues and eigenvectors
c are to be computed.
c u is the unitary matrix used for formation of the eigenvectors.
c nr is the number of rotations.
c the subroutine operates only on the elements of h that are to the
c right of the main diagonal.
implicit real*8(a-h,o-z)
dimension h(3,3), u(3,3), x(3), iq(3)
9 if (iegen) 15,10,15
10 do 14 i=1,n
do 14 j=i,n
if(i-j)12,11,12
11 u(i,j)=1.0d0
go to 14
12 u(i,j)=0.d0
u(j,i) = 0.d0
14 continue
15 nr = 0
if (n-1) 2001,2001,17
c scan for largest off diagonal element in each row
c x(i) contains largest element in ith row
c iq(i) holds second subscript defining position of element
17 nmi1=n-1
hdtest = 0.0d0
do 30 i=1,nmi1
x(i) = 0.d0
ipl1=i+1
do 30 j=ipl1,n
hdtest = hdtest + h(i,j)**2
if (x(i).gt.abs(h(i,j))) go to 30
20 x(i) =abs(h(i,j))
iq(i)=j
30 continue
hdtest = 2.0d0*hdtest
c set indicator for shut-off.rap=2**-27,nr=no. of rotations
rap=1.d-12
do 35 i=1,n
35 hdtest = hdtest + h(i,i)**2
hdtest =rap*(sqrt(hdtest/dfloat(n)) + 1.0d0)
c rap*(rms eigval + 1)
c find maximum of x(i) s for pivot element and
c test for end of problem
40 do 70 i=1,nmi1
if (i-1) 60,60,45
45 if ( xmax- x(i)) 60,70,70
60 xmax=x(i)
ipiv=i
jpiv=iq(i)
70 continue
c return if max.h(i,j)less than(2**-27)absf(h(k,k)-min)
if (hdtest- xmax) 148,1000,1000
148 nr = nr+1
c compute tangent, sine and cosine,h(i,i),h(j,j)
150 hhdifo=0.5d0*(h(ipiv,ipiv)-h(jpiv,jpiv))
hhdifn=dsign(1.0d0,hhdifo)*sqrt(hhdifo**2+h(ipiv,jpiv)**2)
cosine=sqrt(0.5d0*(hhdifo+hhdifn)/hhdifn)
sine=h(ipiv,jpiv)/hhdifn*0.5d0/cosine
hhsum=0.5d0*(h(ipiv,ipiv)+h(jpiv,jpiv))
h(ipiv,ipiv)=hhsum+hhdifn
h(jpiv,jpiv)=hhsum-hhdifn
h(ipiv,jpiv)=0.0d0
c inspect the iqs between i+1 and n-1 to determine
c whether a new maximum value should be computed since
c the present maximum is in the i or j row.
do 350 i=1,nmi1
if(i-ipiv)210,350,200
200 if(i-jpiv)210,350,210
210 if(iq(i)-ipiv)230,240,230
230 if(iq(i)-jpiv)350,240,350
240 htempi = h(i,ipiv)
htempj = h(i,jpiv)
h(i,ipiv) = 0.0d0
h(i,jpiv) = 0.0d0
ipl1=i+1
x(i) =0.d0
c search in depleted row for new maximum
do 320 j=ipl1,n
if (x(i).ge.abs(h(i,j))) go to 320
300 x(i) =abs(h(i,j))
iq(i)=j
320 continue
h(i,ipiv)=htempi
h(i,jpiv)=htempj
350 continue
x(jpiv) =0.d0
x(ipiv) =0.d0
c change the other elements of h
do 530 i=1,n
if(i-ipiv)370,530,420
370 htemp = h(i,ipiv)
h(i,ipiv) = cosine*htemp + sine*h(i,jpiv)
if (x(i).ge.abs(h(i,ipiv))) go to 390
x(i) =abs(h(i,ipiv))
iq(i) = ipiv
390 h(i,jpiv) = -sine*htemp + cosine*h(i,jpiv)
if (x(i).ge.abs(h(i,jpiv))) go to 530
400 x(i) =abs(h(i,jpiv))
iq(i) = jpiv
go to 530
420 if(i-jpiv)430,530,480
430 htemp = h(ipiv,i)
h(ipiv,i) = cosine*htemp + sine*h(i,jpiv)
if (x(ipiv).ge.abs(h(ipiv,i))) go to 450
440 x(ipiv) =abs(h(ipiv,i))
iq(ipiv) = i
450 h(i,jpiv) = -sine*htemp + cosine*h(i,jpiv)
if (x(i)-abs(h(i,jpiv))) 400,530,530
480 htemp = h(ipiv,i)
h(ipiv,i) = cosine*htemp + sine*h(jpiv,i)
if (x(ipiv).ge.abs(h(ipiv,i))) go to 500
x(ipiv) =abs(h(ipiv,i))
iq(ipiv) = i
500 h(jpiv,i) = -sine*htemp + cosine*h(jpiv,i)
if (x(jpiv).ge.abs(h(jpiv,i))) go to 530
x(jpiv) =abs(h(jpiv,i))
iq(jpiv) = i
530 continue
c test for computation of eigenvectors
if(iegen)40,540,40
540 do 550 i=1,n
htemp=u(i,ipiv)
u(i,ipiv)=cosine*htemp+sine*u(i,jpiv)
550 u(i,jpiv)=-sine*htemp+cosine*u(i,jpiv)
go to 40
1000 continue
call schort(u,n)
2000 return
2001 x(1) = h(1,1)
go to 2000
end
c
c
subroutine schort(u,n)
implicit real*8(a-h,o-z)
dimension u(3,3)
do 605 i=1,n
if (i.eq.1) go to 601
im1 = i-1
do 603 j=1,im1
term = 0.0d0
do 602 k=1,n
602 term = term + u(k,i)*u(k,j)
do 603 k=1,n
603 u(k,i) = u(k,i) - term*u(k,j)
601 term = 0.0d0
do 604 k=1,n
604 term = term + u(k,i)**2
term = 1.0d0/sqrt(term)
do 605 k=1,n
605 u(k,i) = u(k,i)*term
return
end
subroutine inpsym(mcu,msu,mru,mctu,mcsu,icomx,mgcsu)
c
c this subroutine reads the data from the karlsruhe integral input
c program (file ft24)
c
implicit real*8(a-h,o-z)
parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
& lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
& lpru=18000,lsu=50,lstu=8,lconu=12,
& lcu=8,lccu=182)
common/supp/nsupp
common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
1 ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
2 eta(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
3 z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
4 la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
5 mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu),iic(lctu,lcsu)
dimension blabel(19)
dimension ipq(256),lfmt(2)
data a0,a1/0.0d0,1.0d0/
itape=24
c file aoinp input for ao basis
c itape=24
open(unit=itape,file='argosin',status='old',
1access='sequential',form='formatted')
ntape=0
inp=itape
inp2=31
c file geomfl
c inp2=31
** open(unit=inp2,file='geomfl',status='unknown',
c** 1access='sequential',form='formatted')
c** rewind inp2
inpn=25
c file inpn is a copy of aoinp where the general
c contractions are resolved into segmented ones
c** open(unit=inpn,file='aoinpn',status='unknown',
c** 1access='sequential',form='formatted')
call inpnew(mcu,msu,mru,mctu,mcsu,icomx,mgcsu,inp,inp2,inpn)
itape=inpn
inp=itape
rewind itape
ipq(1)=0
do 1000 i=1,255
1000 ipq(i+1)=ipq(i)+i
read(itape,890)blabel
c** read(itape,8900) (lfmt(i),i=1,2)
8900 format(19a4)
7900 format(1x,19a4)
read(itape,*)ngen,ns,naords,ncons,ngcs,itol,icut,idummy,idummy,
c & icp, nkoord
1 itgeom
c** if(itgeom.gt.0) read(inp2,999) niter
c** if(itgeom.gt.0) inp=inp2
if(ns.le.msu) goto 2
write(6,970) ns
stop
2 if(naords.le.kaords) goto 4
write(6,971) naords
stop
4 if(ncons.le.mconsu) goto 6
write(6,972) ncons
stop
6 if(ngcs.le.mgcsu) goto 8
write(6,973) ngcs
stop
8 if(nsupp.eq.2) write (6,790) blabel
c if (nkoord .eq. 1) write(6,8770)
c8770 format(/'atomic coordinates were read from file 28'/)
c** if(nsupp.eq.2) write (6,7900) (lfmt(i),i=1,2)
if(nsupp.eq.2) write (6,800) ngen, ns, naords, ncons, ngcs, itol,
1 icut, idummy, nsupp, itgeom
read(itape,901)nst,(nd(ist),ityp(ist),ist=1,nst)
if(nsupp.eq.2) write (6,801) nst, (nd(ist), ityp(ist), ist=1,nst)
if(nst.le.mstu) goto 9
write(6,974) nst
stop
9 read(itape,*) idmx
if(idmx.eq.0) goto 11
do 1346 ixj=1,idmx
1346 read(itape,900)
11 do 10 ist=1,nst
10 nso(ist)=0
do 20 iaords=1,naords
read(itape,*)iru,(la(ir,iaords),ir=1,iru)
if(iru.le.mru) goto 20
write(6,975) iru
stop
20 nir(iaords)=iru
do 30 igcs=1,ngcs
read(itape,*)icsu,ictu,maords(igcs)
if(nsupp.eq.2) write (6,800) icsu, ictu, maords(igcs)
if(icsu.le.mcsu) goto 26
write(6,976) icsu
stop
26 if(ictu.le.mctu) goto 28
write(6,977) ictu
stop
28 do 30 ics=1,icsu
read(itape,*)(iic(ict,ics),ict=1,ictu)
if(nsupp.eq.2) write (6,900) (iic(ict,ics), ict=1,ictu)
do 30 ict=1,ictu
c1=dfloat(iic(ict,ics))
c(ict,ics,igcs)=c1
ccc=abs(c1)
30 if(ccc.ne.a0.and.ccc.ne.a1)c(ict,ics,igcs)=dsign(sqrt(ccc),c1)
c read in exponents and contraction coefficients
do 40 icons=1,ncons
read (itape,'(i4,i3)') iconu,lmnp1(icons)
read(itape,*)(zet(icon,icons), eta(icon,icons), icon=1,iconu)
if(iconu.le.icomx) goto 55
write(6,2000) icons,icomx
2000 format(1x,'too many primitive functions in contraction nr. ',i5,
1 ' icomx = ',i5)
stop
55 continue
do 39 icon=1,iconu
39 if(eta(icon,icons).eq.0.d0.and.iconu.eq.1)
1eta(icon,icons)=1.d0
if(nsupp.eq.2) write (6,803) iconu,lmnp1(icons),(zet(icon,icons),
1 eta(icon,icons), icon=1,iconu)
40 ncon(icons)=iconu
ngenp1=ngen+1
ng=ngenp1
c read in atomic labels, charges, coordinates, and basis sets
isf=0
isfr=0
do 75 is=1,ns
read(inp,904)mtype(is),nf(is),nc(is),chg(is)
if(nsupp.eq.2) write (6,804) mtype(is), nf(is), nc(is), chg(is)
icu=nc(is)
if(icu.le.mcu) goto 46
write(6,979) icu
stop
46 read(inp,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
c if (nkoord .eq. 1)
c & read(28,*)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
if(nsupp.eq.2)write(6,805)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
if(icu.eq.1) go to 60
do 50 ig=2,ng
read(itape,*)(ica(ic,is,ig),ic=1,icu)
50 if(nsupp.eq.2) write (6,800) (ica(ic,is,ig), ic=1,icu)
60 ifu=nf(is)
if(ifu.eq.0) go to 75
do 70 if=1,ifu
isf=isf+1
read(itape,*)mcons(isf),igcs
if(nsupp.eq.2) write (6,800) mcons(isf), igcs
iaords=maords(igcs)
iru= nir(iaords)
mgcs(isf)=igcs
ilmnp1=lmnp1(mcons(isf))
nnt(isf)=ipq(ilmnp1+1)
ntu(isf)=(nnt(isf)*(ilmnp1+2))/3
ntl(isf)=ntu(isf)-nnt(isf)+1
nct(isf)=icu*nnt(isf)
do 70 ir=1,iru
isfr=isfr+1
lai=la(ir,iaords)
nso(lai)=nso(lai)+1
70 continue
75 continue
790 format(1x,19a4)
970 format(14h msu too small,i3)
971 format(17h kaords too small,i3)
972 format(17h mconsu too small,i4)
973 format(16h mgcsu too small,i4)
974 format(15h mstu too small,i3)
975 format(14h mru too small,i3)
976 format(15h mcsu too small,i3)
977 format(15h mctu too small,i3)
978 format(16h mconu too small,i3)
979 format(14h mcu too small,i3)
800 format(1x,26i3)
801 format(i4,12(i3,a3))
802 format(1x,26f3.0)
803 format(i4,i3/(1x,2g14.8))
804 format(1x,a5,2i3,f3.0)
805 format(1x,3f14.8)
890 format(19a4)
900 format(26i3)
901 format(i3,12(i3,a3))
902 format(26f3.0)
903 format(2i3/(2g14.8))
904 format(a3,2i3,f3.0)
905 format(3f14.8)
949 format(6(/), /////68h0primitive a
1o integrals neglected if exponential factor below 10**(-,i2,1h)/63
2h contracted ao and so integrals neglected if value below 10**(-,i
32,1h)/43h symmetry orbital integrals written on unit,i3)
951 format(11h degeneracy,16x,12(i4,4x))
952 format(6h label,21x,12(2x,a4,2x))
953 format(1h1,12x,a4,6h atoms//15h0nuclear charge,f7.2//7h0center,12x
1,1hx,15x,1hy,15x,1hz/(i4,7x,3f16.8))
954 format(/9h0operator,6x,19hcenter interchanges/2x,6h1(id.),6x,24i3)
955 format(i3,6h(gen.),5x,24i3)
956 format(i3,11x,24i3)
958 format(//10x,a2,9h orbitals/18h0orbital exponents,7x,7(1pg15.7))
959 format(25h contraction coefficients,7(1pg15.7))
960 format(10h0sym. orb.,12x,15hcenter, ao type/14x,(18(i3,1h,,a2)))
961 format(i3,a4,i1,6x,(18f6.2))
999 format(20i4)
return
end
subroutine inpnew(mcu, msu, mru, mctu, mcsu, icomx, mgcsu,
1 inp,inp2,inpn)
c
c this subroutine reads the data from the karlsruhe integral input
c program (file ft24) and makes a segmented contraction copy on
c file inpn=25
c
c irun=0 first scan through ao-input-file
c irun=1 create new ao-input-file (inpn)
c
implicit real*8(a-h,o-z)
parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
& lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
& lpru=18000,lsu=50,lstu=8,lconu=12,
& lcu=8,lccu=182,mrcru=8,ione=1,izero=0)
common/supp/nsupp
common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
1 ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
2 etax(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
3 z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
4 la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
5 mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu),iic(lctu,lcsu)
dimension blabel(19), nrcr(lconsu), icnoff(lconsu)
dimension ipq(256),lfmt(2),eta(mrcru,lsu,lconsu)
dimension nfn(lsu)
data a0,a1/0.0d0,1.0d0/
irun=0
ipq(1)=0
do 1000 i=1,255
1000 ipq(i+1)=ipq(i)+i
itape=inp
1100 continue
read(itape,890)blabel
c** read(itape,8900) (lfmt(i),i=1,2)
8900 format(19a4)
7900 format(1x,19a4)
read(itape,*)ngen,ns,naords,ncons,ngcs,itol,icut,idummy,idummy,
c & icp, nkoord
1 itgeom
c** if(itgeom.gt.0) read(inp2,999) niter
c** if(itgeom.gt.0) inp=inp2
if(ns.le.msu) goto 2
write(6,970) ns
stop
2 if(naords.le.kaords) goto 4
write(6,971) naords
stop
4 if(ncons.le.mconsu) goto 6
write(6,972) ncons
stop
6 if(ngcs.le.mgcsu) goto 8
write(6,973) ngcs
stop
8 if(irun.gt.0)write (inpn,790) blabel
c if (nkoord .eq. 1) write(inpn,8770)
c8770 format(/'atomic coordinates were read from file 28'/)
read(itape,901)nst,(nd(ist),ityp(ist),ist=1,nst)
if(irun.gt.0)then
c** write (inpn,7900) (lfmt(i),i=1,2)
write (inpn,900) ngen, ns, naords, nconsn, ngcs, itol,
1 icut, idummy, idummy, izero
write (inpn,901) nst, (nd(ist), ityp(ist), ist=1,nst)
endif
if(nst.le.mstu) goto 9
write(6,974) nst
stop
9 read(itape,*) idmx
if(irun.gt.0)write(inpn,900) idmx
if(idmx.eq.0) goto 11
do 1346 ixj=1,idmx
if(irun.gt.0) write(inpn,900)
1346 read(itape,900)
11 do 10 ist=1,nst
10 nso(ist)=0
do 20 iaords=1,naords
read(itape,*)iru,(la(ir,iaords),ir=1,iru)
if(irun.gt.0) write(inpn,900)iru,(la(ir,iaords),ir=1,iru)
if(iru.le.mru) goto 20
write(6,975) iru
stop
20 nir(iaords)=iru
do 30 igcs=1,ngcs
read(itape,*)icsu,ictu,maords(igcs)
if(irun.gt.0) write (inpn,900) icsu, ictu, maords(igcs)
if(icsu.le.mcsu) goto 26
write(6,976) icsu
stop
26 if(ictu.le.mctu) goto 28
write(6,977) ictu
stop
28 do 30 ics=1,icsu
read(itape,*)(iic(ict,ics),ict=1,ictu)
if(irun.gt.0) write (inpn,900) (iic(ict,ics), ict=1,ictu)
do 30 ict=1,ictu
c1=dfloat(iic(ict,ics))
c(ict,ics,igcs)=c1
ccc=abs(c1)
30 if(ccc.ne.a0.and.ccc.ne.a1)c(ict,ics,igcs)=dsign(sqrt(ccc),c1)
c read in exponents and contraction coefficients
nconsn=0
do 41 icons=1,ncons
c read (itape,*) iconu, lmnp1(icons), ircru
read (itape,*) iconu, lmnp1(icons), ircru
903 format(3i3)
if(ircru.eq.0) ircru=1
c write (6,900) iconu, lmnp1(icons), ircru
if(ircru.gt.mrcru)then
write (6,983) ircru
stop 'program error'
endif
if(iconu.gt.icomx)then
write (6,982) iconu
stop 'program error'
endif
nrcr(icons)=ircru
ncon(icons)=iconu
c
do 40 icon=1,iconu
read(itape,*)zet(icon,icons),(eta(ircr,icon,icons),ircr=1,ircru)
c write (6,84) zet(icon,icons), (eta(ircr,icon,icons),
c 1 ircr=1,ircru)
40 continue
icnoff(icons)=nconsn
do 100 ircr=1,ircru
nconsn=nconsn+1
if(irun.gt.0) write(inpn,900)iconu,lmnp1(icons),ione
do 39 icon=1,iconu
if(eta(ircr,icon,icons).eq.0.d0.and.ircru.eq.1.and.iconu.eq.1)
1eta(ircr,icon,icons)=1.d0
c write(6,2000)ircr,icon,icons,eta(ircr,icon,icons)
if(irun.gt.0) write(inpn,806)zet(icon,icons),eta(ircr,icon,icons)
39 continue
100 continue
nrcr(icons)=ircru
41 continue
ngenp1=ngen+1
ng=ngenp1
c read in atomic labels, charges, coordinates, and basis sets
isf=0
isfr=0
do 75 is=1,ns
if(irun.eq.0) nfn(is)=0
read(inp,904)mtype(is),nf(is),nc(is),chg(is)
if(irun.gt.0) write (inpn,904) mtype(is), nfn(is), nc(is), chg(is)
icu=nc(is)
if(icu.le.mcu) goto 46
write(6,979) icu
stop
46 read(inp,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
c if (nkoord .eq. 1)
c & read(28,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
if(irun.gt.0)write(inpn,905)(x(ic,is),yy(ic,is),z(ic,is),ic=1,icu)
if(icu.eq.1) go to 60
do 50 ig=2,ng
read(itape,*)(ica(ic,is,ig),ic=1,icu)
if(irun.gt.0) write (inpn,900) (ica(ic,is,ig), ic=1,icu)
50 continue
60 ifu=nf(is)
if(ifu.eq.0) go to 75
do 70 if=1,ifu
isf=isf+1
read(itape,*)mcons(isf),igcs
c write(6,900)mcons(isf),igcs
ircru=nrcr(mcons(isf))
irst=icnoff(mcons(isf))
do 200 ircr=1,ircru
if(irun.gt.0) write(inpn,900)irst+ircr,igcs
if(irun.eq.0) then
nfn(is)=nfn(is)+1
endif
200 continue
70 continue
75 continue
if(irun.eq.0) then
irun=1
rewind inp
rewind itape
goto 1100
endif
84 format(1x,5g14.8)
790 format(1x,19a4)
970 format(14h msu too small,i3)
971 format(17h kaords too small,i3)
972 format(17h mconsu too small,i4)
973 format(16h mgcsu too small,i4)
974 format(15h mstu too small,i3)
975 format(14h mru too small,i3)
976 format(15h mcsu too small,i3)
977 format(15h mctu too small,i3)
979 format(14h mcu too small,i3)
982 format(24h symnew: icomx too small,i3)
983 format(24h symnew: mrcru too small,i3)
806 format(2g14.8)
890 format(19a4)
900 format(26i3)
901 format(i3,12(i3,a3))
904 format(a3,2i3,f3.0)
905 format(3f14.8)
949 format(6(/), /////68h0primitive a
1o integrals neglected if exponential factor below 10**(-,i2,1h)/63
2h contracted ao and so integrals neglected if value below 10**(-,i
32,1h)/43h symmetry orbital integrals written on unit,i3)
951 format(11h degeneracy,16x,12(i4,4x))
952 format(6h label,21x,12(2x,a4,2x))
953 format(1h1,12x,a4,6h atoms//15h0nuclear charge,f7.2//7h0center,12x
1,1hx,15x,1hy,15x,1hz/(i4,7x,3f16.8))
954 format(/9h0operator,6x,19hcenter interchanges/2x,6h1(id.),6x,24i3)
955 format(i3,6h(gen.),5x,24i3)
958 format(//10x,a2,9h orbitals/18h0orbital exponents,7x,7(1pg15.7))
959 format(25h contraction coefficients,7(1pg15.7))
960 format(10h0sym. orb.,12x,15hcenter, ao type/14x,(18(i3,1h,,a2)))
999 format(20i4)
c2000 format(3i4,f10.4)
return
end
c
c
subroutine getran(nbfao,frocc,norb,itrnao,idxxt)
c
c this subroutine calculates the symmetry transformtion matrix p
c for the transformation from sao in ao basis
c
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
& lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
& lpru=18000,lsu=50,lstu=8,lconu=12,
& lcu=8,lccu=182)
common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
1 mcxu, ngcs, mdum(45), mccu, mblu, ns, ng,nsf,nst
common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
1 ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
2 ccoef(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
3 z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
4 la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
5 mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu)
common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
dimension frocc(1)
dimension itrnao(1),nsot(14),idxxt(1)
isfxo=0
isfp=0
ibfao=0
root3=sqrt(3.0d0)
root15=sqrt(15.d0)
do 160 is=1,ns
icu=nc(is)
ifu=nf(is)
do 140 ic=1,icu
isfx=isfxo
isf=isfp
if(ifu.eq.0)go to 140
do 120 if=1,ifu
isf=isf+1
icons=mcons(isf)
jtype=lmnp1(icons)
nbu=(jtype*(jtype+1))/2
if(nbu*icu.ne.nct(isf))go to 920
isft=isfx+(ic-1)*nbu
isfx=isfx+nbu*icu
do 100 nb=1,nbu
ibfao=ibfao+1
isft=isft+1
itrnao(isft)=ibfao
idxxt(ibfao)=0
if(jtype.le.2) goto 100
if(jtype.eq.3.and.nb.le.3) idxxt(ibfao)=1
if(jtype.eq.4.and.nb.gt.3.and.nb.le.9) idxxt(ibfao)=1
if(jtype.eq.4.and.nb.le.3) idxxt(ibfao)=2
100 continue
120 continue
140 continue
isfxo=isfx
isfp=isfp+ifu
160 continue
nsotx=0
do 180 ist=1,nst
nsot(ist)=nsotx
nsotx=nsotx+nso(ist)
180 continue
nbfso=nsotx
if(nsotx.ne.nbfso)go to 920
do 200 i=1,nbfao
do 200 j=1,nbfso
p(i,j)=0.0d0
200 continue
isf=0
isfr=0
isftx=0
do 260 is=1,ns
ifu=nf(is)
if(ifu.eq.0)go to 260
do 255 if=1,ifu
isf=isf+1
igcs=mgcs(isf)
iaords=maords(igcs)
iru=nir(iaords)
ictu=nct(isf)
ics=0
do 250 ir=1,iru
isfr=isfr+1
isft=isftx
lai=la(ir,iaords)
nsot(lai)=nsot(lai)+1
ibfso=nsot(lai)
iau=nd(lai)
if(iau.ne.1)go to 920
do 250 ia=1,iau
ics=ics+1
do 248 ict=1,ictu
isft=isft+1
c1=c(ict,ics,igcs)
ibfao=itrnao(isft)
c
c integral routine (r.pitzer) and property package use different
c normalizations of d(xx),d(yy),d(zz),f(xxy),f(xxz),f(yyx),f(yyz),
c f(zzx),f(zzy) -factor=sqrt(3)- , and for f(xxx),f(yyy),f(zzz)
c -factor=sqrt(15)-
c
if(idxxt(ibfao).eq.1) c1=c1*root3
if(idxxt(ibfao).eq.2) c1=c1*root15
p(ibfao,ibfso)=c1
248 continue
250 continue
isftx=isftx+ictu
255 continue
260 continue
call scfinp(nbfao,frocc,norb)
return
920 print 925
925 format(' nct(isf) is not properly set')
stop
end
c
c
subroutine scfinp(nbfao,frocc,nnorb)
c
c this subroutine reads the data as written out by routine outgrd
c of the scf-program (file ft08)
c
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
parameter (mloc=250,lstu=8,lsfru=250)
common/eig/eig(mloc)
common/supp/nsupp
common /big/ y(nomx*nomx),su(nomx*nomx),p(nomx,nomx)
common/scfinf/cm(2,mloc),norbs(lstu),nc(lstu),nlamda(lstu),
1 no(lstu),ityp(lstu)
common /contrc/ chomo,ncont
dimension frocc(1)
common/holl/blabel(19),alabel(19),fmt(4),jl,lomao,iopen,nsym
data ir/ 8/,ncol/10/, ika/2hc(/, a0 /0.0d0/
c number of ireps
mstu=lstu
c lsfru as in integral program
c total number of functions , sum(norb) , lsfru=nomx
msfru=lsfru
c max number of functions/irep , max(norb)
mxbfpi=lsfru
mxcoef=nomx*nomx
c iopen = -1 : ci or mcscf case
c iopen = 0 : closed shell case
c iopen = 1 : open shell case
ipmax=0
icmax=0
idmax=0
if(nsym.le.mstu)go to 41
print 3300,mstu
stop
41 continue
do 106 l=1,nsym
read (ir,1001) ityp(l)
read (ir,6410) lll,norbs(l),nc(l),no(l)
nlamda(l)=lll
norb=norbs(l)
if(norb.eq.0) goto 107
ipmin = ipmax + 1
ipmax = ipmax + norb
if(ipmax.le.msfru)go to 58
write(6,2800) msfru
stop
58 read (ir,6430) (cm(1,kl), kl=ipmin,ipmax)
if(iopen.le.0) go to 52
read (ir,6430) (cm(2,kl), kl=ipmin,ipmax)
go to 56
52 do 54 kl=ipmin,ipmax
54 cm(2,kl)=a0
56 icl=norbs(l)
jcl=nlamda(l)
icmin=icmax+1
icmax = icmax + icl*jcl
if(icmax.le.mxcoef) goto 63
print 2700, mxcoef
stop
63 continue
237 icx = icmin - 1
do 224 ic=1,icl
icxi = icx + 1
icx = icx + jcl
224 read (ir,fmt) (su(jc), jc=icxi,icx)
66 continue
if(iopen.lt.0) goto 67
read(ir,fmt) (eig(ie),ie=ipmin,ipmax)
67 if(nsupp.eq.2) print 1100
if(nsupp.eq.2) print 2000, ityp(l)
do 105 mk = 1, icl, ncol
kin = mk + ncol - 1
if (kin .gt. icl) kin = icl
if(nsupp.eq.2) print 5200,(ika,i, i=mk,kin)
nlow = icmin + jcl*mk - jcl - 1
nhigh = icmin + kin *jcl - 1
do 510 i1 = 1, jcl
nlow = nlow + 1
if(nsupp.eq.2) print 3500, (su(ic), ic=nlow,nhigh,lll)
510 continue
klmin=ipmin+mk-1
klmax=ipmin+kin-1
if(nsupp.eq.2.and.iopen.gt.-1)print 2500,(cm(1,kl),kl=klmin,klmax)
if(nsupp.eq.2.and.iopen.eq.-1)print 2600,(cm(1,kl),kl=klmin,klmax)
if(nsupp.eq.2.and.iopen.gt.0) print 5000,(cm(2,kl),kl=klmin,klmax)
105 continue
goto 106
107 read(ir,6430) dummy
if(iopen.gt.0) read(ir,6430) dummy
read(ir,fmt) (dummy,ie=1,lll)
read(ir,fmt) dummy
106 continue
273 continue
do 109 i=1,ipmax
109 chomo=chomo-cm(1,i)-cm(2,i)
write(6,6700) chomo
print 6500, alabel, blabel
call comtrn(nsym,p,nbfao,nnorb,frocc)
if(ir.ne.5)rewind ir
return
c format statements
1000 format(19a4)
1001 format(a4)
1100 format(47h1symmetry, basis, and configuration information)
2000 format(/1h0,6x,6hirrep ,a4,19x,20horbital coefficients)
2500 format(30h0number of electrons occupying/
1 24h closed shell orbitals =,f6.1,9f11.1)
2600 format(30h0number of electrons occupying/
1 24h natural orbitals =,f8.5,9f11.5)
2700 format(38h0vector components number greater than,i5)
2800 format(1x,' more than ',i5,' mos (lmos), change dimensions in',
1 ' common eig and scfinf',/, 'for dimensions for localization',
2 ' see subroutine output')
3000 format(44h0number of orbitals in a given block exceeds, i4)
3300 format(40h0maximum value of lambda is greater than,i3)
3500 format(22x,10f11.6)
5000 format(30h0number of electrons occupying/
1 22h open shell orbitals =,f8.1,9f11.1)
5200 format(21x,10(6x,a2,i2,1h)))
6410 format(20i4)
6430 format(6f12.9)
6500 format(//,1x,19a4,/,19a4,//)
6700 format(/,' charge of molecule',f6.2)
end
c
c
subroutine comtrn(nsym,p,nbfao,norb,frocc)
c
c this subroutine transforms the eigenvectors from sao in ao basis
c and calculates occupation numbers
c only occupied mos are transformed
c dimension of y = occupied orbitals * total orbitals
c
implicit real*8(a-h,o-z)
parameter(nocmx=50,nbmx=720,nomx=250,lbuff=31375)
parameter (mloc=250,lstu=8)
common/scfinf/cm(2,mloc),norbs(lstu),nc(lstu),nlamda(lstu),
1 no(lstu),ityp(lstu)
common /big/ y(nomx*nomx),su(nomx*nomx),uv(nomx*nomx)
common/supp/nsupp
dimension p(nomx,1),frocc(1),idummy(100)
data idummy/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
1 21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,
2 38,39,40,41,42,43,44,45,46,47,48,49,50,
3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
4 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
5 0,0,0,0,0,0,0,0,0,0/
ipr=0
do 17 i=1,nsym
nn=norbs(i)
17 ipr=ipr+nn*nbfao
maxuv=nomx*nomx
if(ipr.le.maxuv) goto 200
write(6,1000) ipr,maxuv
1000 format(1x,' dimension in comtrn too small, ipr =',i7,' maxuv =',
1 i7)
stop
200 continue
ibfsox=0
max=0
iorb=0
ic=0
delt=1.0d-6
nst=1
if(nsupp.eq.3) write(6,2000)
do 100 lam=1,nsym
nnn=nlamda(lam)
nn=norbs(lam)
if(nn.eq.0) goto 100
do 90 ii=1,nn
iorb=iorb+1
min=max+1
max=max+nnn
do 80 ibfao=1,nbfao
ic=ic+1
ibfso=ibfsox
x=0.0d0
do 70 iq=min,max
ibfso=ibfso+1
x=x+p(ibfao,ibfso)*su(iq)
70 continue
y(ic)=x
80 continue
90 continue
if(nsupp.ne.3) goto 100
write(6,3000) ityp(lam)
call prtblk(y(nst),nbfao,nn,' ao ',' mo ',idummy)
nst=nst+nbfao*nn
100 ibfsox=ibfsox+nnn
c frocc = sum of closed and open shell occupation numbers
norb=iorb
do 140 i=1,norb
140 frocc(i)=cm(1,i)+cm(2,i)
c
2000 format(/,' eigenvectors in ao basis',/)
3000 format(/,' irrep ',a4)
return
end
c
c
subroutine erduw(a,b,na,epslon)
c taken mainly from version 4, aug., 1971, of jacscf, u. of wa.
c matrix diagonalization by the jacobi method.
c a = real symmetric matrix to be diagonalized. it is stored
c by columns with all sub-diagonal elements omitted, so a(i,j) is
c stored as a((j*(j-1))/2+i).
c b = matrix to be multiplied by the matrix of eigenvectors.
c na = dimension of the matrices.
c epslon is the convergence criterion for off-diagonal elements.
c modified nov 82 r.a.
c programm keeps track of non-diagonal norm and adjusts thresh
c dynamically. the eventually quadratic convergence of jacobi is
c exploited to reduce thresh faster at the end.
c
implicit real*8(a-h,o-z)
dimension a(1), b(1)
data a0, a1s2 /0.0d0, 0.5d0/
if(na.eq.1) return
loopc=64
sumnd=a0
sum=a0
ij=1
do 24 i=1,na
do 16 j=1,i
term=a(ij)*a(ij)
sum=sum+term
if(i.eq.j) go to 16
sumnd=sumnd+term
16 ij=ij+1
24 continue
thrshg=sqrt((sum+sumnd)/na)*epslon
small=sumnd*epslon*na*na
32 if (sumnd.lt.small) go to 35
thresh=sqrt(sumnd+sumnd)/na
go to 37
35 thresh=thresh*0.03d0
37 thresh=dmax1(thresh,thrshg)
n=0
ij=2
jj=1
do 112 j=2,na
jj=jj+j
jm1=j-1
ii=0
do 104 i=1,jm1
ii=ii+i
if(abs(a(ij)).lt.thresh) go to 104
n=n+1
sumnd=sumnd-a(ij)*a(ij)
sum=a1s2*(a(jj)+a(ii))
term=a1s2*(a(jj)-a(ii))
amax=dsign(sqrt(term*term+a(ij)*a(ij)),term)
c=sqrt((amax+term)/(amax+amax))
s=a(ij)/(c*(amax+amax))
a(ii)=sum-amax
a(jj)=sum+amax
a(ij)=a0
im1=i-1
if(im1) 40,56,40
40 ki=ii-i
kj=jj-j
do 48 k=1,im1
ki=ki+1
kj=kj+1
term=c*a(ki)-s*a(kj)
a(kj)=s*a(ki)+c*a(kj)
48 a(ki)=term
56 if(jm1.eq.i) go to 72
ip1=i+1
ik=ii+i
kj=ij
do 64 k=ip1,jm1
kj=kj+1
term=c*a(ik)-s*a(kj)
a(kj)=s*a(ik)+c*a(kj)
a(ik)=term
64 ik=ik+k
72 if(j.eq.na) go to 88
jp1=j+1
ik=jj+i
jk=jj+j
do 80 k=jp1,na
term=c*a(ik)-s*a(jk)
a(jk)=s*a(ik)+c*a(jk)
a(ik)=term
ik=ik+k
80 jk=jk+k
88 ki=im1*na
kj=jm1*na
do 96 k=1,na
ki=ki+1
kj=kj+1
term=c*b(ki)-s*b(kj)
b(kj)=s*b(ki)+c*b(kj)
96 b(ki)=term
104 ij=ij+1
112 ij=ij+1
loopc=loopc-1
if(loopc) 120,1024,120
120 if(thresh.le.thrshg.and.n.eq.0) return
go to 32
1024 write (6,2048)
stop
2048 format(12h error erduw)
end
c
c
subroutine outpak (matrix,nrow,nctl,isw)
c...........version = 09/05/73/03
c.......................................................................
c
c outpak prints a real*8 symmetric matrix stored in row-packed lower
c
c triangular form (see diagram below) in formatted form with numbered
c
c rows and columns. the input is as follows:
c
c matrix(*)...........packed matrix
c
c nrow................number of rows to be output
c
c nctl................carriage control flag: 1 for single space,
c 2 for double space,
c 3 for triple space.
c
c the matrix elements are arranged in storage as follows:
c
c
c 2
c 4 5
c 7 8 9
c 11 12 13 14
c 16 17 18 19 20
c 22 23 24 25 26 27
c
c and so on.
c
c outpak is set up to handle 8 columns/page with a 8f15.8 format
c
c for the columns. if a different number of columns is required, change
c
c formats 1000 and 2000, and initialize kcol with the new number of
c
c columns.
c
c author: nelson h.f. beebe, quantum theory project, university of
c florida, gainesville
c
c.......................................................................
real*8 matrix,column,columm
integer begin,asa,blank,ctl
dimension matrix(1),asa(3)
data kcol/7/, column/8h atom nr/, asa/4h , 4h0 , 4h- /,
& blank/4h /, zero/0.d+00/, columm/8h lmo nr/
ctl = blank
if ((nctl.le.3).and.(nctl.gt.0)) ctl = asa(nctl)
c
c last is the last column number in the row currently being printed
c
last = min0(nrow,kcol)
c
c begin is the first column number in the row currently being printed.
c
ncol=1
c.....begin non standard do loop.
begin=1
1050 ncol = 1
if(isw.eq.1) write (6,1000) (column,i,i = begin,last)
if(isw.eq.2) write (6,1000) (columm,i,i = begin,last)
do 40 k = begin,nrow
ktotal = (k*(k-1))/2 + begin - 1
do 10 i = 1,ncol
if (matrix(ktotal+i) .ne. zero) go to 20
10 continue
go to 30
20 iplk=ktotal+ncol
kstart=ktotal+1
if(isw.eq.1) write (6,2000) ctl,k,(matrix(ii),ii=kstart,iplk)
if(isw.eq.2) write (6,3000) ctl,k,(matrix(ii),ii=kstart,iplk)
30 if (k .lt. (begin+kcol-1)) ncol = ncol + 1
40 continue
last = min0(last+kcol,nrow)
begin=begin+ncol
if (begin.le.nrow) go to 1050
1000 format (/14x,6(3x,a8,i4),4x,a6,i4)
2000 format (a1,'atom nr',i4,2x,8f15.8)
3000 format (a1,' lmo nr',i4,2x,8f15.8)
return
end
c
c
subroutine prtblk(c,nbfn,nomx,lrow,lcol,mos)
c....................................................................
c subroutine prtblk prints orbital coefficients by symmetry blocks.
c it accepts as argumemts the coefficient matrix c, the number of
c functions in the symmetry block nbfn, and the number of orbitals
c that have been printed previoulsy nskip.
c.....................................................................
implicit real*8(a-h,o-z)
dimension c(nbfn,nomx),mos(1)
data ncol/8/
jlast=min0(nomx,ncol)
do 400 jstrt=1,nomx,ncol
write(6,100)(lcol,mos(j),j=jstrt,jlast)
100 format(/8x,8(6x,a4,i4,1x))
do 300 i=1,nbfn
c
c check for non-zero element before printing row
c
do 120 j=jstrt,jlast
if(c(i,j).ne.0.d0)go to 150
120 continue
c
c exit from loop means all c's are zero. skip the write statement
c
go to 300
150 continue
write(6,200)lrow,i,(c(i,j),j=jstrt,jlast)
200 format(1x,a4,i4,8f15.6)
300 continue
jlast=min0(jlast+ncol,nomx)
400 continue
return
end
subroutine orbval(eta,nbmx,norbmx,nomx,nstart,nstop,contrk,
1ftype,frocc,c,d,dinx,cin,v,ddsq,kpmax,norbpl,mdens,nocc,y,iplot)
implicit real*8 (a-h,o-z)
character*4 ftype
common/contrc/chomo,ncont
dimension y(2),c(kpmax,3),d(kpmax,3),dinx(kpmax)
dimension cin(kpmax,norbmx),v(kpmax),ddsq(kpmax),a(3),ftype(1)
dimension eta(nbmx,25),frocc(1),
1nstart(1),nstop(1),contrk(nomx,6),nocc(2)
norb1=norbpl+1
do 15 k=1,kpmax
do 15 i=1,norb1
15 cin(k,i)=0.0
do 50 i=1,ncont
istart=nstart(i)
istop=nstop(i)
iprim=0
do 51 k=1,kpmax
51 dinx(k)=0.0
do 40 ii=istart,istop
iprim=iprim+1
do 11 l=1,3
11 a(l)=eta(ii,l)
do 100 k=1,kpmax
ddsq(k)=0.0
do 25 j=1,3
d(k,j)=c(k,j)-a(j)
25 ddsq(k)=ddsq(k)+d(k,j)**2
ddsq(k)=-ddsq(k)*eta(ii,4)
if(ddsq(k).lt.-290.) go to 40
ddsq(k)=exp(ddsq(k))*contrk(i,iprim)
100 v(k)=1.0
if(ftype(ii).eq.'s ') go to 68
if(ftype(ii).eq.'x ') go to 2
if(ftype(ii).eq.'y ') go to 3
if(ftype(ii).eq.'z ') go to 4
if(ftype(ii).eq.'xx ') go to 5
if(ftype(ii).eq.'yy ') go to 6
if(ftype(ii).eq.'zz ') go to 7
if(ftype(ii).eq.'xy ') go to 8
if(ftype(ii).eq.'xz ') go to 9
if(ftype(ii).eq.'yz ') go to 10
2 do 102 k=1,kpmax
102 v(k)=d(k,1)
go to 68
3 do 103 k=1,kpmax
103 v(k)=d(k,2)
go to 68
4 do 104 k=1,kpmax
104 v(k)=d(k,3)
go to 68
5 do 105 k=1,kpmax
105 v(k)=d(k,1)**2
go to 68
6 do 106 k=1,kpmax
106 v(k)=d(k,2)**2
go to 68
7 do 107 k=1,kpmax
107 v(k)=d(k,3)**2
go to 68
8 do 108 k=1,kpmax
108 v(k)=d(k,1)*d(k,2)
go to 68
9 do 109 k=1,kpmax
109 v(k)=d(k,1)*d(k,3)
go to 68
10 do 110 k=1,kpmax
110 v(k)=d(k,2)*d(k,3)
68 do 69 k=1,kpmax
69 dinx(k)=v(k)*ddsq(k)+dinx(k)
40 continue
do 70 nn=1,norbpl
n=nocc(nn)
kk=ncont*(n-1)+i
do 80 k=1,kpmax
80 cin(k,nn)=dinx(k)*y(kk)+cin(k,nn)
70 continue
50 continue
if(mdens.gt.0) go to 150
do 850 n=1,norbpl
c write(iplot) nocc(n)
850 write(iplot,99)(cin(k,n),k=1,kpmax)
99 format(9f8.5)
if(mdens.eq.0) return
150 do 90 n=1,norbpl
do 90 k=1,kpmax
cin(k,n)=cin(k,n)*cin(k,n)
90 cin(k,norb1)=cin(k,norb1)+cin(k,n)*frocc(n)
do 870 n=1,norb1
c write(iplot) nocc(n)
870 write(iplot,99)(cin(k,n),k=1,kpmax)
c do 860 k=1,kpmax
c 860 write(6,861) cin(k,norbpl),(c(k,i),i=1,3)
861 format(2x,4f16.8)
return
end
subroutine plotin(nmaxx,nmaxy,kpmax,plane,norbpl,mdens,
1nocc,norbmx,iplot,iat)
implicit real*8 (a-h,o-z)
dimension plane(2),nocc(2)
common /vec/ vecx(3),vecy(3),center(3)
c
c plane cartesian coordinates for each grid point
c norbpl number of orbitals to be plotted
c mdens =0 orbital values are calculated
c =1 orbital density and sum of densities is calculated
c nocc indices of orbitals to be plotted
c
c vecx,vecy vectors defining plane for which density is
c calculated
c center center of the plane
c dxmax boundary of drawing in x direction
c dymax boundary of drawing in y direction
c nx number of points in x direction
c ny number of points in y direction
c
read(5,100)vecx
100 format(8f10.0)
read(5,100)vecy
read(5,100)center
read(5,105)dxmax,dymax,nx,ny
105 format(2f10.0,2i4)
write(6,200)
200 format('0plot input'//' vectors defining plane of plot'//
1' original'/)
write(6,210)vecx,vecy
210 format(/' x direction: ',3f10.3,' y direction:',3f10.3)
write(6,220)
220 format(' normalized')
sumx=0.
sumy=0.
do 1000 i=1,3
sumx=sumx+vecx(i)**2
sumy=sumy+vecy(i)**2
1000 continue
sumx=sqrt(sumx)
sumy=sqrt(sumy)
do 1010 i=1,3
vecx(i)=vecx(i)/sumx
vecy(i)=vecy(i)/sumy
1010 continue
write(6,210)vecx,vecy
write(6,225)nx,ny
225 format(' number of points: x direction',i4,
1' y direction',i4)
write(6,230)center
230 format(' origin: ',3f10.5)
write(6,240)dxmax,dymax
240 format(/' boundaries: +/- xmax ',f10.5,' +/-ymax ',f10.5/)
if(nx.gt.nmaxx)nx=nmaxx
if(ny.gt.nmaxy)ny=nmaxy
if(mod(nx,2).eq.0)nx=nx-1
if(mod(ny,2).eq.0)ny=ny-1
write(6,250)nx,ny
250 format(' number of grid points: nx ',i4,' ny ',i4)
write(iplot,99) nx,ny
write(iplot,98) dxmax,dymax
99 format(2i4)
98 format(2f10.5)
kpmax=nx*ny
distx=2.*dxmax/(nx-1)
disty=2.*dymax/(ny-1)
call gridxy(plane,nx,ny,distx,disty,vecx,vecy,center,dxmax,dymax)
read(5,260)norbpl,mdens
260 format(20i4)
read(5,260)(nocc(i),i=1,norbpl)
if(mdens.eq.0)write(6,280)norbpl
280 format(1x,i4,' orbitalcontours will be plotted'/
1' orbital numbers:')
if(mdens.eq.1)write(6,285)norbpl
285 format(1x,i4,' orbitaldensities will be plotted'/
1' orbital numbers:')
if(norbpl.le.norbmx)go to 1020
norbpl=norbmx
write(6,270)norbpl
270 format(' number of orbitals to be plotted reduced to',i4)
1020 continue
write(iplot,99) norbpl
write(6,290)(nocc(i),i=1,norbpl)
290 format(1x,30i4)
read(5,260) iat
write(iplot,99) iat
if(iat.eq.0) write(6,300)
300 format(/,' the atoms will be not ploted')
if(iat.eq.1) write(6,301)
301 format(/,' the atoms will be ploted together with the orbitals')
if(iat.eq.2) write(6,302)
302 format(/,' the atoms will be ploted on separat filed')
return
end
subroutine gridxy(plane,nx,ny,distx,disty,vecx,vecy,center)
implicit real*8 (a-h,o-z)
dimension plane(nx,ny,3),vecx(3),vecy(3),center(3),start(3)
nx1=nx/2
ny1=ny/2
do 10 i=1,3
start(i)=center(i)-nx1*distx*vecx(i)-ny1*disty*vecy(i)
10 continue
do 40 k=1,3
do 30 j=1,ny
do 20 i=1,nx
plane(i,j,k)=start(k)+(j-1)*distx*vecx(k)+(i-1)*disty*vecy(k)
20 continue
30 continue
40 continue
return
end
subroutine molcor(iplot)
implicit real*8(a-h,o-z)
parameter(laords=18,lconsu=60,lgcsu=18,lru=20,lcsu=24,
& lctu=30,lsfu=200,lgu=9,lsfru=250,lnsfu=6000,
& lpru=18000,lsu=50,lstu=8,lconu=12,
& lcu=8,lccu=182)
common /ntgr/ kaords, mconsu, mstu, msfu, mgu, msfru, mnsfu, mpru,
1mcxu,ngcs,mdum(45),mccu,mblu,ns,ng,nsf,nst,ncons
common/bl/nf(lsu),nc(lsu),ncon(lconsu),lmnp1(lconsu),
1 ntl(lsfu),ntu(lsfu),nnt(lsfu),mcons(lsfu),ica(lcu,lsu,lgu),
2 eta(lsu,lconsu),zet(lconu,lconsu),x(lcu,lsu),yy(lcu,lsu),
3 z(lcu,lsu),nir(laords),nso(lstu),nd(lstu),nblpr(lstu),
4 la(lru,laords),nct(lsfu),maords(lgcsu),mgcs(lsfu),
5 mtype(lsu),c(lctu,lcsu,lgcsu),chg(lsu),ityp(lstu),iic(lctu,lcsu)
common/vec/ vecx(3),vecy(3),center(3)
parameter(maxatm=100)
dimension xatom(maxatm),yatom(maxatm)
natom=0
do 75 is=1,ns
icu=nc(is)
do 76 ic=1,icu
natom=natom+1
if(natom.gt.maxatm) stop 'maxatm is not eneught'
xatom(natom)=(x(ic,is)-center(1))*vecx(1)
+ +(yy(ic,is)-center(2))*vecx(2)
+ +(z(ic,is)-center(3))*vecx(3)
yatom(natom)=(x(ic,is)-center(1))*vecy(1)
+ +(yy(ic,is)-center(2))*vecy(2)
+ +(z(ic,is)-center(3))*vecy(3)
76 continue
75 continue
write(iplot,99) natom,(xatom(i),yatom(i),i=1,natom)
99 format(i4,/,(9f8.5))
return
end
|
