*deck deck=servec expand=true block data implicit real*8(a-h,o-z) include 'lib/lib.com' c common/discc/ yed(16) c character*4 yed common/titel/ title(10) common /cntrl/ lprint,maxdim,ncore logical lprint data lprint/.true./ data title/9*' ','servec '/ c data yed/'ed0','ed1','ed2','ed3','ed4','ed5','ed6','ed7', c 1 'mt0','mt1','mt2','mt3','mt4','mt5','mt6','mt7'/ end program servec c program servec(input,output,tape5=input,tape6=output) c c c Vector Service Program Manual c c c
c 

The Vector Service Program

c c c

Tabel of Contents

c
c
c
c
c..             atmol vector - service program
c..      -- j.h. v. lenthe  december 1982 / may 1987 / may 1988 ---
c
c***********************************************************************
c                 i n p u t   -  description
c***********************************************************************
c
c-----------------------------------------------------------------------
c
c

Notation Conventions

c c .. istring = string of integer numbers finished by 'end' c 'to' convention is allowed i.e. 2 to 5 = 2 3 4 5 c .. tapnam = atmol file (ed0-ed7,mt0-mt7) c next parameters are iblock (start-block c for dumpfile) and isect (section-number) c .. 'a'/'b' = choice between textstrings 'a' and 'b' c .. 'flop' = text-string (directive is always text) c .. ('aa') = optional text-string c .. ** = after ** comment is given c .. ex. = ex. denotes an example (starting in column 1) c----------------------------------------------------------------------- c c ** note atmol input is free format / maximal 72 columns ** c ** tapnam etc. must be on same card ** c ** istring may go over many cards ** c c----------------------------------------------------------------------- c c

The Directives

c c *************** each directive is immediately obeyed **************** c c c c The TITLE Directive c c title ** give on next line the new title c ex. title c ex. this is the title that will appear on vector-files c c The READ Directive c c read tapnam ('notran') ** read vectors c ** if 'notran' is specified the vectors c ** are not transformed to the original c ** basis / all ctrans (adapt) information c ** is lost however c ex. read ed3 1 1 c c read 'input' ndim nmdim ** read vectors from input freeform c ex read input 3 2 c read 'free' ndim nmdim ** read vectors from input using read * c c read 'format' ** read vectors and occupations etc. formatted c ** for exact form of input *see* the print vectors c ** 'format' directive. this option is used to c ** transport vectors as text-file c c ex. read format c ex. ..................... c ex. ............................ c c c The WRITE Directive c c write tapnam ** write vectors (see read) c c c The INIT Directive c c init tapnam ** produce new (empty) dumpfile c ex. init ed3 1 c c The CONVERT Directive c c convert gamess tapnam ** convert a gamess dumpfile to atmol-format c convert atmol tapnam ** convert an atmol dumpfile to gamess-format c ** mainly vector-sections are affected c ex. convert gamess ed3 1 c c c The S Directive c c s ndim tapnam 'print' ** read overlap integrals ; ndim = dimension c ** if ndim = 0 the current dimension is used c ** if print specified the s-matrix is printed c ex. s 0 ed3 1 c c The SCREEN Directive c c screen 'crit' ** make all s-matrix elements < crit exactly zero c ** if crit is omitted the current criterium is used c c The PRINT Directive c c print 'on'/'off' ** switches print flag c print 'vectors' nv ('occ'/'eig') (npr) ** print first nv vectors c ** default nv = nmdim c ** if 'occ' => give occupations c ** if 'eig' => give orbital energies c ** npr (def.10) = x columns c ex. print off c 2ex. print vectors 7 c c print 'vectors' nv 'format' format ** print vectors formatted c ** according to format(format). nv=0 means nv = nmdim c ** default format = 5f16.9 c ** this option is used to transport vectors as text-file c ** see also *read format* / the output looks as follows c *** 'vectors' ndim nv format c *** ((vc(i,j),i=1,ndim),j=1,nv) (format(format)) c *** 'occ' nn ** nn = 0 for no occupations c *** (occ(i),i=1,nn) (format(format)) c *** 'eig' nn ** nn = 0 for no eigen-values c *** (eig(i),i=1,nn) (format(format)) c *** 'end of format-print' c ** the input for read format should look the same c c ex. print vectors 0 format 5e16.10 c c The SCHMIDT Directive c c schmidt istring ** schmidt orthogonalize the orbitals specified c ex. schmidt 1 5 7 3 end c c schmidt 'set' istring istring ** schmidt first istring onto second c ** (no internal orthogonalisations) c ex. schmidt set 1 4 end 2 3 end c c The LOWDIN Directive c c lowdin istring ** lowdin orthogonalise the orbitals specified c ex. lowdin 2 3 4 5 end c c ** accuracy of schmidt : crit / accuracy of lowdin : crit*10 c c The NORMALISE Directive c c normalise istring ** normalise the orbitals specified c c The CHECK Directive c c check ('print') ** check vectors on orthonormality within crit. c ** if in the overlap-matrix of the vectors an c ** element > crit is found the matrix is printed c ** if 'print' is selected it is printed anyway c c The SDIAG Directive c c sdiag ('set') ('istring') ** diagonalise s-matrix and leave its eigen- c ** vectors as current vector-set, ordered by c ** decreasing eigen-value / needs s-matrix c ** see 's' and 'crit' directives c ** ** if set is specified the s-matrix c ** for that set of unit vectors is diagp- c ** nalised and the eigenvectors are of c ** the original dimension c c The LOCMP Directive c c locmp iset nset ** localise within the set of nset orbitals c ** starting at orbital iset. c ** v. magnasco, a. perico jcp 47,971 (1967) c istring 0 ** a string of numbers finished by a 0 indicates c ** the ao' to localise an mo on. when all mo's c ** to be localised are specified, give 'end' c 'end' ** # mo's to be localised .le. nset c ** an moperm may be useful before applying locmp c ** *** locmp needs s-integrals and vectors so c ** *** the 's' and a 'read' directive must c ex. locmp 3 7 c ex. 1 4 5 0 c ex. 1 2 3 0 c ex. end c c The MOPERM Directive c c moperm istring ** permute mo's c ex. moperm 1 to 10 13 12 end c c The AOPERM Directive c c aoperm istring ** permute ao's c ex. aoperm 17 18 end c c The COMBINE Directive c c combine tapnam tapnam ** identical to atmolscf combine option c ** except that no orthogonalisation is done c ** combine first and second vector - set c ** next cards define the origin of the ao's ** c ifr ito 'a'/'b' ii ** ao's ifr to ito (inclusive) in new set c ... ... ... .. ** are ao's ii onwards of vectorset a/b c ... ... ... .. ** go on until all ao's are assigned c 'end' c ** next cards define the origin of the mo's ** c ifr ito 'a'/'b' ii ** same as for ao's / finish with 'end' c 'end' c ex. combine ed6 1 1 ed4 200 3 c ex. 1 5 a 1 c ex. 6 10 b 1 c ex. 11 20 a 6 c ex. 21 26 b 6 c ex. end c ex. 1 4 a 1 c ex. 5 7 b 1 c ex. 8 18 a 5 c ex. 19 26 b 4 c ex. end c c The EXTRA Directive c c extra istring ** add ao's and corresponding unit-vectors c ** (in existing vectors zeros are inserted) c ** warning : ndim **must** equal nmdim c ** e.g. one performed a dz-calculation on h2 using ao's c ** 1,2 and 3,4 . now one adds polarisation functions (p) c ** on each h-atom (p's will be 3,4,5 and 8,9,10). the start- c ** orbitals for this new calculation are produced (starting c ** from the old 4*4 set) by putting : c ex. extra 3 4 5 8 9 10 end c c The SUMMARY Directive c c summary ** print summary of files and parameters c c The TRANSFORM Directive c c transform tapnam ('notran') ** transform nmdim*ndim vectors in core c ** with a nmdim*nmdim set on tapnam c ** notran as in read c c The OCCUP Directive c c occup nval (occ(i),i=1,nval) ** read occupation numbers (free form) c ex. occup 3 2.0 2.0 1.0 c ** the occupation-numbers will appear on any tape4 or dump- c ** file written , following this directive c c The EIGEN Directive c c eigen nval ** read nval orbital-energies on following c (eig(i),i=1,nval) ** according to format(5e15.5) c ex. eigen 2 c ex. -20.5673456 -11.43456 c c The CRIT Directive c c crit acr npower ** set criterium to acr.10**(-npower) c ** (default 1.0e-14) c ex. crit 1.0 20 c c The NDIM Directive c c ndim nd ** reduce ao-dimension to nd (loose last ao's) c ex. ndim 22 c c The NMDIM Directive c c nmdim nmd ** reduce mo-dimension to nmd (loose last mo's) c c The LOG Directive c c log 'on'/'of' ** switch printflag for printing of directive-cards c ** in interactive sessions log on is more forgiving c c The ADDVEC Directive c c addvec tapnam ** add vectors (same dimension as current vectors) c ** from specified file ,as last orbitals c ex. addvec ed3 1 2 c addvec 'unit' n1 n2 .. .. ** add specified unit vectors c ** (i.e. with 1.0 in positions n1,n2, .. c ** respectively ) c ex. addvec unit 11 7 1 c c The SYMMETRY Directive c c symmetry nao nirr n1 n2 n3 .. ** 'symmetry' selection of mo's c ** nao = number of ao's in each mo to look at (max 8) c ** nirr= number of representations to select (max 12) c ** n1,n2,n3 etc : dimensions of subsets within mo's c ** next cards define selection (nirr cards) ** c (iao(i),ich(i),i=1,nao) ** iao = number of ao c .. .. .. .. .. .. .. .. ** ich = 1,0 or -1 = sign of ao c .. .. .. .. .. .. .. .. ** in mo (parity is ok) c ex. symmetry 3 2 3 5 c ex. 1 1 3 0 12 -1 c ex. 1 0 4 0 10 1 c c symmetry 'h' tapnam ** the symmetry is determined from the c ** (1-electron) h-matrix on 'tapnam' c ** following directives may be issued / finish with 'end' c 'sets' is1 is2 ... ** dimensions of subsets in orbital-space c 'rename' ir1 ir2 ... ** new numbers of the representations c 'end' ** c ex. symmetry h tape3 c ex. order 1 4 2 3 c ex. sets 1 4 8 1 c ex. end c ** this sorts the orbitals for water in a dz basis , keeping c ** the 1s and the 1s* , the occupied and the virtual spaces c ** apart / the a2 symmetry does not occur in the occupied c ** space ,is therefore originally assigned # 4 and is reordered c ** to position 2 c c ******************************************************************* c **** a directive : crit 5 5 : is recommended before symmetry **** c **** selection, since orbitals are usually only accurate to **** c **** 5 figures **** c ******************************************************************* c c The MAXDIM Directive c c maxdim 'low'/'medium'/'high' / set maximum dimension (default low) c ** low : all actions possible c ** medium : no orthogonalisation possible c ** high : and no combine possible c ex. maxdim high c c The STOP / FINISH Directive c c stop/finish ** ends calculation c ex. stop c c*********************************************************************** c c ** e r r o r n u m b e r s ** c

Error Monitoring

c c 10 directive unknown c c 13 filename not recognised c 14 'maxdim' specification not recognised c 15 number of selection ao's > 8 or < 1 .. or .. c number of representations >12 or < 1 'symmetry' c 16 subset specification wrong in 'symmetry' c 17 ao specification wrong in 'symmetry' c 18 total number of mo's in 'schmidt set' > nmdim c 19 no or a too small s-matrix available for orthog.,check,sdiag c 20 wrong permutation c 21 linear dependency in 'schmidt' c 22 linear dependency in 'schmidt set' c c 32 dimension of extra set of vectors not equal to the dimension c of the resident set in 'addvec' c 33 new number of mo's > 255 in 'addvec' c 34 dimension of vector-set read from atmol-file > maxdim c 38 newbas .ne. nbas (contraction used in atmol scf) c 39 error on read format / 'occ' or 'eig' not found c c 40 ndim .ne. nmdim at start of extra c 42 dumpfile name not recognised ('combine') c 43 vector-set dimensions wrong in 'combine' c 44 total dimension > maxdim in 'combine' c 45 dimension of available vector-set not sufficient for localise c c 61 invalid dumpfile c 62 atmol io-error c 63 section number out of range (>190 or <0) c 64 requested section does not exist c 65 incorrect section type c 67 dumpfile starting block < 1 ('combine') c 71 adding section would cause size of dumpfile to exceed maximum c 72 section already exists with length less then new one c c 171 ao specification error in 'combine' c 172 mo specification error in 'combine' c c 201 convert not followed by gamess or atmol c 202 convert dimension vector-set > maxdim c 203 convert conversion problem/inconsistent # elements in block c c*********************************************************************** c
c c c implicit real*8(a-h,o-z) common /cntrl/ lprint,maxdim,ncore logical lprint common // cr(786432) data ntri/6/,kcr/ 1/ c kcr = 1 call getcor(ncore) c c... (nosve) call prep99(ncore) c mm = ncore/ntri maxdim = (dsqrt(1.0d0+8.0d0*mm)-1.0d0)/2.0d0 maxdim = min0(maxdim,255) c call driver (cr(kcr),cr(kcr+maxdim**2), 1 cr(kcr+maxdim**2+maxdim*(maxdim+1)/2), 2 cr(kcr+2*maxdim**2+maxdim*(maxdim+1)/2)) c end c subroutine izero(ii,n) cc... integer zero - set c implicit real*8(a-h,o-z) c dimension ii(n) c do 10 i=1,n c10 ii(i) = 0 c return c end subroutine getcor(ncore) c c... simulate ibm getcor c implicit real*8(a-h,o-z) common // flop(786432) ncore = 196608 return end subroutine driver(vc,s,sq1,tri1) c c... vector-service program driver c implicit real*8(a-h,o-z) dimension vc(*),sq1(*),s(*),tri1(*) dimension comm(6) common /titel/ title(10) c common /cntrl/ lprint,maxdim,ncore common /blkin/ ip(255),ipb(255),flop(35) common /local/ occ(255),hulp(255),eig(255),scrat(255) logical lprint logical lprinl character*8 mon,moff,mormat,mocc,meig,mset,mectors,minput, 1 mblan,mmoper,maoper,munit,mh,mrewin,mretur,notran, 1 matmol,mamess,itest,ites,mfree c data lprinl/.true./,ndims,ndim,nmdim/3*0/,iocc,ieig/2*0/ data crit/1.0d-14/ data comm(4),comm(5)/'servec',' '/ c data mon,moff,mormat,mocc,meig,mset,mectors,minput,mfree, 1 mblan,mmoper,maoper,munit,mh,mrewin,mretur,notran 1 /'on','off','format','occ','eig','set','vectors','input','free', 1 ' ','moperm','aoperm','unit','h','rewind','return','notran'/ data matmol,mamess/'atmol','gamess'/ c maxdio = maxdim c call tidajt(comm(2),comm(3),comm(6),comm(1),ittt) c write(6,2) comm(3),comm(2),comm(1),ittt,maxdim 2 format('1',10x,'atmol vector service program',/, 1 11x,'version may 1988 ',/ 2 11x,'time ',a10,' date ',a10,' jobid ',a10,/, 3 11x,'jobtime',i7,' secs maxdim ',i5,/) sss = cpulft(1) call prep98 write(6,3) sss 3 format(' enter program at ',f8.2,' seconds',/) c 1 call pass(ipp,lprinl) go to (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160, 1 1000,1000,170,180,190,200,210,240,250,270,280, 1 290,300,310,320) 2 ipp c... title 10 read(5,11) title 11 format(10a8) if (lprint) write(6,12) title 12 format(' current title ',10a8,/) go to 1 c... read 20 call inpa(ites) if (ites.ne.minput.and.ites.ne.mfree) go to 21 c... read input call rdvcin(vc,ndim,nmdim,ites) go to 1 c... read format 21 if (ites.ne.mormat) go to 25 call form(vc,ndim,nmdim,occ,eig,iocc,ieig,-1) go to 1 c... read tapnm 25 call goback call tapnm(ifile,iblock,isect,.false.) c... check for notran keyword call inpa(ites) ntra = 0 if (ites.eq.notran) ntra = 1 call getq(vc,eig,occ,ndim,nmdim,ntra,ieig,iocc,isect) go to 1 c... write 30 call tapnm(ifile,iblock,isect,.false.) call putq(comm,title,eig,occ,ndim,ndim,nmdim, 1 ieig,iocc,vc,isect,iddd) go to 1 c... init 40 call tapnm(ifile,iblock,isect,.true.) go to 1 c... s 50 isectd = 192 call inpi(ndims) if (ndims.eq.0) ndims = ndim call tapnm(idump,iblkd,isectd,.false.) c... read s integrals call getsh(s,ndims,1,isectd) c... print ? call inpa(itest) if (itest.eq.'print') then write(6,51) ndims 51 format(/' S-matrix of dimension',i4,' as read ..') call prtrs(s,ndims) end if go to 1 c... print 60 call inpa(itest) if (itest.eq.mon) lprint=.true. if (itest.eq.moff) lprint=.false. if (itest.ne.mectors) go to 63 call inpi(nn) if (nn.le.0) nn = nmdim if = 0 call inpa(itest) if (itest.eq.mormat) go to 61 if (itest.ne.mocc.and.itest.ne.meig) call goback if (itest.eq.mocc) if = iocc if (itest.eq.meig) if = -ieig call inpi(nprvc) if (nprvc.le.0) nprvc = 10 call prvc(vc,ndim,nn,occ,eig,if,10) go to 63 c... print vectors nv 'format' format 61 call form(vc,ndim,nn,occ,eig,iocc,ieig,1) c 63 go to 1 c... schmidt 70 iort=-1 call inpa(itest) if (itest.eq.mset) go to 75 if (itest.ne.mblan) call goback go to 85 c. set 75 iort = 0 call rdistr(ip,n1) call rdistr(ip(n1+1),n2) if (n2.eq.255) n2 = nmdim - n1 if ((n1+n2).gt.nmdim) call error(18) go to 87 c... lowdin 80 iort= 1 c 85 call rdistr(ip,nn) nn = min0(nn,nmdim) 87 if (ndims.lt.ndim) call error(19) call coperm(ip,hulp,nmdim) call pinver(ip,ipb,nmdim) call change(ip,occ,eig,vc,1,ndim,hulp,ndim,nmdim) if (iort.eq.-1) call normvc(vc,s,sq1,ndim,nn,crit) if (iort.eq.0) 1 call schmids(vc,vc(n1*ndim+1),s,hulp,n1,n2,ndim,crit) if (iort.eq.1) call lowdins(vc,s,tri1,sq1,hulp,scrat,nn,ndim,crit) call change(ipb,occ,eig,vc,1,ndim,hulp,ndim,nmdim) go to 1 c... moperm 90 call rdistr(ip,nn) call coperm(ip,hulp,nmdim) if (lprint) write(6,95) mmoper,(ip(i),i=1,nmdim) if (lprint) write(6,999) call change(ip,occ,eig,vc,1,ndim,hulp,ndim,nmdim) go to 1 95 format(' ',a6,'utation :',(t15,25i4)) c... aoperm 100 call rdistr(ip,nn) call coperm(ip,hulp,ndim) if (lprint) write(6,95) maoper,(ip(i),i=1,ndim) if (lprint) write(6,999) call change(ip,scrat,scrat,vc,ndim,1,hulp,nmdim,ndim) go to 1 c... combine 110 nnn = maxdio*(maxdio+3)/2+1 maxsq1 = maxdio**2+maxdio*(maxdio+5)/2 - 4*maxdim call combin(vc,sq1,ndim,nmdim, 1 ip,ipb,occ,eig,hulp, 2 tri1(nnn),tri1(nnn-maxdim),tri1(nnn-2*maxdim), 3 tri1(nnn-3*maxdim),maxsq1) ieig = 0 iocc = 1 do 111 i=1,nmdim 111 occ(i) = 2.0d0 write(6,112) 112 format(' **** all occupations set to 2.0 ****') go to 1 c... extra 120 call extra(vc,occ,eig,hulp,ndim,nmdim) go to 1 c... summary 130 call secsum call whtps write(6,131) lprint,maxdim,ndim,nmdim 131 format(/' lprint ',l5,' maxdim ',i4,' ndim ',i4,' nmdim ',i4,/) go to 1 c... occup 140 call zero(occ,255) call inpi(n) do 141 i=1,n 141 call readf(occ(i)) iocc = 1 if (lprint) write(6,146) (occ(i),i=1,nmdim) if (lprint) write(6,999) go to 1 155 format(5e15.5) 146 format(/' occupations :',(t18,8f13.6)) 147 format(/' eigenvalues :',(t18,8f13.6)) c... eigen 150 call zero(eig,255) call inpi(n) read(5,155) (eig(i),i=1,n) ieig = 1 if (lprint) write(6,147) (eig(i),i=1,nmdim) if (lprint) write(6,999) go to 1 c... crit 160 call inpf(acrit) call inpi(icrit) crit = acrit*10.0d0**(-icrit) if (lprint) write(6,161) crit 161 format(' criterion set to ',e12.5,/) go to 1 c... ndim 170 call inpi(ndimi) ko = ndim + 1 kn = ndimi + 1 do 175 i=2,nmdim call fmove(vc(ko),vc(kn),ndimi) ko = ko + ndim kn = kn + ndimi 175 continue ndim = ndimi go to 1 c... log 180 call inpa(itest) if (itest.eq.mon) lprinl = .true. if (itest.eq.moff) lprinl = .false. go to 1 c... nmdim 190 call inpi(nmdim) if (nmdim.le.0) ndim = 0 go to 1 c... addvec 200 call inpa(itest) if (itest.eq.munit) go to 205 c... addvec from dumpfile or tape4 call goback call tapnm(ifile,iblock,isect,.false.) call getq(vc(nmdim*ndim+1),eig(nmdim+1), 1 occ(nmdim+1),nd2,nmd2,0,ii,ii,isect) if (nd2.ne.ndim) call error(32) nmdim = nmdim + nmd2 if (nmdim.gt.255) call error(33) go to 1 c... addvec / unit vectors 205 call zero(hulp,ndim) call inpi(nn) if (nn.le.0) go to 206 hulp(nn) = 1.0d0 call fmove(hulp,vc(nmdim*ndim+1),ndim) nmdim = nmdim + 1 go to 205 c 206 if (lprint) write(6,209) ndim,nmdim 209 format(' addvec results in ndim=',i4,' / nmdim=',i4,/) go to 1 c... symmetry 210 call inpa(ites) if (ites.ne.mh) go to 215 c... symmetry check on h-matrix isec = 192 call tapnm(ifile,iblk,isec,.false.) call getsh(tri1,ndim,2,isec) call symh(tri1,sq1,vc,occ,eig,hulp,nmdim,ndim,crit,lprint) go to 1 c... symmetry check on vector-coefficients 215 call goback call symmetr(vc,eig,occ,ip,hulp,ipb,ndim,nmdim,crit) go to 1 c... maxdim 240 call smaxd(maxdim,ncore,maxdio) go to 1 c... check 250 if (ndims.lt.ndim) call error(19) call checks(vc,s,tri1,hulp,ndim,nmdim,crit) go to 1 c... locmp 270 if (ndims.lt.ndim) call error(19) call inpi(ill) call inpi(nll) write(6,271) ill,nll 271 format(/' // localise in set starting at mo ',i3, ',of dim. ',i3) if (ill+nll-1.gt.nmdim) call error(45) call locmp(vc(ndim*(ill-1)+1),s,hulp,sq1,nll,ndim,50,crit) go to 1 c... convert 280 call inpa(ites) if (ites.eq.mamess) then call conver(1,sq1,sq1,maxdim) else if (ites.ne.matmol) call error(201) call conver(-1,sq1,sq1,maxdim) end if go to 1 c... sdiag 'set' 'istring' 290 ndim = ndims if (ndim.le.0) call error(19) call inpa(ites) if (ites.eq.mset) then call rdistr(ip,nn) write(6,293) nn 293 format(/' SDIAG SET dimension ',i5) write(6,95) mset,(ip(i),i=1,nn) kk = 0 do 194 i=1,nn do 194 j=1,i kk = kk + 1 imi = max0(ip(i),ip(j)) jmj = min0(ip(i),ip(j)) ks = imi*(imi-1)/2 + jmj tri1(kk) = s(ks) 194 continue else nn = ndim call fmove(s,tri1,ndim*(ndim+1)/2) end if call jacobs(tri1,vc,eig,hulp,nn,crit,3) ieig = 1 iocc = 0 nmdim = ndim write(6,291) (eig(i),i=1,nn) 291 format(//' eigenvalues of s-matrix ',/,(7e11.3)) write(6,292) 292 format(/' eigenvectors of s-matrix are the current vectors ') c c... make nn * ndim vectors if set was specified c if (ites.eq.mset) then nmdim = nn call fmove(vc,sq1,nn*nn) call zero(vc,ndim*ndim) do 195 i=1,nmdim do 195 j=1,nmdim vc(ip(j)+(i-1)*ndim)=sq1(j+(i-1)*nmdim) 195 continue end if go to 1 c... transform 300 call tapnm(ifile,iblock,isect,.false.) c... check for notran keyword call inpa(ites) ntra = 0 if (ites.eq.notran) ntra = 1 call getq(sq1,eig,occ,nd2,nmd2,ntra,ieig,iocc,isect) if (nd2.ne.nmd2) call caserr('transf. matrix must be square') if (nd2.ne.nmdim) call caserr('dim. transform = nmdim') call matmpl(vc,sq1,tri1,nmdim,ndim) ndims = 0 go to 1 c... screen crit 310 if (ndims.eq.0) call caserr('s-matrix must be present for screen') call inpf(acrit) call inpi(icrit) if (acrit.ne.0.0d0) then critsc = acrit*10.0d0**(-icrit) else critsc = crit end if write(6,311) critsc 311 format(/' S-matrix screened for elements < ',1pe12.5) smax = 0.0d0 nmax = 0 do 312 i=1,ndims*(ndims+1)/2 if (dabs(s(i)).lt.critsc) then if (dabs(s(i)).gt.dabs(smax)) smax = s(i) if (s(i).ne.0.0d0) nmax = nmax + 1 s(i) = 0.0d0 end if 312 continue write(6,313) nmax,smax 313 format(' # elements zeroed',i6,' max elemt zeroed ',1pe12.5) go to 1 c... normalise 320 if (ndims.lt.ndim) call caserr('s-matrix needed for normalise') call rdistr(ip,nn) nn = min0(nn,nmdim) call coperm(ip,hulp,nmdim) call pinver(ip,ipb,nmdim) call change(ip,occ,eig,vc,1,ndim,hulp,ndim,nmdim) do 321 i=1,nn 321 call normvc(vc(ndim*(i-1)+1),s,sq1,ndim,1,crit) call change(ipb,occ,eig,vc,1,ndim,hulp,ndim,nmdim) go to 1 c... stop 1000 call clear call shut sss = cpulft(1) write(6,1001) sss 1001 format(' exit program at ',f8.2,' seconds',/,1x) 999 format(1x) c return end subroutine smaxd(maxdim,ncore,mxdio) c c... set new maxdim c implicit real*8(a-h,o-z) dimension name(4) data name/' ','low','medium','high'/ c call inpa(itest) c do 10 i=1,4 if (itest.eq.name(i)) go to 20 10 continue call error(14) c 20 go to (1,1,3,4),i 1 maxdim = mxdio write(6,61) maxdim 61 format(' maxdim restored to old value ',i3,/) return 3 maxdim = mxdio + dsqrt(dble(mxdio)) maxdim = min0(maxdim,255) write(6,63) maxdim 63 format(' maxdim set to',i4,' no orthog possible',/) return 4 maxdim = dsqrt(dble(ncore)) maxdim = min0(maxdim,255) write(6,64) maxdim 64 format(' maxdim set to',i5,' no orthog or combine or trans',/) return end subroutine putq(comm,tit,eig,def,norb,norbn,ncolu,ieig, *ideff,q,mpos,iblkq) c... standard e.vector outputting routine(+ header blocks) implicit real*8(a-h,o-z) dimension q(*),comm(*),tit(*),eig(*),def(*) common/sector/numd,iblkd common/blkin/com(19),title(10),value(255),occ(255), *nbasis,newbas,ncol,ivalue,iocc,ndum c************************************************************* c common/tran/ilifc(256),ntran(256),itran(600),ctran(600), c *iftran c logical iftran c dimension ilifc(511) c equivalence (ilifc(1),value(1)) c equivalence (iftran,value(180)) c*** set iftran to .true. by setting all of value to true c************************************************************* if(mpos.gt.190)call error(63) nbsq=norb*ncolu n544 = (5-1)/nipw() + 1 + 539 n1713 = (1113-1)/nipw() + 1 + 600 c.. j=(nbsq-1)/511+1 + (n544-1)/511+1 + (n1713-1)/511+1 j=(nbsq-1)/511+1 + 2 + (n1713-1)/511+1 call secput(mpos,3,j,iblk) iblkq=iblk+ 2 + (n1713-1)/511+1 call fmove(comm,com,19) call fmove(tit,title,10) call fmove(eig,value,ncolu) call fmove(def,occ,ncolu) nbasis=norb newbas=norbn ncol=ncolu ivalue=ieig iocc=ideff call wrt3(com(1),n544,iblk,numd) j=iblk + 2 c*************************************************************** c... write formally ilifc thoug only iftran=.true. matters n = nipw()*511 call setstl(n,.true.,value) nn = (n1713-1)/511 do 10 i=1,nn 10 call put(value,511,numd) n = n1713 - nn*511 call put(value,n,numd) c**** no contraction **** because only .true.'s written c call wrt3(ilifc(1),n1713,j,numd) c*************************************************************** c call wrt3(q(1),nbsq,iblkq,numd) c call clear call revind c return end subroutine getq(q,eig,occ,nbas,ncolu,ntra,ieig,iocc,j) c servec getq routine / ntra : if 1 => no backtransform of mo's implicit real*8(a-h,o-z) logical iftran,mugg,mug,ipunch,jchek,iprin,nopr dimension q(*),eig(*),occ(*) common /cntrl/ lprint,maxdim,ncore logical lprint common/sector/numd,iblkd common/blkin/com(19),tit(10),deig(255),docc(255), *nba,new,ncol,jeig,jocc common/discc/ ied(16) character*4 ied c if(j.gt.190)call error(63) if (lprint) write(6,100)j,iblkd,ied(numd) 100 format(/4x,'vectors restored from section',i4, *' of dumpfile starting at block',i6,' of ',a4) call secget(j,3,k) n544 = (5-1)/nipw() + 1 + 539 n1713 = (1113-1)/nipw() + 1 + 600 call rdedx(com(1),n544,k,numd) k=k+ (n544-1)/511+1 if (lprint) write(6,200)ncol,nba,(com(7-i),i=1,6),tit nbas=nba newb=new if (nbas.ne.newb) write(6,601) nbas,newb 601 format(' ** note nbas ',i5,' ne newbas ',i5) ncolu = ncol ieig=jeig iocc=jocc call fmove(docc,occ,ncol) call fmove(deig,eig,ncol) if(nbas.gt.maxdim) call error(34) 200 format(/' header block information :'/,i4,' * ',i4, *' vectors created under acct. ',a8/ *a7,'vectors created by ',a8,' atmol3 program at ', *a8,' on ',a8,' in the job ',a8/ *' with the title: ',10a8,/) 500 nw=newb*ncol k=k+ (n1713-1)/511+1 call rdedx(q,nw,k,numd) c if (ntra.ne.1) go to 700 c... no transformation nbas = newb nbas = newb write(6,602) nbas 602 format(/' ** no backtransform : ndim = ',i5) return c... transformation / ctrans or adapt 700 k = k - ( (n1713-1)/511+1 ) call tdownn(q,q(nbas*ncolu+1),ncolu,nbas,newb,k,numd) c return end subroutine tdownn(qnew,q,ncol,nbasis,newbas,iblk,numd) implicit real*8(a-h,o-z) dimension qnew(*),q(*) logical iftran common/blkin/xa(256) common/tran/ilifc(256),ntran(256),itran(600),ctran(600),iftran, * ndum c c... read ctran - info c n1713 = (1113-1)/nipw() + 1 + 600 call rdedx(ilifc,n1713,iblk,numd) c if(iftran) return c c... copy vectors from qnew to q c call fmove(qnew,q,newbas*ncol) c do 731 i=1,ncol m= (i-1)*newbas call zero(xa,nbasis) do 733 j=1,newbas n=ntran(j) do 733 k=1,n l=ilifc(j) +k 733 xa(itran(l))=ctran(l)*q(m+j)+xa(itran(l)) 731 call fmove(xa(1),qnew((i-1)*nbasis+1),nbasis) return end subroutine extra(vc,occ,eig,hulp,ndim,nmdim) implicit real*8(a-h,o-z) dimension vc(*),occ(*),eig(*),hulp(*) common /cntrl/ lprint,maxdim,ncore logical lprint common /blkin/ jug(255),ind1(255) c if (ndim.ne.nmdim) call error(40) call rdistr(jug,mextra) 2 if(mextra)9,8,9 9 if (lprint) write(6,10) 10 format(/' list of extra basis functions') if (lprint) write(6,11)(jug(i),i=1,mextra) 11 format(30i4) if (lprint) write(6,11) goto 12 8 write(6,13) 13 format(/' no extra basis functions specified',/) return c... add unit matrix and permute ao's 12 call extr(vc,vc,nmdim,nmdim+mextra) call izero(ind1,nmdim+mextra) call zero(occ(nmdim+1),mextra) call zero(eig(nmdim+1),mextra) call imove(jug,ind1(nmdim+1),mextra) nmdim = nmdim + mextra ndim = nmdim call coperm(ind1,jug,ndim) call pinver(ind1,jug,ndim) call change(jug,hulp,hulp,vc,ndim,1,ind1,nmdim,ndim) c return end subroutine combin(q,scr,nbasis,newbas, 1 lga,lgb,lgab,ilifa,ilifb,iga,igb,igab,itab, 2 maxsq1) implicit real*8(a-h,o-z) logical iftran,nopr,mugg,mug,iprin,jchek,ipunch logical lga,lgb,lgab dimension iab(2),q(*),scr(*) common/cntrl/ lprint,maxdim,ncore logical lprint dimension lga(255),lgb(255),lgab(255),ilifa(255), *ilifb(255),iga(255),igb(255),igab(255),itab(255) common/work/jrec,jump common/blkin/com(19),tit(10),deig(510),nbb,newb,ncol, *isp(2) c common/disc/ispppp(5),ipos(32) common/discc/ nam(16) character*4 nam,ytest,itest,iab,iend,mole,itab data iab/'a','b'/ data iend/'end'/,mole/'mole'/ n544 = 539 + (5-1)/nipw()+1 n1713 = (1113-1)/nipw() + 1 + 600 kexit=1 c... read in a-set vectors 555 call inpa4(ytest) do 1 ia=1,16 if(ytest.eq.nam(ia))goto 2 1 continue callerror(42) 2 callinpi(la) if(la.lt.1)callerror(67) callsecini(la,ia) call inpi(ja) if(ja.gt.190)callerror(63) if (lprint) write(6,100) iab(kexit),ja,la,ytest 100 format(/7x,a1,'-set vectors restored from section',i4, *' of dumpfile starting at block',i6,' of ',a4) callsecget(ja,3,k) call rdedx(com,n544,k,ia) if (lprint) write(6,200)(com(6-i),i=1,5),tit 200 format(/' header block information:'/ *a7,'vectors created by ',a8,' atmol3 program at ', *a8,' on ',a8,' in the job ',a8/ *' with the title: ',10a8) if(maxdim.lt.nbb.or.ncol.ne.newb) *callerror(43) nw=nbb*ncol k=k+2+ (n1713-1)/511+1 goto (666,777),kexit 666 call rdedx(scr,nw,k,ia) newa=newb nba=nbb nwb=nw c... read in b-set vectors kexit=2 goto555 777 newbas = newa+newb nbasis = nba + nbb if ((nwb+nw).gt.maxsq1) call error(44) call rdedx(scr(nwb+1),nw,k,ia) la=0 do 4 ia=1,newbas ilifa(ia)=(ia-1)*nba ilifb(ia)=(ia-1)*nbb+nwb lga(ia)=.false. lgb(ia)=.false. 4 lgab(ia)=.false. c... process lcbf params 888 callinput callinpa4(itest) if(itest.eq.iend)goto999 c if (itest.ne.mole) go to 889 cc... =molecule= option c call commol(q,nn,newa,newb) c la = la + nn c go to 888 c889 continue jrec=jrec-1 callinpi(m) callinpi(n) callinpa4(itest) callinpi(k) if(m.lt.1.or.n.gt.newbas.or.k.lt.1.or.m.gt.n) *callerror(171) if(itest.eq.iab(2))goto8 if(itest.ne.iab(1))callerror(171) do7ia=m,n if(k.gt.newa.or.lga(k).or.lgab(ia))callerror(171) lga(k)=.true. lgab(ia)=.true. igab(ia)=k itab(ia)=iab(1) iga(k)=ia k=k+1 7 la=la+1 goto888 8 do9 ia=m,n if(k.gt.newb.or.lgb(k).or.lgab(ia))callerror(171) lgb(k)=.true. lgab(ia)=.true. igab(ia)=k itab(ia)=iab(2) igb(k)=ia k=k+1 9 la=la+1 goto888 999 if(la.ne.newbas)callerror(171) if (.not.lprint) go to 43 write(6,300) 300 format(/' origins of lcbf of ab system'/,1x) c nbatch = (newbas-1)/24 + 1 do 10 kk=1,nbatch ke = min0(kk*24,newbas) kb = (kk-1)*24 + 1 write(6,301) (k,k=kb,ke) write(6,302) (itab(k),k=kb,ke) write(6,303) (igab(k),k=kb,ke) 10 continue 301 format(/' ab lcbf ',24i5) 302 format(' subsystem',24(4x,a1)) 303 format(' lcbf ',24i5) c c... construct trial mos 43 nbsq = nbasis**2 la=0 do 44 ia=1,newbas lga(ia)=.false. lgb(ia)=.false. 44 lgab(ia)=.false. c... process mo specification lines callzero(q,nbsq) 8888 callinput callinpa4(itest) if(itest.eq.iend)goto9999 jrec=jrec-1 callinpi(m) callinpi(n) callinpa4(itest) callinpi(k) if(m.lt.1.or.n.gt.newbas.or.k.lt.1.or.m.gt.n) *callerror(172) if(itest.eq.iab(2))goto88 if(itest.ne.iab(1))callerror(172) do77 ia=m,n if(k.gt.newa.or.lga(k).or.lgab(ia))callerror(172) lga(k)=.true. lgab(ia)=.true. igab(ia)=k itab(ia)=iab(1) la=la+1 ja=ilifa(k) nn= (ia-1)*nbasis k=k+1 do77mm=1,newa 77 q(nn+iga(mm))=scr(ja+mm) goto8888 88 do99ia=m,n if(k.gt.newb.or.lgb(k).or.lgab(ia))callerror(172) lgb(k)=.true. lgab(ia)=.true. igab(ia)=k itab(ia)=iab(2) la=la+1 ja=ilifb(k) nn= (ia-1)*nbasis k=k+1 do99mm=1,newb 99 q(nn+igb(mm))=scr(ja+mm) goto8888 9999 if(la.ne.newbas)callerror(172) if (.not.lprint) return write(6,400) 400 format(/' origins of trial mos for ab system'/,1x) c nbatch = (newbas-1)/24 + 1 do 110 kk=1,nbatch ke = min0(kk*24,newbas) kb = (kk-1)*24 + 1 write(6,304) (k,k=kb,ke) write(6,302) (itab(k),k=kb,ke) write(6,305) (igab(k),k=kb,ke) 110 continue 304 format(/' ab mo ',24i5) 305 format(' mo ',24i5) write(6,306) 306 format(1x) c return end subroutine rdvcin(vc,ndim,nmdim,ites) c c... read vectors from input c implicit real*8(a-h,o-z) common /cntrl/ lprint,maxdim,ncore dimension vc(*) character*8 ites,mfree data mfree/'free'/ c call inpi(ndim) call inpi(nmdim) if (ndim.gt.maxdim.or.nmdim.gt.maxdim) call error(35) nn = ndim*nmdim c if (ites.eq.mfree) then read *,(vc(i),i=1,nn) else do 10 i=1,nn 10 call readf(vc(i)) end if c return end subroutine coperm(ip,flop,n) c c make a permutation complete c implicit real*8(a-h,o-z) dimension ip(*),flop(*) logical flop c do 10 i=1,n 10 flop(i) = .false. c c check which numbers are already mentioned c do 20 i=1,n if (ip(i).le.0) go to 20 if (flop(ip(i))) go to 100 flop(ip(i)) = .true. 20 continue c c fill in the missing numbers c k = 0 do 40 i=1,n if (ip(i).gt.0) go to 40 30 k = k + 1 if (flop(k)) go to 30 ip(i) = k 40 continue c return c c error... a number appears twice c 100 write(6,1) (ip(i),i=1,n) 1 format(/1x,'****** a number appears twice in the permutation' x,/,(t10,40(i2,1h,),/)) call error(20) c end subroutine change(idm,dm,em,vc,nspc,nspr,vcin,ndim,nmdim) implicit real*8(a-h,o-z) dimension dm(*),em(*),idm(*),vc(*),vcin(*) c c perform permutations of vectors and occupation-numbers/ eigenval c according to idm / nspc,nspr : column,row-spacing c do 100 k=1,nmdim if (idm(k).eq.k) go to 100 c l = k ll = (l-1)*nspr + 1 c c save vector k c do 50 i=1,ndim 50 vcin(i) = vc(ll+(i-1)*nspc) din = dm(l) ein = em(l) c n = idm(l) 60 nn = (n-1)*nspr + 1 do 70 i=1,ndim 70 vc(ll+(i-1)*nspc) = vc(nn+(i-1)*nspc) dm(l) = dm(n) em(l) = em(n) c idm(l) = l l = n ll = nn n = idm(l) if (idm(n).ne.n) go to 60 c do 80 i=1,ndim 80 vc(ll+(i-1)*nspc) = vcin(i) dm(l) = din em(l) = ein idm(l) = l c 100 continue return end subroutine normvc( vc, s, w, ndim,nmdim, critlow) c c... schmidt othogonalisation of nmdim vectors of length ndim c c.. vc ndim*nmdim c.. s ndim*(ndim+1)/2 c.. w ndim*2 c implicit real*8(a-h,o-z) dimension vc(*),s(*),w(*) c if (nmdim.le.0) return c if (ndim .gt. 1) go to 100 if (vc(1) * vc(1) .le. critlow) goto 160 vc(1) = 1/dsqrt(s(1)) return 100 np1 = ndim+1 ivci = 1 do 150 i=1,nmdim call cntrc(s,vc(ivci),w,ndim,critlow) t = 0 tnorm = 0 ivcj = 1 ivcw = ndim do 120 j=1,i tnorm = tnorm - t*t t = 0 do 110 jj=1,ndim t = vc(ivcj)*w(jj) + t 110 ivcj = ivcj + 1 ivcw = ivcw + 1 120 w(ivcw) = -t if (tnorm + t .le. critlow) goto 160 tnorm = 1/dsqrt(tnorm+t) w(ivcw) = 1 do 140 k = 1,ndim t = 0 ivcw = np1 do 130 j=k,ivci,ndim t = vc(j)*w(ivcw) + t 130 ivcw = ivcw + 1 vc(ivci) = t*tnorm 140 ivci = ivci + 1 150 continue return 160 write (6,620) 620 format(/1x, ' initial vector-set is linearly dependent within thre xshold (critlow)') print 621,i 621 format(/' detected at vector ',i6) call prvc(vc,ndim,nmdim,vc,vc,0,10) call error(21) return end subroutine schmids(v1, v2, s, w, n1, n2, n, critlow) c c schmidt orthogonalization of n1 vectors(v1) of dimension n on c n2 vectors(v2) of dimension n c the resulting vector-set is not internally orthogonalized c results are returned in v1 c c version with metric matrix c implicit real*8(a-h,o-z) dimension v1(*),v2(*),s(*),w(*) c if (n1.le.0.or.n2.le.0) return c k1 = -n do 60 i1=1,n1 k1 = k1 + n call cntrc(s, v1(k1 + 1), w, n, critlow) k2 = -n do 60 i2=1,n2 k2 = k2 + n res = 0.0d0 do 50 j=1,n 50 res = res + w(j)*v2(k2+j) do 60 j=1,n 60 v1(k1+j) = v1(k1+j) - res*v2(k2+j) c c... normalize set v1 c k1 = 0 do 100 i1=1,n1 call cntrc(s,v1(k1+1),w,n,critlow) ss = ddot(n,v1(k1+1),1,w,1) if (ss.lt.critlow) go to 200 tnorm = 1.0d0/dsqrt(ss) do 90 j=1,n 90 v1(k1+j) = v1(k1+j)*tnorm 100 k1 = k1 + n c return c 200 write(6,620) critlow,i1 620 format(/' vectorset linear dependent within ',e12.5, 1 ' detected at vector no ',i5) call error(22) c return end subroutine cntrc(s, c, v, n, eps) c c cntrc computes the inproducts of the vector c with the rows c of the triangular matrix s c the results are stored in the vector v c implicit real*8(a-h,o-z) dimension s(*),c(*),v(*) do 10 mu=1,n 10 v(mu) = 0 ix=0 do 50 nu=1,n cnu=c(nu) if (dabs(cnu) .lt. eps) goto 50 munu=ix do 20 mu=1,nu munu=munu+1 20 v(mu)=v(mu)+s(munu)*cnu if(nu.eq.n) go to 50 nu1=nu+1 munu=munu+nu do 30 mu=nu1,n v(mu)=v(mu)+s(munu)*cnu 30 munu=munu+mu 50 ix=ix+nu return end c function sumvec(a,b,n) c implicit real*8(a-h,o-z) c dimension a(*),b(*) c s = 0.0d0 c do 10 i=1,n c 10 s = s + a(i)*b(i) c sumvec = s c return c end subroutine lowdins(c,sao,s,cmsq,lambda,vecmt,n,m,epsss) c c subroutine lowdin performs a lowdin (symmetric) orthonormalization c of a given set of n vectors of dimension m c input are the vector-array c and the triangular primitive c overlap-matrix sao c the results are returned in the vector-input array c c arrays c c : n*m c sao : m*(m+1)/2 c s : n*(n+1)/2 c cmsq : n*n c lambda : n c vecmt : m c c implicit real*8(a-h,o-z) dimension c(*),sao(*),s(*),lambda(*),cmsq(*),vecmt(*) real*8 lambda c c form the overlap-matrix s c call fmos(s, sao, c, vecmt, n, m, epsss) k = n*(n+1)/2 c c diagonalize the s-matrix c call jacobi(s,cmsq,n,k,epsss) c c form s**(-.5) c c make a row of the negative squareroot of the eigenvalues c k = 0 do 70 i=1,n k = k + i 70 lambda(i) = 1/dsqrt(s(k)) c c perform the matrix-multiplication c k = 0 do 90 i=1,n kcmin = i - n do 80 j=1,n kcmin = kcmin + n 80 vecmt(j) = cmsq(kcmin)*lambda(j) do 90 j=1,i kl = j - n res = 0.0d0 do 85 ij =1,n kl = kl + n 85 res = res + vecmt(ij)*cmsq(kl) k = k + 1 90 s(k) = res c c bring the lowdin-transformation-matrix on rectangle-form c k = 0 l = 0 iend = n - 1 do 110 i=1,iend do 100 j=1,i k = k + 1 l = l + 1 100 cmsq(l) = s(k) jbeg = i + 1 kl = k + i do 110 j=jbeg,n l = l + 1 cmsq(l) = s(kl) 110 kl = kl + j do 120 i=1,n k = k + 1 l = l + 1 120 cmsq(l) = s(k) c c transformation of the vectors c call matmpl(c,cmsq,vecmt,n,m) c return end subroutine fmos(smo, sao, c, v, nmdim, ndim, epsss) c c compute the s-matrix on mo-basis for nmsim mo's (c) c of dimension ndim with ao-overlap-marix sao c implicit real*8(a-h,o-z) dimension smo(1),sao(1),c(1),v(1) c ks = 1 ki = 1 c do 50 i=1,nmdim call cntrc(sao, c(ki), v, ndim, epsss) kj = 1 do 40 j=1,i smo(ks) = ddot(ndim,c(kj),1,v,1) ks = ks + 1 40 kj = kj + ndim 50 ki = ki + ndim c return end subroutine jacobi(a,b,n,nnhalf,eps) implicit real*8(a-h,o-z) dimension a(nnhalf), b(n,n) real*8 ii,jj,ij integer contr nmin1 = n - 1 do 10 i=1,n do 10 j=1,n 10 b(i,j) = 0.d0 do 20 i=1,n 20 b(i,i) = 1.d0 if (n.le.1) return 100 u = 0 l = 2 do 110 i = 2,n do 120 j = 2,i u = u + a(l)**2 l = l+1 120 continue l = l+1 110 continue u = dsqrt(u/nnhalf) if (u.lt.eps) goto 300 i = 1 j = 1 kii = 1 kjj = 1 kij = 1 contr = 0 goto 200 202 if (j.eq.nmin1) goto 400 i = j+1 j = i kii = kjj + j kjj = kii kij = kii 200 if (i.eq.n) go to 202 kij = kij + i i = i + 1 kii = kii + i ij = a(kij) if (dabs(ij).lt.u) goto 200 contr = contr + 1 ii = a(kii) jj = a(kjj) vm = (jj-ii)/2.d0 vn = dsqrt(ij*ij+vm*vm) cos = (vn+dabs(vm))/(2.d0*vn) co = dsqrt(cos) si = ij/(dsign(2.d0,vm)*vn*co) sin = si*si imin = i-1 jmin = j-1 iplus = i + 1 jplus = j+1 ik = kii - imin kj = kjj - jmin if (j.eq.1) goto 211 do 210 kv = 1,jmin v = a(ik) w = a(kj) a(ik) = co*v - si*w a(kj) = si*v + co*w ik = ik+1 kj = kj+1 210 continue 211 akjj = ii * sin + jj * cos + 2 * si * co * ij a(kjj) = akjj a(kij) = 0 ik = ik + 1 kj = kj + j if(j.eq.imin) goto 213 do 212 kv = jplus , imin v = a(ik) w = a(kj) a(ik) = co*v - si*w a(kj) = si*v + co*w ik = ik+1 kj = kj+kv 212 continue 213 a(kii) = ii + jj - akjj ik = ik + i kj = kj + i if(i.eq.n) goto 215 do 214 kv = iplus , n v = a(ik) w = a(kj) a(ik) = co*v - si*w a(kj) = si*v + co*w ik = ik+kv kj = kj+kv 214 continue 215 do 216 kv = 1,n v = b(kv,i) w = b(kv,j) b(kv,i) = co*v - si*w b(kv,j) = si*v + co*w 216 continue goto 200 400 u = u / 6 if(contr.lt.4) u = u/2 if (u.lt.eps) go to 100 contr = 0 i = 1 j = 1 kii = 1 kjj = 1 kij = 1 goto 200 300 return end subroutine matmpl(v,t,vi,n,m) c c matmpl performs the transformation of a n*m vector-matrix v with c a n*n transformation matrix t, both lineairly stored. c results are returned in array v. c one intermidiate array vi of length n is required c implicit real*8(a-h,o-z) dimension v(*),t(*),vi(*) jend = n*m do 30 l=1,m k = 0 do 20 i=1,n res = 0.0d0 do 10 j=l,jend,m k = k + 1 10 res = res + t(k)*v(j) 20 vi(i) = res k = 0 do 30 j=l,jend,m k = k + 1 30 v(j) = vi(k) c return end subroutine getsh(sh,ndim,ichoice,isonel) c c read one electron integrals from dumpfile (atmol) c implicit real*8(a-h,o-z) common/sector/num3,iblk3 c common/disc/irep,iout,iin,iun,iblkk,ipos(16),name(16) c common /cntrl/ lprint,maxdim,ncore logical lprint c dimension sh(*) c*ap* for nav=1 : 72 // nav=2 : 77 c*ap* parameter (nints = 72, nav=1) parameter (nints = 77, nav=2) common/blkin/pp(4),y(6*nints),in(nints),nword,dumdum(6) c if (ichoice.lt.0.or.ichoice.gt.2) call error(9999) c naint =ndim*(ndim+1)/2 numblk=(naint-1)/nints+1 call secget(isonel,2,iblkd) call search(iblkd,num3) jst= 0 if(ichoice.eq.2) jst=nints*2 do 500 i=1,numblk call find(num3) call get(pp,nw) do 400 j=1,nword sh(in(j)) = y(jst + j) 400 continue 500 continue c return end subroutine tapnm(ifile,iblock,isect,init) c c... read a file/tape name/ block number and section c... init = .true. initialise dumpfile / .false. existing dumpfile c implicit real*8(a-h,o-z) logical init c common/disc/irep,iout,iin,iun,iblkk,ipos(16),name(16) common/discc/ name(16) character*4 name,nnam,ytest common/work/ jrec,jump common /cntrl/ lprint,maxdim,ncore logical lprint data mrewin,mblan/'rewind',' '/ c iblock = 1 c call inpa4(ytest) c do 10 i=1,16 if (ytest.eq.name(i)) go to 100 10 continue call error(13) c 100 ifile = i if (jrec.lt.jump) call inpi(iblock) if (init) call inidum(ifile,iblock) call secini(iblock,ifile) if (jrec.lt.jump) call inpi(isect) c return end subroutine inidum(num3,iblk3) c... start a new dumpfile implicit real*8(a-h,o-z) common /sector/ numd,iblkd,revise,rv,aiiii(204),maxblk,kblkla logical revise,rv c... **note** gamess sector has 2 revises for padding c revise = .false. call zero(aiiii,204) kblkla = 1 maxblk = 7777 numd = num3 iblkd = iblk3 c call revind c return end subroutine goback c... go back 1 variable on a card implicit real*8(a-h,o-z) common/work/jrec,jump c... jrec = jrec - 1 c... return end subroutine pass(ii,print) implicit real*8(a-h,o-z) logical print c c... read password and see if you recognise it c common /workc/ icard character*80 icard parameter (nww=31) character*8 iww(nww),itest data iww/'title','read','write','init','s','print', 1 'schmidt','lowdin','moperm','aoperm','combine','extra', 2 'summary','occup','eigen','crit','stop','finish', 3 'ndim','log','nmdim','addvec','symmetry', 4 'maxdim','check', 'locmp','convert', 5 'sdiag','tranform','screen','normalis'/ c 10 if (.not.print) write(6,3) 3 format(' *servec* ',$) call reada(itest) if (print) write(6,1) icard 1 format(/1x,a80,/) do 20 i=1,nww if (itest(1:4).eq.iww(i)(1:4)) go to 100 20 continue if (print) go to 50 write(6,1) icard write(6,2) 2 format(' the directive above is not recognised / please retry '/) call input go to 10 50 call error(10) c 100 ii = i c return end subroutine rdistr(istr,nn) c c... read a string numbers (including 'to') until end is encountered c implicit real*8(a-h,o-z) dimension istr(*) logical to character*8 itest,iend,ito data iend,ito/'end','to'/ c nn = 0 to = .false. call izero(istr,255) c 10 call reada(itest) if (itest.eq.iend) go to 100 if (to) go to 20 if (itest.ne.ito) go to 20 to = .true. go to 10 c 20 call goback call inpi(ne) if (to) go to 30 nn = nn + 1 istr(nn) = ne go to 10 c 30 nb = istr(nn) nn = nn - 1 do 40 ii=nb,ne nn = nn + 1 40 istr(nn) = ii c to = .false. go to 10 c 100 if (to) nn = 255 c return end subroutine pinver(ip,ipb,nd) c c... invert permutation ip => ipb c implicit real*8(a-h,o-z) dimension ip(nd),ipb(nd) c do 10 i=1,nd 10 ipb(ip(i)) = i c return end subroutine extr(vco,vcn,nold,new) c c... add unit matrix to set of old vectors / vcn may be vco c implicit real*8(a-h,o-z) dimension vco(nold,nold),vcn(new,new) c if (new.le.nold) return nold1 = nold + 1 c... copy vco to vcn do 20 i=1,nold do 10 j=1,nold vcn(nold1-j,nold1-i) = vco(nold1-j,nold1-i) 10 continue 20 continue c... add unit matrix do 40 i=nold1,new do 30 j=1,new vcn(i,j) = 0.0d0 vcn(j,i) = 0.0d0 30 continue 40 vcn(i,i) = 1.0d0 c return end subroutine prvc (c,n,m,occ,eps,ich,nc) implicit real*8(a-h,o-z) c c c routine for printing of column-vectors and eigenvalues. c c c : matrix of coefficients c n : dimension of c(n,m) and eps(n) c m : x columns to be printed c eps: vector of eigenvalues c occ: vector of occupation-numbers c n : dimension of c(n,n) and eps(n) c ich: 1 => occ / 0 => column-numbers / -1 => eig c nc : maximum number of columns on one row c c dimension c(n,m), eps(m) , occ(m) dimension sub(2) data sub /'-----','-----'/ ncc = nc if ( ncc.gt.10 ) ncc=10 nbl = (m-1)/ncc nlast = m-ncc*nbl nbl = nbl+1 nc1 = 0 do 60 ib=1,nbl if ( ib .eq. nbl ) ncc=nlast nc0 = nc1+1 nc1 = nc1+ncc if (ich.lt.0) write(6,10) ( eps(ic), ic=nc0,nc1 ) if (ich.gt.0) write(6,10) ( occ(ic), ic=nc0,nc1 ) if (ich.eq.0) write(6,11) ( ic, ic=nc0,nc1 ) 10 format(////8x,10f12.6) 11 format(////8x,10(4x,i4,4x)) write(6,20) ( sub, i=1,ncc ) 20 format(8x,10(2x,2a5)) write(6,30) 30 format(1h ) do 50 ia=1,n write(6,40) ia, ( c(ia,ic), ic=nc0,nc1 ) 40 format(3x,i3,2x,10f12.6) 50 continue 60 continue write(6,30) c return end subroutine symmetr(vc,eig,occ,ip,hulp,ndsub,ndim,nmdim,crit) c c... select mo's by 'symmetry'-relations c implicit real*8(a-h,o-z) dimension vc(ndim,nmdim),eig(nmdim),occ(nmdim),ip(ndim),hulp(ndim) dimension ipos(8,12),ich(8,12),nsbs(12),ndsub(nmdim) c common/cntrl/ lprint,maxdim,ncore logical lprint c call inpi(nao) call inpi(nirr) if (nao.gt.8.or.nao.le.0.or.nirr.gt.12.or.nirr.le.0)call error(15) c nsub = 0 nn = 0 10 nsub = nsub + 1 call inpi(ndsub(nsub)) if (ndsub(nsub).le.0) ndsub(nsub) = nmdim - nn nn = nn + ndsub(nsub) if (nn.lt.nmdim) go to 10 if (nn.gt.nmdim) call error(16) c do 30 irr=1,nirr call input do 20 iao=1,nao call inpi(ipos(iao,irr)) if (ipos(iao,irr).lt.0.or.ipos(iao,irr).gt.nmdim) 1 call error(17) call inpi(ich(iao,irr)) 20 continue 30 continue c if (.not.lprint) go to 50 write(6,602) nao,nirr,nsub do 40 irr=1,nirr 40 write(6,603) irr,(ipos(iao,irr),ich(iao,irr),iao=1,nao) 602 format(/' symmetry selection on',i3,' ao s for',i3, 1 ' representations in',i4,' subsets',//, 2 ' representation ao;ich etc ....') 603 format(1x,i7,4x,'/',8(i5,i5,3x,'/')) 50 continue c c... loop over subsets c kkv = 0 do 100 is=1,nsub nd = ndsub(is) call izero(ip,nd) kk = 0 do 80 irr=1,nirr do 60 kv=1,nd if (icomps(vc(1,kkv+kv),ipos(1,irr),ich(1,irr),nao,crit) 1 .eq.0) go to 60 kk = kk + 1 ip(kk) = kv 60 continue 80 nsbs(irr) = kk call coperm(ip,hulp,nd) c if (lprint) write(6,604) is,(ip(i)+kkv,i=1,nd) c604 format(' permutation for subset ',i4,/,(5x,16i4)) call change(ip,occ(kkv+1),eig(kkv+1),vc(1,kkv+1),1,ndim, 1 hulp,ndim,nd) do 70 i=2,nirr 70 nsbs(nirr-i+2) = nsbs(nirr-i+2) - nsbs(nirr-i+1) if (lprint) write(6,600) is,nd,(nsbs(i),i=1,nirr) 600 format(' subset ',i3,' (dim ',i3,' ) repr.dim : ',12i4) 100 kkv = kkv + nd c write(6,601) 601 format(1x) c return end integer function icomps(v,ipos,ich,nao,crit) c c... see if v is in representation described by ipso/ich c implicit real*8(a-h,o-z) dimension v(*),ipos(*),ich(*) c icomps = 0 ipar = 0 c do 100 i=1,nao if (ich(i)) 10,20,30 c... -1 10 if (ipar.eq.0) go to 15 if (ipar*v(ipos(i)).lt.-crit) go to 100 return 15 ipar = 1 if (v(ipos(i)).gt.crit) ipar = -1 go to 10 c... 0 20 if (dabs(v(ipos(i))).lt.crit) go to 100 return c... +1 30 if (ipar.eq.0) go to 35 if (ipar*v(ipos(i)).gt.crit) go to 100 return 35 ipar = 1 if (v(ipos(i)).lt.-crit) ipar = -1 go to 30 100 continue c icomps = 1 c return end subroutine checks(vc,s,smo,hulp,ndim,nmdim,crit) c c c... check orthogonality of vectors vc and print s-matrix (mo-basis) c... if not a unit-matrix within threshold crit c implicit real*8(a-h,o-z) dimension vc(*),s(*),smo(*),hulp(*) c data mprin/'print'/ c call fmos(smo,s,vc,hulp,nmdim,ndim,crit) c c prepare diagonal (substract 1.0 ) and check c k = 0 kkd = 0 ddev = 0.0d0 do 10 i=1,nmdim k = k + i smo(k) = smo(k) - 1.0d0 if (dabs(smo(k)).lt.crit.or.dabs(smo(k)).le.ddev) go to 10 ddev = dabs(smo(k)) kkd = i 10 continue c c check off-diagonal c kkoffx = 0 kkoffy = 0 offdev = 0.0d0 if (nmdim.le.1) go to 100 k = 0 do 20 i=2,nmdim k = k + 1 je = i - 1 do 20 j=1,je k = k + 1 if (dabs(smo(k)).lt.crit.or.dabs(smo(k)).le.offdev) go to 20 offdev = dabs(smo(k)) kkoffx = i kkoffy = j 20 continue c 100 continue c... print non-orth-warnings if (kkd.ne.0) write(6,1) kkd,ddev 1 format(' ********** check ********** ',/, 1 ' max. deviation from normality at',i4,' of ',e12.5) if (kkoffx.ne.0) write(6,2) kkoffx,kkoffy,offdev 2 format(' ********** check ********** ',/, 1 ' max. deviation from orthogonality at',i4,' /',i4, 2 ' of ',e12.5) c c print if there's a 'print' on the directive c call inpa(itest) if (itest.ne.mprin) return c write(6,3) nmdim,ndim 3 format(/' the (s-i) matrix for ',i4,' vectors (dim',i4,') is ') c call prtrs(smo,nmdim) c return end subroutine prtrs(s,n) c c triangle print for s-matrices c implicit real*8(a-h,o-z) dimension s(*) k = 1 do 10 i=1,n je = k + i - 1 write(6,1) i,(s(j),j=k,je) 10 k = k + i c 1 format(i5,(t6,10e12.5)) c return end subroutine symh(ihao,hmo,vc,occ,eig,isy,nmdim,ndim,crit,lprint) c c symmetry division using h-matrix blocking c ihao : h-matrix (ao-basis) + scratch space (ndim*(ndim+1)/2) c hmo : scratch (nmdim*(nmdim+1)/2) c vc : vector-set c occ/eig occupations + eigenvalues (must be permuted as well) c isy : scratch vector to store symmetries (ndim) c implicit real*8(a-h,o-z) dimension vc(ndim,nmdim),ihao(*),hmo(*),isy(ndim),occ(*),eig(*) dimension nd(20),ior(20) logical lprint c character*8 mend,mrenam,msets,ites data mend,mrenam,msets /'end','rename','sets'/ c... transform h-matrix to mo-basis c call fmos(hmo,ihao,vc,isy,nmdim,ndim,crit) c c... sort symmetries out c call symdiv(hmo,nmdim,nsym,isy,crit) c if (lprint) write(6,1) (isy(i),i=1,nmdim) 1 format(/' symmetry selection on h-matrix blocking / ', 1 ' global assignment :',/,(20i3)) if (lprint) write(6,7) nset = 1 nd(1) = nmdim c c... try to read more instructions c 5 call input call inpa(ites) c if (ites.eq.mend) go to 100 if (ites.eq.mrenam) go to 10 if (ites.eq.msets) go to 20 call goback go to 100 c c... rename c 10 nn = 0 11 nn = nn + 1 call inpi(ior(nn)) if (ior(nn).ne.0) go to 11 if (nn.ne.nsym+1) call error(24) c... change names of representations do 12 i=1,nmdim 12 isy(i) = ior(isy(i)) c if (lprint) write(6,2) (isy(i),i=1,nmdim) 2 format(' new global assignment :',/,(20i3)) if (lprint) write(6,7) c go to 5 c c... sets c 20 nset = 0 nn = 0 21 nset = nset + 1 if (nset.gt.20) call error(123) call inpi(nd(nset)) nn = nn + nd(nset) if (nd(nset).gt.0) go to 21 nd(nset) = nmdim - nn if (nd(nset).eq.0) nset = nset - 1 go to 5 c c... end or not recognized word c 100 if (lprint) write(6,3) nset,(nd(i),i=1,nset) 3 format(' division of orbital space into',i4,' sets / dimensions ', 1 10i4) kk = 0 ks = 0 do 39 i=1,nset nnd = nd(i) do 35 j=1,nsym ior(j) = 0 do 33 k=1,nnd if (isy(k+kk).ne.j) go to 33 ior(j) = ior(j) + 1 ks = ks + 1 ihao(ks) = kk + k 33 continue 35 continue if (lprint) write(6,6) i,(ior(j),j=1,nsym) 6 format(' subset',i3,' dim. irrep. repr. :',20i4) 39 kk = kk + nnd c if (lprint) write(6,4) (ihao(i),i=1,nmdim) 4 format(1x,/,' total symmetry permutation :',/,(5x,30i4)) if (lprint) write(6,7) 7 format(1x) c call change(ihao,occ,eig,vc,1,ndim,hmo,ndim,nmdim) c return end subroutine symdiv(hmat,ndim,nrep,isy,critsym) c c... find symmetry classification from 1-electron matrix c... results (symmetry-numbers) are returned in isy c implicit real*8(a-h,o-z) dimension hmat(*),isy(ndim) c logical ch ch(i,j) = dabs(hmat(j*(j-1)/2+i)).le.critsym c call izero(isy,ndim) c nrep = 0 c 10 continue do 100 i=1,ndim if (isy(i).ne.0) go to 30 nrep = nrep + 1 isy(i) = nrep 30 if (i.eq.ndim) return jb = i+1 do 50 j=jb,ndim if (ch(i,j)) go to 50 if (isy(j).le.0.or.isy(j).eq.isy(i)) go to 40 call comsym(isy(i),isy(j),isy,ndim,nrep) go to 50 40 isy(j) = isy(i) 50 continue 100 continue c c return end subroutine comsym(ir1,ir2,iro,norb,nrep) c c... symmetry ir1 and ir2 are identical really c... change the highest to the lowest one c implicit real*8(a-h,o-z) dimension iro(norb) c irf = max0(ir1,ir2) irt = min0(ir1,ir2) nrep = nrep - 1 c 10 do 20 i=1,norb 20 if (iro(i).eq.irf) iro(i) = irt c c... shift everything to get a consecutive list of symm.numbers again c if (irf.gt.nrep) return irt = irf irf = irf + 1 go to 10 c end subroutine form(vc,ndim,nmdim,occ,eig,iocc,ieig,ipr) c c... print (ipr=1) or read (ipr=-1) a complete vector file 'formatted' c... this option is intended for transport of vectors as text c implicit real*8(a-h,o-z) dimension vc(*),occ(*),eig(*) character*8 dumbo(3),mblan,mvecto,mocc,meig,dumb,ites data dumbo/'(',' ',')'/,dumb/'5f16.9'/,mblan/' '/ data mvecto,mocc,meig/'vectors','occ','eig'/ c if (ipr.lt.0) go to 1000 c... print / read format first call inpa(dumbo(2)) if (dumbo(2).eq.mblan) dumbo(2) = dumb c c... begin printing c write(6,1) 1 format(//' read format') write(6,2) ndim,nmdim,dumbo(2) 2 format(/' vectors ',2i10,5x,a8) nnnd = ndim*nmdim write(6,dumbo) (vc(i),i=1,nnnd) nn = nmdim if (iocc.le.0) nn = 0 write(6,3) nn 3 format(/' occ ',i10) if (nn.gt.0) write(6,dumbo) (occ(i),i=1,nmdim) nn = nmdim if (ieig.le.0) nn = 0 write(6,4) nn 4 format(/' eig ',i10) if (nn.gt.0) write(6,dumbo) (eig(i),i=1,nmdim) write(6,5) 5 format(/' end of format-print') c c... end of print c return c c... read c 1000 call reada(ites) if (ites.ne.mvecto) go to 1000 call inpi(ndim) call inpi(nmdim) call inpa(dumbo(2)) if (dumbo(2).eq.mblan) dumbo(2) = dumb nnnd = ndim*nmdim read(5,dumbo) (vc(i),i=1,nnnd) call reada(ites) if (ites.ne.mocc) call error(39) call inpi(nn) if (nn.gt.0) read(5,dumbo) (occ(i),i=1,nmdim) iocc = 0 if (nn.gt.0) iocc = 1 call reada(ites) if (ites.ne.meig) call error(39) call inpi(nn) if (nn.gt.0) read(5,dumbo) (eig(i),i=1,nmdim) ieig = 0 if (nn.gt.0) ieig = 1 write(6,6) 6 format(/' end of format-read') c return end subroutine locmp(u,sover,iloc,g,ntot,norb,nrmx,var) c implicit real*8(a-h,o-z) integer g(*),iloc(*) real*8 u(*),sover(*) character*8 ites,iend c c localization c v.magnasco,a.perico jcp 47,971 (1967) c c ntot = size of subset c norb = dimension of orbital basis c nrmx=maximum number of cycles to localize c var =accuracy of localization c c u orbitals (ntot*norb) c sover overlap matrix (norb*(norb+1)/2) c iloc scratch array (norb) c g scratch array (ntot*norb) c data simin/1.0d-7/,iend/'end'/ c nmx=ntot*norb do 10 i=1,nmx 10 g(i) = 0 nrot = 0 nmoc = 0 c write(6,601) 601 format(/' magnasco/perico localisation ') c 20 call reada(ites) if (ites.ne.iend) then call goback nmoc = nmoc + 1 max = 0 30 max = max + 1 call readii(iloc(max)) if (iloc(max).ne.0) go to 30 max = max - 1 write(6,602) nmoc,(iloc(j),j=1,max) 602 format(/,1x,i3,(2x,40i3)) idx=norb*(nmoc-1) do 40 j=1,max jdx=idx+iabs(iloc(j)) 40 g(jdx)=isign(1,iloc(j)) go to 20 end if c... print if it fits on the screen easily if (norb.le.10) then il = 0 do 50 j=1,nmoc is = il + 1 il = il + norb 50 write(6,603) (g(i), i = is,il) 603 format(1x,(/1x,40i3)) end if c c... start rotating the orbitals c 60 sd = 0.0d0 nrot=nrot+1 do 90 i=1,nmoc idx=norb*(i-1) iplus=i+1 if (iplus.gt.ntot) go to 100 do 90 j=iplus,ntot jdx=norb*(j-1) sa=0.0d0 sb=0.0d0 ndx=0 do 70 m=1,norb im=idx+m jm=jdx+m do 70 n=1,m ndx=ndx+1 in=idx+n jn=jdx+n if (j.le.nmoc) then fac=g(im)*g(in)-g(jm)*g(jn) else fac=g(im)*g(in) end if if (fac.ne.0.0d0) then fac=fac*sover(ndx) if (m.ne.n) fac=fac+fac sa=sa+fac*(u(im)*u(in)-u(jm)*u(jn)) sb=sb+fac*(u(im)*u(jn)+u(jm)*u(in)) end if 70 continue sc= dsqrt(sa*sa+sb*sb) if (sc.ne.0.0d0) then cosine=(sa+sc)/(sc+sc) if(cosine.lt.0.0d0) cosine=0.0d0 sine=1.0d0-cosine if (sine.ge.simin) then sine = dsign(dsqrt(sine),sb) else sine = sb/(sa+sa) sine = sine*(1-1.5d0*sine**2) end if cosine = dsqrt(cosine) do 80 m=1,norb im=idx+m jm=jdx+m sa=u(im) sb=u(jm) uim = sa*cosine + sb*sine ujm = -sa*sine + sb*cosine sd = sd + dabs(u(im)-uim) + dabs(u(jm)-ujm) u(im) = uim u(jm) = ujm 80 continue end if 90 continue 100 continue if (var.lt.sd.and.nrot.lt.nrmx) goto 60 c c... end of calc c write (6,604) nrot,sd 604 format(5h0****,' number of rotations :',i3,' convergence : ', 1 e12.5,5h ****) if (nrot.gt.nrmx) write(6,605) 605 format(/' ##### warning no convergence ##### '/) write (6,606) 606 format(/' ******* orbitals localised *****'/) c return end subroutine conver(ic,q,iq,maxdim) c*** convert dumpfile from gamess to atmol (ic=1) or vice versa (-1) implicit real*8 (a-h,p-w),integer (i-n),logical (o) implicit character *8 (z),character *1 (x) implicit character *4 (y) dimension q(*),iq(*) c*** q must be big enough to hold the vectors // maxdim include 'sizes' parameter (mxorb3=maxorb*3) c*** c*** junk/tran gamess/atmol format for type-3 section c*** c common/junkg1/zcom(19),ztitle(10) common/junkg1/com(19),title(10) common/junkg2/value(maxorb),pop(maxorb),escf, *nbasis,newbas,ncol,ivalue,ipop,nextra c.. nextra may disappear again common/junka/ coma(19),titla(10),valua(255),poa(255),nbasia(6) common/trang/ilifc(maxorb),ntran(maxorb),itran(mxorb3), * ctran(mxorb3),otran,odum common/trana/ilifa(256),ntraa(256),itraa(600),ctraa(600),otraa, * odum2 c*** common/discc/ yed(16) character*4 yed character*8 gatmol(2) data gatmol/'gamess','atmol'/ c*** n544 = 539 + (5-1)/nipw()+1 n1713 = (1113-1)/nipw() + 1 + 600 c call tapnm(numd,iblkd,isecd,.false.) if (ic.gt.0) then write(6,1) yed(numd),iblkd,gatmol(1),gatmol(2) else write(6,1) yed(numd),iblkd,gatmol(2),gatmol(1) end if 1 format(/' ** convert dumpfile on ',a4,' at block ',i5,' from ', 1 a8,' to ',a8) c*** c*** nav = # integers/real supposed to be 1 c*** c*** first copy ed3 to a scratch (fortran) file c*** open(99,status='scratch',form='unformatted',recl=10000) capp open(99,file='martyn',form='unformatted') do 10 i=1,204 call sectsg(i,iret,ipos,iclass,ilen) if (iret.gt.0) then if (iclass.eq.3.and.ic.gt.0) then c... for gamess look out for phony type3 sections call search(ipos,numd) call fget2(numd,q,nw) if (nw.ne.29) iclass = 3333 end if write(99) i,iclass,ilen call search(ipos,numd) do 5 j=1,ilen call fget2(numd,q,nw) write(99) nw,(q(k),k=1,nw) 5 continue end if 10 continue write(99) 0,0,0 c*** c*** now rewrite dumpfile in new format c*** rewind 99 call inidum(numd,iblkd) c 20 read(99) mpos,iclass,ilen if (iclass.eq.3.and.ilen.le.6.or.iclass.eq.3333) write(6,3) mpos 3 format(' section ',i3,' not converted ** no vectors ** ') if (iclass.eq.3.and.ilen.gt.6) then if (ic.gt.0) then c... c... gamess ==> atmol c... c... header-block (29 + maxorb*2 + 6 + 1 ==> 29 + 255*2 + 5 call read99(coma,29) call read99(value,maxorb*2+6/nipw()+1) c.. fill atmol-block call fmove(value,valua,255) call fmove(pop,poa,255) call fmove(nbasis,nbasia,6/nipw()) c.. write(6,2) mpos,nbasis,titla 2 format(' convert section ',i3,' dimension',i4,' title ',10a8) c... ctrans (maxorb*2 + mxorb3*2 + 1 => 256*2 + 600*2 +1) call read99(ilifc,2*maxorb/nipw()+ * mxorb3/nipw()+mxorb3+1) c.. fill atmol-type block call fmove (ilifc,ilifa,256/nipw()) call fmove (ntran,ntraa,256/nipw()) call fmove (itran,itraa,600/nipw()) call fmove (ctran,ctraa,600) otraa = otran c... read orbitals in q nbsq= nbasis*ncol if (nbsq.gt.maxdim*maxdim) call error(202) call read99(q,nbsq) c.... write type3 section to dumpfile ilen = 2 + 4 + (nbsq-1)/511 + 1 call secput(mpos,3,ilen,iblk) call wrt3(coma,n544,iblk,numd) call wrt3s(ilifa,n1713,numd) call wrt3s(q,nbsq,numd) else c*** c*** atmol => gamess c*** c... c... header-block c... call read99(coma,n544) c... load gamess block call fmove(coma,com,29) call fmove(valua,value,255) call fmove(poa,pop,255) call fmove(nbasia,nbasis,6/nipw()) escf = 0.0d0 nextra = 0 c.. write(6,2) mpos,nbasis,titla c... ctrans call read99(ilifa,n1713) call fmove (ilifa,ilifc,256/nipw()) call fmove (ntraa,ntran,256/nipw()) call fmove (itraa,itran,600/nipw()) call fmove (ctraa,ctran,600) otran = otraa c... orbitals nbsq= nbasis*ncol if (nbsq.gt.maxdim*maxdim) call error(202) call read99(q,nbsq) c.... write type3 section to ed3 ilen = 1 + (maxorb*2+6/nipw()-1)/511 + 1 + 1 (maxorb*2/nipw()+mxorb3/nipw()+mxorb3+1-1)/511 2 + 1 + (nbsq-1)/511 + 1 call secput(mpos,3,ilen,iblk) call wrt3(com,29,iblk,numd) call wrt3s(value,maxorb*2+6/nipw()+1,numd) call wrt3s(ilifc,maxorb*2/nipw()+mxorb3/nipw() 1 +mxorb3+1,numd) call wrt3s(q,nbsq,numd) end if else if (mpos.eq.0) then c*** end of tape99=martyn call clear call revind close(99) return c*ap* c*ap* temporary fix also change type 2 c*ap* else if (iclass.eq.2) then call secput(mpos,iclass,ilen,iblock) call search(iblock,numd) write(6,4) mpos 4 format(' convert 1-integral section ',i5) do 25 j=1,ilen read(99) nn,(q(i),i=1,nn) if (ic.gt.0) then c.. gamess => atmol iq(1010) = iq(1020) call put(q,505,numd) else c.. atmol => gamess (**not recommended**) iq(1020) = iq(1010) iq(1021) = 1 iq(1022) = 1 call put(q,511,numd) end if 25 continue c else c... normal section iclass=3333 translates to iclass=3 if (iclass.eq.3333.and.ic.le.0) iclass=3 call secput(mpos,iclass,ilen,iblock) call search(iblock,numd) do 30 j=1,ilen read(99) nn,(q(i),i=1,nn) call put(q,nn,numd) 30 continue end if c*** go to 20 c*** end subroutine read99(q,nt) c*** read exactly nt numbers of tape99 (check) implicit real*8(a-h,o-z) dimension q(nt) kk = 0 10 read(99) nn,(q(kk+i),i=1,nn) kk = kk + nn if (kk.lt.nt) go to 10 if (kk.ne.nt) then write(6,1) nt,kk,nn 1 format(' tried to read ',i5,' read ',i5,' last record ',i5) call error(203) end if c*** return end subroutine sectsg(mpos,iretrn,ipos,iclass,ilen) c*** secget that tests if a section is there c*** and returns the vital statistics if it's there c*** *** note difference in cray/atmol format for //sector// c*** so for cray atmol-format **cray** implicit real*8 (a-h,p-w),integer (i-n),logical (o) implicit character *8 (z),character *1 (x) implicit character *4 (y) common/sector/num3,iblk3,orevis(2), 1 apos(204),maxb,kblkla call upack3(apos(mpos),ipos,iclass,ilen) iretrn=0 if (ipos.gt.0) then iretrn = 1 ipos = ipos + iblk3 end if c*** return end subroutine jacobs(s,vc,eig,iscr,n,crit,isort) c... c... diagonalise matrix and yield sorted vectors/values c... isort : 0,1 : nop / 2 increasing / 3 decreasing c... implicit real*8(a-h,o-z) dimension s(*),vc(*),eig(*),iscr(*) c... c... diagonalise c... nn = n*(n+1)/2 call jacobi(s,vc,n,nn,crit) c... c... gather eigenvalues in s / eig c... nn = 0 do 10 i=1,n nn = nn + i eig(i) = s(nn) 10 continue c... c... sort ( s = scratch) c... call sort(eig,iscr,s,n,isort) c... c... permute vc/eig c... call change(iscr,eig,eig,vc,1,n,s,n,n) c return end subroutine sort(val,ip,hulp,n,isort) c c... determine permutation to sort val c... isort = 2 : increasing / 3 decreasing / else unit c implicit real*8(a-h,o-z) logical hulp dimension val(n),ip(n),hulp(n) c do 10 i=1,n 10 hulp(i) = .true. c do 30 ii=1,n do 14 i=1,n 14 if (hulp(i)) go to 15 15 ipi = i do 20 i=1,n if (hulp(i).and.( 2 (isort.eq.2.and.val(i).lt.val(ipi)) .or. 3 (isort.eq.3.and.val(i).gt.val(ipi)) ) ) then ipi = i end if 20 continue ip(ii) = ipi hulp(ipi) = .false. 30 continue c return end