PROGRAM QGRAD c c Takes S-matrix and cartesian gradient, c calculates normal mode gradient c IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NAT0=900,npar0=50) parameter (nsy0=13) DIMENSION x(3*NAT0),gx(3*NAT0),w(3*NAT0),io(3*NAT0), 1am(3*NAT0),gn(3*NAT0),dq(3*NAT0),ity(NAT0) ,gxold(3*NAT0), 2sx(3*NAT0),tv(3*NAT0),ZM(3*NAT0),C(3,NAT0), 3sq(3*NAT0),gnold(3*NAT0),EDIAG(3*NAT0),bijsame(npar0,2), 4xq(NAT0*3),cq(NAT0),pot(NAT0),potc(NAT0),ijsame(npar0,2) character*3 pg(nsy0) common/cx/s(3*NAT0,3*NAT0),fx(3*NAT0,3*NAT0),A(3*NAT0,3*NAT0) common/cx1/TEM(3*NAT0,3*NAT0),fq(3*NAT0,3*NAT0), 2wt(3*NAT0),st(3*NAT0,3*NAT0),qt(3*NAT0),gt(3*NAT0) common/fixpar/apar(npar0),bpar(npar0),gg(3*NAT0),p(4,3), 1ipar(npar0,4),kpar(npar0),pp(4,3,4,3),fxpar(3*NAT0,3*NAT0), 1p1(4,3),p2(4,3),pp1(4,3,4,3),pp2(4,3,4,3) logical*4 lupdate,lff,lm c call initoutput call loadg(NAT0,gx,nat,io,glim,e,wmin,wmax,lm) call loadpar(npar,ipar,kpar,npar0,bpar,apar,nsame,ijsame, 1bijsame) call loadx(NAT0,x,nat,ity,nq,xq,cq) call projectgr(gx,nat,x) call setp(x,nat,nq,xq,cq,pot) call findsymm(x,NAT0,nat,ne,pg,tol,ity) if(npar.gt.0)call fixeg(npar,e,ipar,kpar,npar0,bpar,apar,NAT0,x, 1gg,nat,p,pp,ijsame,bijsame,nsame,p1,p2,pp1,pp2) if(nq.gt.0)call fixegc(e,x,gg,nat,pot,potc,xq,cq,nq) call checkenergy(e,sx,x,ity,nat,NAT0,cq,nq,pot,xq) call loadhesopts(lupdate,lff) call loadm(NAT0,am,nat) call wx(NAT0,x,nat,ity) if(lff)then call loadff(NAT0,fx,nat) if(lupdate)call updateff(NAT0,fx,gx,nat,gxold,sx,tv) if(npar.gt.0)call fixff(npar,fx,ipar,kpar,npar0,gg,gx, 1bpar,apar,NAT0,x,nat,p,fxpar,pp,nsame,ijsame,bijsame, 2p1,p2,pp1,pp2) if(nq.gt.0)call fixffc(nat, 1fxpar,NAT0,x,gg,gx,fx,potc,pot,xq,cq,nq) call symmetrizef(NAT0,x,fx,nat,ne,pg,ity) call makes(NAT0,fx,s,nat,w,nm,am,A,TEM,ity,x,fq,ZM,C,EDIAG) else call loads(NAT0,s,nat,w,nm) endif call writeg(NAT0,gx,nat,npar,gg) call dimensionlesss(NAT0,s,am,nat,nm) call symmetrizeg(NAT0,x,gx,nat,ne,pg,ity) call normalg(NAT0,gn,gx,s,nat,nm,am) if(.not.lff.and.lupdate) 1call updatefq2(NAT0,fq,gn,gnold,sq,tv,nm,w) call writeng(NAT0,gn,nm) call report(NAT0,gn,w,io,nm,dq,e,glim,lupdate,lff,fq,lm, 1wmin,wmax,fx,wt,st,EDIAG,qt,gt,x,nat,ity,ne,pg,s) call newsx(NAT0,s,dq,nm,nat,io,am,sx) call symmetrizex(NAT0,sx,x,nat,ne,pg,ity) call savestep(NAT0,nat,sx) call newx(NAT0,x,nat,sx) call ginp(NAT0,x,nat,ity,cq,nq,pot,xq) call bye('QGRAD finished OK') stop end c ====================================================================== subroutine initoutput character*80 ts ist=0 open(3,file='QGRAD.OUT') 2 read(3,3333,end=200,err=200)ts 3333 format(a80) if(ts(15:24).eq.'- Step -')ist=ist+1 goto 2 200 backspace 3 write(3,3000)ist 3000 format(/,' ------------ Step --------------',/, 1 ' ------------ ',i6,' --------------',/) maxstep=-1 open(10,file='Q.OPT') 11 read(10,3333,end=999,err=999)ts if(ts(1:7).eq.'MAXSTEP')read(10,*)maxstep goto 11 999 close(10) if(maxstep.ge.0.and.ist.ge.maxstep)call bye('MAXSTEP overflow.') return end c ====================================================================== subroutine checkenergy(e,sx,x,ity,nat,NAT0,cq,nq,pot,xq) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) dimension sx(3*NAT0),x(3*NAT0),ity(NAT0),cq(*),pot(*),xq(*) logical*4 lex,lcct character*80 ts c c dll=0.05d0 dll=0.001d0 lcct=.false. open(10,file='Q.OPT') 11 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:9).eq.'del e max')read(10,*)dll if(ts(1:4).eq.'lcct')read(10,*)lcct goto 11 999 close(10) c inquire(file='E.SCR',exist=lex) if(lex)then open(10,file='E.SCR',form='unformatted') read(10)eold close(10) de=e-eold write(3,3009)e,eold,de,dll 3009 format(' Current energy ',f15.8,/, 1 ' Old energy ',f15.8,/, 2 ' -------------- ',/, 3 ' Difference ',f15.8,' tol ',f12.6/) c if((e.gt.eold.and.npar.eq.0.and.nq.eq.0).or.de.gt.dll)then if(lcct)then write(3,3000) 3000 format(' lcct = TRUE, no energy control') else if(de.gt.dll)then call msg(' Step will be reduced.') call loadstep(NAT0,nat,sx) call msg('Old step loaded') do 1 i=1,3*nat 1 sx(i)=sx(i)/sqrt(2.0d0) call msg('Step reduced') call savestep(NAT0,nat,sx) call msg('Step saved') call loadx(NAT0,x,nat,ity,nq,xq,cq) call msg('Coordinates loaded') bohr=0.529177249d0 do 2 i=1,3*nat 2 x(i)=x(i)+sx(i)*(1.0d0-sqrt(2.0d0))*bohr call msg('New coordinates produced') call ginp(NAT0,x,nat,ity,cq,nq,pot,xq) call msg('New input produced') call bye('Qgrad finished') endif endif endif open(10,file='E.SCR',form='unformatted') write(10)e close(10) return end c ====================================================================== subroutine loadstep(NAT0,nat,sx) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION sx(3*NAT0) c open(10,file='SX.SCR',form='unformatted') read(10)(sx(i),i=1,3*nat) close(10) return end c ====================================================================== subroutine savestep(NAT0,nat,sx) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION sx(3*NAT0) c open(10,file='SX.SCR',form='unformatted') write(10)(sx(i),i=1,3*nat) close(10) open(10,file='AX.SCR') do 1 ia=1,nat 1 write(10,100)(sx(3*(ia-1)+j),j=1,3) 100 format(3f12.6) close(10) call msg('Cartesian step written into SX.SCR.') return end c ====================================================================== subroutine updatefq2(NAT0,fq,gn,gnold,sq,tv,nm,w) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION fq(3*NAT0,3*NAT0),gn(3*NAT0),gnold(3*NAT0), 1sq(3*NAT0),tv(3*NAT0),w(3*NAT0) logical*4 lex c call msg('Internal update called') inquire(file='FN.SCR',exist=lex) if(lex)then open(10,file='FN.SCR',form='unformatted') read(10)((fq(i,j),i=1,nm),j=1,nm) close(10) call msg('Internal Hessian^-1 loaded from FN.SCR.') open(10,file='QGRADIENT.SCR',form='unformatted') read(10)(gnold(i),i=1,nm) close(10) call msg('Last gradient loaded from QGRADIENT.SCR.') open(10,file='SQ.SCR',form='unformatted') read(10)(sq(i),i=1,nm) close(10) call msg('Last normal step loaded from SQ.SCR.') c c BFGS update of the inverse Hessian do 2 i=1,nm sum=0.0d0 do 3 j=1,nm 3 sum=sum+fq(i,j)*sq(j) 2 tv(i)=sum c shs=0.0d0 sf=0.0d0 do 4 i=1,nm sf = sf+sq(i)*(gn(i)-gnold(i)) 4 shs=shs+sq(i)*tv(i) c ems=1.0d-10 do 1 i=1,nm fi=gn(i)-gnold(i) ti=tv(i) do 1 j=1,nm fj=gn(j)-gnold(j) tj=tv(j) c plus here because local Q-coordinate is different from global if(abs( sf).gt.ems)fq(i,j)=fq(i,j)+fi*fj/sf 1 if(abs(shs).gt.ems)fq(i,j)=fq(i,j)-ti*tj/shs c call msg('Internal Hessian updated.') else c SI constants: amu=1.6605402D-27 a0=0.529177249d-10 hartree=4.3597482d-18 clight=2.99792458d8 pi=3.141592653589d0 c constant that transfers cm-1 to au const=0.01d0/(2.0d0*pi*clight)*sqrt(hartree/(a0**2*amu)) do 5 i=1,nm do 5 j=1,nm 5 fq(i,j)=0.0d0 do 6 i=1,nm fqii=(w(i)/const)**2 if(w(i).lt.0.0d0)fqii=-fqii 6 fq(i,i)=fqii call msg('Internal Hessian initiated.') endif open(10,file='FN.SCR',form='unformatted') write(10)((fq(i,j),i=1,nm),j=1,nm) close(10) call msg('Internal Hessian written to FN.SCR.') return end c ====================================================================== subroutine updateff(NAT0,fx,gx,nat,gxold,sx,tv) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION fx(3*NAT0,3*NAT0),gx(3*NAT0),gxold(3*NAT0), 1sx(3*NAT0),tv(3*NAT0) logical*4 lex c inquire(file='GRADIENT.SCR',exist=lex) if(lex)then open(10,file='GRADIENT.SCR',form='unformatted') read(10)(gxold(i),i=1,3*nat) close(10) call loadstep(NAT0,nat,sx) c c BFGS update sf=0.0d0 do 2 i=1,3*nat sum=0.0d0 do 3 j=1,3*nat 3 sum=sum+fx(i,j)*sx(j) tv(i)=sum 2 sf=sf+sx(i)*(gx(i)-gxold(i)) shs=0.0d0 do 4 i=1,3*nat 4 shs=shs+sx(i)*tv(i) ems=1.0d-15 do 1 i=1,3*nat0 dfi=gx(i)-gxold(i) do 1 j=1,3*nat0 dfj=gx(j)-gxold(j) if(abs(sf).gt.ems)fx(i,j)=fx(i,j)-dfi*dfj/sf 1 if(abs(shs).gt.ems)fx(i,j)=fx(i,j)-tv(i)*tv(j)/shs if(abs(sf).gt.ems.or.abs(shs).gt.ems)then call msg('Force constants updated.') call WRITEFF(NAT0,nat,fx,'FILE.FC') else call msg('Zero grad/step; no update') endif else call msg('Force constants not updated - first point.') endif return end c ====================================================================== subroutine loadff(NAT0,fx,nat) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION fx(3*NAT0,3*NAT0) character*80 filename,ts c filename='FILE.FC' open(10,file='Q.OPT') 11 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:11).eq.'force field')read(10,1000)filename goto 11 999 close(10) c open(20,file=filename) N=3*nat N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(fx(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) c to comply with new4 transfer ff in H/bohr^2 to mdyn/A^2 or something DO 3 I=1,N DO 3 J=I,N 3 fx(J,I)=fx(J,I) DO 31 I=1,N DO 31 J=I+1,N 31 fx(I,J)=fx(J,I) close(20) call msg('Cartesian FF read in from '//filename//'.') RETURN END c ====================================================================== subroutine loadhesopts(lupdate,lff) logical*4 lupdate,lff character*80 ts c lupdate=.true. lff=.true. open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:14).eq.'hessian update')read(10,*)lupdate if(ts(1:16).eq.'cartesian update')read(10,*)lff goto 1 999 close(10) write(3,*)' Lupdate: ',lupdate write(3,*)' Lff : ',lff c return end c ====================================================================== function notemptyline(ts) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) LOGICAL*4 notemptyline CHARACTER*(*) ts c l=len(ts) notemptyline=.false. do 1 i=1,l 1 if(ts(i:i).ne.' ')notemptyline=.true. return end c ====================================================================== subroutine chkchek(ts,lchk) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) CHARACTER*(*) ts CHARACTER*16 s c l=len(ts) do 1 i=1,l-15 do 2 ii=1,16 ia=ichar(ts(i+ii-1:i+ii-1)) if(ia.ge.65.and.ia.le.90)ia=ia+32 2 s(ii:ii)=char(ia) 1 if(s.eq.'guess=checkpoint')lchk=1 return end c ====================================================================== subroutine chkoniom(ts,loniom) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) CHARACTER*(*) ts CHARACTER*16 s c l=len(ts) do 1 i=1,l-4 do 2 ii=1,5 ia=ichar(ts(i+ii-1:i+ii-1)) if(ia.ge.65.and.ia.le.90)ia=ia+32 2 s(ii:ii)=char(ia) 1 if(s.eq.'oniom')loniom=1 return end c ====================================================================== subroutine ginp(NAT0,x,nat,ity,cq,nq,pot,xq) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MENDELEV=89) LOGICAL*4 lnewcharges,ldalton,notemptyline,lturbo,lgauss, 1ltinker,lcct CHARACTER*80 filename,ts CHARACTER*5 s5 DIMENSION x(3*NAT0),ity(3*NAT0),cq(*),pot(*),xq(*),ib(7) CHARACTER*2 sy(MENDELEV) data sy/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne', 3'Na','Mg','Al','Si',' P',' S','Cl','Ar', 4' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn', 4 'Ga','Ge','As','Se','Br','Kr', 5'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd', 5 'In','Sn','Sb','Te',' I','Xe', 6'Cs','Ba','La', 6 'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho', 6 'Er','Tm','Yb','Lu', 6'Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg', 6 'Tl','Pb','Bi','Po','At','Rn', 7'Fr','Ra','Ac'/ c filename='G98.INP' lnewcharges=.false. ldalton=.false. lcct=.false. ltinker=.false. lturbo=.false. lgauss=.true. open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:14).eq.'gaussian input')read(10,1000)filename if(ts(1:11).eq.'new charges')read(10,*)lnewcharges if(ts(1:4).eq.'dalt')read(10,*)ldalton if(ts(1:4).eq.'lcct')read(10,*)lcct if(ts(1:4).eq.'tink')read(10,*)ltinker if(ts(1:4).eq.'turb')read(10,*)lturbo goto 1 999 close(10) c if(lcct.or.ldalton.or.lturbo.or.ltinker)lgauss=.false. if(ltinker)filename='tinker.xyz' if(ldalton)filename='d.mol' if(lcct)filename='FILE.X' if(lturbo)filename='coord' open(10,file=filename) open(11,file='INP.NEW') c ccccccccccc cct: cccccccccccccccccccccccccccccccccccc if(lcct)then read(10,1000)ts write(11,1000)ts read(10,*)nat write(11,*)nat do 6 ia=1,nat read(10,*)qt,qt,qt,qt,(ib(ii),ii=1,7),qt 6 write(11,1200)ity(ia),(x(3*(ia-1)+ix),ix=1,3),(ib(ii),ii=1,7),qt 1200 format(i5,3f12.6,7i5,f8.4) endif c cccccccccccccccccccccccccccccccccccccccccccccccccccccccc C DDDDDDDDDDD Dalton: DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD if(ldalton)then read(10,1000,end=899,err=899)ts write(11,1000)ts iab=1 if(ts(1:9).eq.'ATOMBASIS')iab=0 if(iab.eq.1)then read(10,1000,end=899,err=899)ts write(11,1000)ts endif do 901 i=1,3 read(10,1000,end=899,err=899)ts 901 write(11,1000)ts ia=0 902 read(10,1000,end=899,err=899)ts idp=0 do 903 i=1,80 903 if(ts(i:i).eq.'.')idp=idp+1 if(idp.eq.3)then ia=ia+1 write(11,7701)ts(1:3),(x(3*(ia-1)+ix)/0.529177d0,ix=1,3) 7701 format(a3,3f20.15) if(ia.lt.nat)goto 902 else write(11,1000)ts goto 902 endif endif C DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD C IIIIIIIIII Tinker: IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII if(ltinker)then read(10,1000,end=899,err=899)ts write(11,1000)ts do 905 ia=1,nat read(10,1000,end=899,err=899)ts 905 write(11,7702)ts(1:11),(x(3*(ia-1)+ix),ix=1,3), 1 ts(48:80) 7702 format(a11,3f12.6,a33) endif C IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII C TTTTTTTTTTTT Turbomole: TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT if(lturbo)then write(11,1001) read(10,*) 1001 format('$coord') do 9037 ia=1,nat read(10,7741)s5 7741 format(66x,a5) 9037 write(11,7731)(x(3*(ia-1)+ix)/0.529177d0,ix=1,3),s5 7731 format(f20.14,2f22.14,2x,a5) write(11,1002) read(10,*) read(10,*) 1002 format('$user-defined bonds',/,'$end') endif C TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT if(lgauss)then c GGGGGGGGGG Gaussian : GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG c c rewrite header lchk=0 loniom=0 2 read(10,1000,end=899,err=899)ts call chkchek(ts,lchk) call chkoniom(ts,loniom) if(notemptyline(ts))then write(11,1000)ts goto 2 else if(lchk.eq.0)write(11,1111) 1111 format('guess=checkpoint') write(11,1000)ts endif c c name 21 read(10,1000,end=899,err=899)ts write(11,1000)ts if(notemptyline(ts))goto 21 c c charge and multiplicity read(10,1000,end=899,err=899)ts write(11,1000)ts c c Cartesian coordinates - replace by new guess do 3 ia=1,nat if(loniom.eq.1)then read(10,1000)ts write(11,7700)ts(1:19),(x(3*(ia-1)+ix),ix=1,3),ts(56:80) 7700 format(a19,3f12.6,a25) else read(10,1000)ts write(11,1100)sy(ity(ia)),(x(3*(ia-1)+ix),ix=1,3) 1100 format(a2,2x,3f15.8) endif 3 continue c c fix point charges to keep constant potentials at nuclei: if(nq.gt.0)then read(10,1000,end=899,err=899)ts write(11,1000)ts do 4 i=1,nq 4 read(10,*)(xq(3*(i-1)+ix),ix=1,3),cq(i) if(lnewcharges)call genchar(nat,nq,pot,x,xq,cq) do 5 i=1,nq 5 write(11,119)(xq(3*(i-1)+ix),ix=1,3),cq(i) 119 format(4f12.6) endif c c optional parameters until end of file: 31 read(10,1000,end=899,err=899)ts write(11,1000)ts goto 31 endif c GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG c 899 close(10) close(11) call msg('New input written in INP.NEW') return end c ====================================================================== subroutine findsymm(x,NAT0,nat,ne,pg,tol,ity) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) parameter (nsy0=13) CHARACTER*3 pg(nsy0) CHARACTER*80 ts,pp DIMENSION x(3*NAT0),ity(NAT0) DIMENSION a2(nsy0,3,3) CHARACTER*3 sl(nsy0) c call filla(nsy0,a2,sl) tol=1.0d-2 pp(1:3)='all' open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:8).eq.'SYMMETRY')read(10,1000)pp if(ts(1:5).eq.'S-TOL')read(10,*)tol if(ts(1:5).eq.'NOSYM')then close(10) return endif goto 5 999 close(10) c c try all symmetries ne=0 do 3 is=1,nsy0 if(pp(1:3).eq.'all'.or.pp(1:3).eq.sl(is))then iy=1 do 1 ia=1,nat ic=0 do 2 ja=1,nat dist=sqrt((ax(is,x,a2,nsy0,ia,1)-x(3*(ja-1)+1))**2+ 1 (ax(is,x,a2,nsy0,ia,2)-x(3*(ja-1)+2))**2+ 2 (ax(is,x,a2,nsy0,ia,3)-x(3*(ja-1)+3))**2) c if(ia.eq.1.and.ja.eq.3)write(6,*)sl(is),dist 2 if(dist.lt.tol.and.ity(ia).eq.ity(ja))ic=ic+1 1 if(ic.eq.0)iy=0 if(iy.eq.1)then ne=ne+1 pg(ne)=sl(is) endif endif 3 continue c write(3,3000)ne,(pg(i),i=1,ne) 3000 format(i2,' symmetry elements: ',10(a3,1x)) c return end c ====================================================================== function symviol(i,s,x,nat,NAT0,ity,ne,pg) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) LOGICAL*4 symviol CHARACTER*80 ts parameter (nsy0=13) CHARACTER*3 pg(nsy0) DIMENSION x(3*NAT0),s(3*NAT0,3*NAT0),ity(NAT0) DIMENSION a2(nsy0,3,3) CHARACTER*3 sl(nsy0) c call filla(nsy0,a2,sl) c symviol=.false. c tol=1.0d-2 open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.'S-TOL')read(10,*)tol if(ts(1:5).eq.'NOSYM')then close(10) return endif goto 5 999 close(10) c c determine, if the mode is totally symmetric do 1 ie=1,ne do 2 is=1,nsy0 if(sl(is).eq.pg(ie))then do 3 ia=1,nat do 4 ja=1,nat dist=sqrt((ax(is,x,a2,nsy0,ia,1)-x(3*(ja-1)+1))**2+ 1 (ax(is,x,a2,nsy0,ia,2)-x(3*(ja-1)+2))**2+ 2 (ax(is,x,a2,nsy0,ia,3)-x(3*(ja-1)+3))**2) if(dist.lt.tol.and.ity(ia).eq.ity(ja))then diss=sqrt((a2(is,1,1)*s(3*(ia-1)+1,i)+ 1 a2(is,1,2)*s(3*(ia-1)+2,i)+ 1 a2(is,1,3)*s(3*(ia-1)+3,i)-s(3*(ja-1)+1,i))**2+ 1 (a2(is,2,1)*s(3*(ia-1)+1,i)+ 1 a2(is,2,2)*s(3*(ia-1)+2,i)+ 1 a2(is,2,3)*s(3*(ia-1)+3,i)-s(3*(ja-1)+2,i))**2+ 1 (a2(is,3,1)*s(3*(ia-1)+1,i)+ 1 a2(is,3,2)*s(3*(ia-1)+2,i)+ 1 a2(is,3,3)*s(3*(ia-1)+3,i)-s(3*(ja-1)+3,i))**2) if(diss.gt.tol)then symviol=.true. goto 997 endif endif 4 continue 3 continue endif 2 continue 1 continue 997 return end c ====================================================================== subroutine symmetrizex(NAT0,sx,x,nat,ne,pg,ity) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION x(3*NAT0),sx(3*NAT0) character*80 ts parameter (nsy0=13) character*3 pg(nsy0) DIMENSION ity(NAT0) DIMENSION a2(nsy0,3,3) CHARACTER*3 sl(nsy0) c call filla(nsy0,a2,sl) tol=1.0d-2 open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.'S-TOL')read(10,*)tol if(ts(1:5).eq.'NOSYM')then close(10) return endif goto 5 999 close(10) c do 10 ie=1,ne do 11 is=1,nsy0 if(pg(ie).eq.sl(is))then do 1 ia=1,nat do 1 ja=ia,nat dist=sqrt((ax(is,x,a2,nsy0,ia,1)-x(3*(ja-1)+1))**2+ 1 (ax(is,x,a2,nsy0,ia,2)-x(3*(ja-1)+2))**2+ 2 (ax(is,x,a2,nsy0,ia,3)-x(3*(ja-1)+3))**2) if(dist.lt.tol.and.ity(ia).eq.ity(ja))then sx(3*(ja-1)+1)=ax(is,sx,a2,nsy0,ia,1) sx(3*(ja-1)+2)=ax(is,sx,a2,nsy0,ia,2) sx(3*(ja-1)+3)=ax(is,sx,a2,nsy0,ia,3) endif 1 continue endif 11 continue 10 continue c call msg('The step has been symmetrized.') return end c ====================================================================== subroutine symmetrizef(NAT0,x,fx,nat,ne,pg,ity) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION fx(3*NAT0,3*NAT0),x(3*NAT0) character*80 ts parameter (nsy0=13) character*3 pg(nsy0) DIMENSION ity(NAT0) DIMENSION a2(nsy0,3,3) CHARACTER*3 sl(nsy0) c call filla(nsy0,a2,sl) c tol=1.0d-2 open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.'S-TOL')read(10,*)tol if(ts(1:5).eq.'NOSYM')then close(10) return endif goto 5 999 close(10) c tol2=tol**2 do 10 ie=1,ne do 11 iy=1,nsy0 if(pg(ie).eq.sl(iy))then do 1 ia=1,nat ii=3*(ia-1) xi=ax(iy,x,a2,nsy0,ia,1) yi=ax(iy,x,a2,nsy0,ia,2) zi=ax(iy,x,a2,nsy0,ia,3) iat=ity(ia) do 1 iap=ia,nat iip=3*(iap-1) xip=ax(iy,x,a2,nsy0,iap,1) yip=ax(iy,x,a2,nsy0,iap,2) zip=ax(iy,x,a2,nsy0,iap,3) iatp=ity(iap) do 1 ja=1,nat jj=3*(ja-1) d2 =(xi -x(jj +1))**2+(yi -x(jj +2))**2+(zi -x(jj +3))**2 if(d2 .lt.tol2.and.iat .eq.ity(ja ))then do 112 jap=1,nat jjp=3*(jap-1) d2p=(xip-x(jjp+1))**2+(yip-x(jjp+2))**2+(zip-x(jjp+3))**2 if(d2p.lt.tol2.and.iatp.eq.ity(jap))then do 111 ix=1,3 ist=1 if(ia.eq.iap)ist=ix do 111 ixp=ist,3 sss=0.0d0 do 2 iix=1,3 do 2 jjx=1,3 2 sss=sss+a2(iy,ix,iix)*fx(ii+iix,iip+jjx)*a2(iy,ixp,jjx) c write(6,600)jj+ix,jjp+ixp,sss,fx(jj+ix,jjp+ixp) c00 format(3(2i3,2f8.3)) fx(jj+ix,jjp+ixp)=sss 111 fx(jjp+ixp,jj+ix)=sss endif 112 continue endif 1 continue endif 11 continue 10 continue c call msg('Cartesian force field has been symmetrized.') return end c ====================================================================== subroutine symmetrizeg(NAT0,x,gx,nat,ne,pg,ity) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION gx(3*NAT0),x(3*NAT0) character*80 ts parameter (nsy0=13) character*3 pg(nsy0) DIMENSION ity(NAT0) DIMENSION a2(nsy0,3,3) CHARACTER*3 sl(nsy0) c call filla(nsy0,a2,sl) c tol=1.0d-2 open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.'S-TOL')read(10,*)tol if(ts(1:5).eq.'NOSYM')then close(10) return endif goto 5 999 close(10) c do 10 ie=1,ne do 11 is=1,nsy0 if(pg(ie).eq.sl(is))then do 1 ia=1,nat do 1 ja=ia,nat dist=sqrt((ax(is,x,a2,nsy0,ia,1)-x(3*(ja-1)+1))**2+ 1 (ax(is,x,a2,nsy0,ia,2)-x(3*(ja-1)+2))**2+ 2 (ax(is,x,a2,nsy0,ia,3)-x(3*(ja-1)+3))**2) if(dist.lt.tol.and.ity(ia).eq.ity(ja))then gx(3*(ja-1)+1)= ax(is,gx,a2,nsy0,ia,1) gx(3*(ja-1)+2)= ax(is,gx,a2,nsy0,ia,2) gx(3*(ja-1)+3)= ax(is,gx,a2,nsy0,ia,3) endif 1 continue endif 11 continue 10 continue c call msg('Cartesian gradient has been symmetrized.') return end c ====================================================================== subroutine newsx(NAT0,s,dq,nm,nat,io,am,sx) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION dq(3*NAT0),io(3*NAT0),s(3*NAT0,3*NAT0), 1am(3*NAT0),sx(3*NAT0) c do 1 i=1,3*nat xi=0.0d0 do 2 j=1,nm 2 if(io(j).eq.1)xi=xi-s(i,j)*dq(j) 1 sx(i)=xi/am(i) call msg('New cartesian step produced.') return end c ====================================================================== subroutine newx(NAT0,x,nat,sx) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION x(3*NAT0),sx(3*NAT0) c bohr=0.529177249d0 do 1 i=1,3*nat 1 x(i)=x(i)+sx(i)*bohr call msg('New coordinates produced.') return end c ====================================================================== subroutine report(NAT0,gn,w,io,nm,dq,e,glim,lupdate,lff,fq,lm, 1wmin,wmax,fx,wt,st,EDIAG,qt,gt,x,nat,ity,ne,pg,s) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) LOGICAL*4 symviol DIMENSION gn(3*NAT0),w(3*NAT0),io(3*NAT0),dq(3*NAT0), 1fq(3*NAT0,3*NAT0),fx(3*NAT0,3*NAT0),x(3*NAT0),ity(NAT0), 2wt(3*NAT0),st(3*NAT0,3*NAT0),EDIAG(3*NAT0),qt(3*NAT0),gt(3*NAT0), 3s(3*NAT0,3*NAT0) CHARACTER*7 state(2) CHARACTER*80 ts parameter (nsy0=13) CHARACTER*3 pg(nsy0) logical*4 lupdate,lff,lchkg,lch2,lrfo,lm data state/' fixed','relaxed'/ c qmax=0.8d0 lchkg=.false. lch2=.true. lrfo=.true. open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:6).eq.'allRFO')read(10,*)lrfo if(ts(1:4).eq.'qmax')read(10,*)qmax if(ts(1:9).eq.'checkgrad')read(10,*)lchkg if(ts(1:14).eq.'checktotalgrad')read(10,*)lch2 goto 5 999 close(10) c c SI constants: amu=1.6605402D-27 a0=0.529177249d-10 hartree=4.3597482d-18 clight=2.99792458d8 pi=3.141592653589d0 c constant that transfers cm-1 to au const=0.01d0/(2.0d0*pi*clight)*sqrt(hartree/(a0**2*amu)) c c produce new step: if(lm)then do 4 i=1,nm 4 if(w(i).ge.wmin.and.w(i).le.wmax)io(i)=0 endif c gmax=0.0d0 gnorm=0.0d0 do 14 i=1,nm gni=gn(i) if(io(i).eq.1)then if(abs(gni).gt.gmax)gmax=abs(gni) gnorm=gnorm+gni**2 endif 14 continue gnorm=sqrt(gnorm) c cosinus=0.0d0 dqn=0.0d0 IF(lupdate.and..not.lff)THEN c do 8 i=1,nm do 8 j=1,nm 8 fx(i,j)=fq(i,j) CALL TRED12(3*NAT0,nm,fx,st,wt,2,IERR,EDIAG) IF(IERR.NE.0)call bye(' Diagonalization error') c c Gradient in the eigenspace do 15 i=1,nm gt(i)=0.0d0 do 15 j=1,nm 15 gt(i)=gt(i)+st(j,i)*gn(j) c c call msg('Movement in the eigenspace') ems=1.0d-5 do 16 i=1,nm qt(i)=0.0d0 c if(symviol(i,st,x,nat,NAT0,ity,ne,pg))then c write(3,3335)i c else write(3,3336)i dd=wt(i)+sqrt(wt(i)**2+4.0d0*gt(i)**2) if(dd.gt.ems)qt(i)=-2.0d0*gt(i)/dd c endif 16 continue c call msg('Internal Hessian eigenvalues (cm-1):') do 6 i=1,nm t=wt(i) wt(i)=const*sqrt(abs(t)) 6 if(t.lt.0.0d0)wt(i)=-wt(i) write(3,308)(wt(i),i=1,nm) 308 format(8f10.2) c call msg('Movement in the normal mode space') do 11 i=1,nm dqi=0.0d0 if(symviol(i,s,x,nat,NAT0,ity,ne,pg))then write(3,3335)i else write(3,3336)i if(io(i).eq.1)then do 3 j=1,nm 3 dqi=dqi+st(i,j)*qt(j) cosinus=cosinus+dqi*gn(i) dqn=dqn+dqi**2 endif if(abs(dqi).gt.qmax)then dqi=dqi/abs(dqi)*qmax write(3,3334)i endif endif 11 dq(i)=dqi c ELSE do 2 i=1,nm dqi=0.0d0 gni=gn(i) if(io(i).eq.1)then Qi=(w(i)/const)**2 if(w(i).lt.0.0d0)Qi=-Qi if(lrfo.or.Qi.lt.0.0d0)then dd=Qi+sqrt(Qi**2+4.0d0*gni**2) if(dd.ne.0.0d0)dqi=-2.0d0*gni/dd else if(Qi.gt.1.0d-9)dqi=-gni/Qi endif c c check that the step is in the direction of gradient if(dqi*gni.gt.0.0d0.and.lchkg)then dqi=-qmax write(3,3333)i 3333 format(' Mode ',i5,' step not along gradient, step reversed') endif c c check, if symmetry violation would arise if(symviol(i,s,x,nat,NAT0,ity,ne,pg))then write(3,3335)i 3335 format(' Mode ',i5,' symmetry violation, do not move') dqi=0.0d0 else write(3,3336)i 3336 format(' Mode ',i5,' symmetric') endif c if(abs(dqi).gt.qmax)then dqi=dqi/abs(dqi)*qmax write(3,3334)i 3334 format(' Mode ',i5,' maximum step exceeded') endif c cosinus=cosinus+dqi*gni dqn=dqn+dqi**2 endif 2 dq(i)=dqi ENDIF dqn=sqrt(dqn) if(dqn*gnorm.gt.0.0d0)cosinus=cosinus/(dqn*gnorm) write(3,3009)cosinus 3009 format('Cos grad/dq ',f15.8) if(cosinus.gt.0.0d0.and.lch2)then call msg('Step reversed') do 1009 i=1,nm 1009 dq(i)=-dq(i) endif c write(3,100) 100 format(/,' ','Mode',' Frequency',' ',' Status', 1' Gradient', 2' Second Derivative', 3' Predicted Shift') write(3,101) 101 format(80(1h-)) do 1 i=1,nm if(lupdate.and..not.lff)then qi=fq(i,i) else qi=(w(i)/const)**2 endif 1 write(3,102)i,w(i),state(io(i)+1),gn(i),qi,dq(i) 102 format(i5,f12.2,1x,a7,3f18.8) write(3,101) write(3,*) write(3,103)e,gnorm,gmax,glim 103 format(' Energy : ',f15.8,' Norm of gradient: ',f15.8,/, 1 ' Max grad : ',f15.8,' Limit norm grad : ',f15.8,/) if(gnorm.gt.glim)then call msg('Not converged') else call bye('Converged') stop endif open(10,file='SQ.SCR',form='unformatted') write(10)(dq(i),i=1,nm) close(10) call msg('Normal step written into SQ.SCR.') return end c ====================================================================== subroutine writeng(NAT0,gn,nm) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION gn(3*NAT0) c open(10,file='QGRADIENT.SCR',form='unformatted') write(10)(gn(i),i=1,nm) close(10) call msg('Normal gradient saved to QGRADIENT.SCR.') return end c ====================================================================== subroutine normalg(NAT0,gn,gx,s,nat,nm,am) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION gx(3*NAT0),gn(3*NAT0),s(3*NAT0,3*NAT0),am(3*NAT0) c do 1 i=1,nm sum=0.0d0 do 2 j=1,3*nat 2 sum=sum+gx(j)*s(j,i)/am(j) 1 gn(i)=sum call msg('Normal gradient made') return end c ====================================================================== subroutine wx(NAT0,x,nat,ity) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION x(3*NAT0),ity(NAT0) c write(3,300) 300 format(/,' Cartesian coordinates:', 1/,25x,'Z-Matrix orientation:', 2/,1x,69(1h-), 3/,' Center Atomic Atomic',14x,'Coordinates (Angstroms)', 4/,' Number Number Type',14x,'X Y Z', 5/,1x,69(1h-)) do 1 i=1,nat 1 write(3,3001)i,ity(i),0,(x(3*(i-1)+j),j=1,3) 3001 format(i5,i11,i14,4x,3f12.6) write(3,301) 301 format(1x,69(1h-)) return end c ====================================================================== subroutine loadx(NAT0,x,nat,ity,nq,xq,cq) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION x(3*NAT0),ity(NAT0),xq(*),cq(*) character*80 filename,ts character*4 s4 character*2 sy logical*4 lex,lcharges,ldalton,lturbo,lgauss,ltinker,lcct c nq=0 lcharges=.false. ldalton=.false. lcct=.false. ltinker=.false. lturbo=.false. lgauss=.true. filename='G98.OUT' open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:11).eq.'coordinates')read(10,1000)filename if(ts(1:4).eq.'turb')read(10,*)lturbo if(ts(1:4).eq.'dalt')read(10,*)ldalton if(ts(1:4).eq.'lcct')read(10,*)lcct if(ts(1:4).eq.'tink')read(10,*)ltinker if(ts(1:6).eq.'CHARGE')read(10,*)lcharges goto 1 999 close(10) if(ldalton)filename='d.out' if(lcct)filename='FILE.X' if(ltinker)filename='tinker.xyz' if(lturbo)filename='coord' if(lturbo.or.ldalton.or.ltinker.or.lcct)lgauss=.false. c inquire(file=filename,exist=lex) if(.not.lex)goto 9999 open(10,file=filename,status='old') c cctn: if(lcct)then read(10,*) read(10,*) do 10 ia=1,nat 10 read(10,*)ity(ia),(x(3*(ia-1)+ix),ix=1,3) endif C DALTON: if(ldalton)then 21 read(10,1000,end=899,err=899)ts if(ts(26:48).eq.'Molecular geometry (au)')then read(10,*) read(10,*) do 199 ia=1,nat read(10,1000)ts sy=ts(1:2) read(ts(10:80),*)(x(3*(ia-1)+ix),ix=1,3) ity(ia)=0 if(sy.eq.' H')ity(ia)=1 if(sy.eq.' N')ity(ia)=7 if(sy.eq.' O')ity(ia)=8 if(sy.eq.' C')ity(ia)=6 do 199 ix=1,3 199 x(3*(ia-1)+ix)=x(3*(ia-1)+ix)*0.529177d0 goto 899 endif goto 21 endif C TINKER if(ltinker)then read(10,*) do 8 ia=1,nat read(10,1000)ts sy=ts(8:9) read(ts(12:47),*)(x(3*(ia-1)+ix),ix=1,3) ity(ia)=0 if(sy.eq.' H')ity(ia)=1 if(sy.eq.' N')ity(ia)=7 if(sy.eq.' O')ity(ia)=8 if(sy.eq.' C')ity(ia)=6 8 if(ity(ia).eq.0)write(6,*)'unknown atom' endif C TURBOMOLE if(lturbo)then read(10,*) do 9 ia=1,nat read(10,1009)(x(3*(ia-1)+ix),ix=1,3),sy 1009 format(f20.14,2f22.14,5x,a2) ity(ia)=0 if(sy.eq.' h')ity(ia)=1 if(sy.eq.' n')ity(ia)=7 if(sy.eq.' o')ity(ia)=8 if(sy.eq.' c')ity(ia)=6 do 9 ix=1,3 9 x(3*(ia-1)+ix)=x(3*(ia-1)+ix)*0.529177d0 endif C GAUSSIAN: if(lgauss)then if(filename.eq.'G98.OUT')then c gaussian output 2 read(10,1000,end=899,err=899)ts if(ts(27:44).eq.'Input orientation:'.or. 1 ts(26:46).eq.'Z-Matrix orientation:')then do 4 i=1,4 4 read(10,*) do 5 i=1,nat read(10,*)ity(i),ity(i),x(3*(i-1)+1),(x(3*(i-1)+ix),ix=1,3) 5 if(ity(i).lt.0) 1 read(10,*)ity(i),ity(i),x(3*(i-1)+1),(x(3*(i-1)+ix),ix=1,3) endif if(lcharges.and.ts(2:15).eq.'Point Charges:')then i=0 88 read(10,112)s4 112 format(1x,a4) if(s4.eq.'XYZ=')then backspace 10 i=i+1 if(i.gt.NAT0)call bye('Too many charges.') read(10,10012)(xq(3*(i-1)+ix),ix=1,3),cq(i) 10012 format(5x,3f10.4,3x,f10.4) goto 88 endif nq=i goto 899 endif goto 2 else c FILE.X format read(10,*) read(10,*)natx if(natx.ne.nat)call bye('loadx natx<>nat') do 6 ia=1,nat 6 read(10,*)ity(ia),(x(3*(ia-1)+ix),ix=1,3) endif endif 899 close(10) call msg('Coordinates loaded from '//filename//'.') if(lcharges.and.nq.eq.0)then call bye('No charges found in '//filename//'.') else call msg('Charges:') do 7 i=1,nq 7 write(3,301)i,(xq(3*(i-1)+ix),ix=1,3),cq(i) 301 format(i5,3f10.4,' Q = ',f10.4) endif return 9999 call bye('Coordinates not found in '//filename//'.') end c ====================================================================== subroutine setp(x,nat,nq,xq,cq,pot) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION x(*),xq(*),cq(*),pot(*) logical lex c if(nq.lt.1)return inquire(file='NUCPOT.SCR',exist=lex) if(lex)then open(11,file='NUCPOT.SCR',form='unformatted') read(11)(pot(ia),ia=1,nat) close(11) else call msg('Setting potential at nuclei.') call gcp(nat,pot,x,xq,cq,nq) open(11,file='NUCPOT.SCR',form='unformatted') write(11)(pot(ia),ia=1,nat) close(11) endif return end c ====================================================================== function exe(ts) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) character*(*) ts ns=len(ts) do 1 i1=1,ns 1 if(ts(i1:i1).eq.'=')goto 2 2 do 3 i2=i1,ns 3 if(ts(i2:i2).eq.'A')goto 4 4 read(ts(i1+1:i2-1),*)e exe=e return end c ====================================================================== subroutine writeg(NAT0,gx,nat,npar,gg) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION gx(3*NAT0),gg(3*NAT0) c open(10,file='GRADIENT.SCR',form='unformatted') if(npar.gt.0)write(10)(gx(i)-gg(i),i=1,3*nat) if(npar.eq.0)write(10)(gx(i),i=1,3*nat) close(10) call msg('Cartesian gradient saved to GRADIENT.SCR (binary).') return end c ====================================================================== subroutine loadg(NAT0,gx,nat,io,glim,e,wmin,wmax,lm) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION gx(3*NAT0),io(3*NAT0) parameter (imax=10) character*80 filename,ts character*1 tci(imax),ftci(imax),tcie(imax) logical*4 lm,ldalton,lturbo,lgauss,ltinker,lcct,lexc, 1lfollow integer*4 aci(imax),bci(imax),faci(imax),fbci(imax), 1acie(imax),bcie(imax) real*8 cci(imax),fcci(imax),ccie(imax) logical rfix c load=0 loade=0 lm=.false. ldalton=.false. lcct=.false. ltinker=.false. lturbo=.false. lexc=.false. lgauss=.true. wmin=-1.0d30 wmax= 1.0d30 rfix=.false. do 11 i=1,3*NAT0 11 io(i)=1 glim=0.0001d0 open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.'gradx')read(10,1000)filename if(ts(1:11).eq.'fixed range')then read(10,*)istart,iend do 131 i=istart,iend 131 io(i)=0 close(11) endif if(ts(1:11).eq.'modes fixed')then read(10,*)ix do 12 i=1,ix read(10,*)im 12 io(im)=0 endif if(ts(1:4).eq.'wmin')then read(10,*)wmin lm=.true. endif if(ts(1:4).eq.'rfix')read(10,*)rfix if(ts(1:4).eq.'dalt')read(10,*)ldalton if(ts(1:4).eq.'dalt')read(10,*)ldalton if(ts(1:4).eq.'lcct')read(10,*)lcct if(ts(1:4).eq.'tink')read(10,*)ltinker if(ts(1:4).eq.'turb')read(10,*)lturbo if(ts(1:4).eq.'lexc')read(10,*)lexc if(ts(1:4).eq.'wmax')then read(10,*)wmax lm=.true. endif if(ts(1:6).eq.'glimit')read(10,*)glim goto 1 999 close(10) if(lturbo.or.ldalton.or.ltinker.or.lcct)lgauss=.false. c c CCTN: if(lcct)then open(10,file='FILE.X') read(10,*) read(10,*)nat close(10) open(10,file='FILE.GR') do 6 ia=1,nat 6 read(10,*)(gx(3*(ia-1)+ix),ix=1,3) close(10) e=0.0d0 load=1 loade=1 endif C DALTON: if(ldalton)then open(10,file='d.out',status='old') 21 read(10,1000,end=899,err=899)ts if(ts(3:30).eq.'Total number of coordinates:')then read(ts(31:40),*)nat nat=nat/3 if(nat.gt.nat0)call bye('Too many atoms') endif if(ts(6:23).eq.'Nuclear repulsion:')then write(3,1000)ts read(ts(27:50),*)enuc read(10,1000)ts write(3,1000)ts read(ts(27:50),*)e e=e+enuc loade=1 endif if(ts(26:48).eq.'Molecular gradient (au)')then c try to jump over the symmetry excersizes: ia=1 901 read(10,1000,end=899)ts idp=0 do 902 i=1,80 902 if(ts(i:i).eq.'.')idp=idp+1 if(idp.ne.3)goto 901 read(ts(8:80),*)(gx(3*(ia-1)+ix),ix=1,3) do 903 ia=2,nat read(10,1000)ts 903 read(ts(8:80),*)(gx(3*(ia-1)+ix),ix=1,3) do 904 ix=1,3*nat 904 gx(ix)=-gx(ix) load=1 endif goto 21 endif c C TINKER: if(ltinker)then open(10,file='tinker.xyz') read(10,*)nat close(10) if(nat.gt.nat0)call bye('Too many atoms') open(10,file='GRAD.LST') c read energy kcal/mol: read(10,1000,end=899,err=899)ts read(ts,*)e c kcal/mol to hartree: fac=627.509469d0 e=e/fac loade=1 write(3,1000)ts c read gradient: do 905 ia=1,nat read(10,*)gx(3*(ia-1)+1),(gx(3*(ia-1)+ix),ix=1,3) c make forces from gradient and convert to kcal/mol/A: do 905 ix=1,3 905 gx(3*(ia-1)+ix)=-gx(3*(ia-1)+ix)/fac load=1 close(10) endif c C TURBOMOLE: if(lturbo)then open(10,file='gradient',status='old') read(10,*) read(10,1001)e 1001 format(32x,f18.10) loade=1 c c How many atoms: ia=0 211 read(10,*,end=8991,err=8991)xx,xx,xx ia=ia+1 if(ia.gt.nat0)call bye('Too many atoms') goto 211 8991 nat=ia/2 rewind 10 read(10,*) read(10,*) do 4 ia=1,nat 4 read(10,*) do 41 ia=1,nat 41 read(10,*)(gx(3*(ia-1)+j),j=1,3) do 9041 ix=1,3*nat 9041 gx(ix)=-gx(ix) close(10) load=1 endif C GAUSSIAN: if(lgauss)then c current excited state, and its scf and energy: ice=0 eeau=0.0d0 escf=0.0d0 inowf=0 faci(1)=0 fbci(1)=0 fcci(1)=0.0d0 if(lexc)then inquire(file='FOLLOW.SCR',exist=lfollow) if(lfollow)then call rwfollow(1,inowf,faci,fbci,fcci,ftci) spmax=0.0d0 ispmax=0 write(3,*)'Excited state from FOLLOW.SCR' endif endif filename='G98.OUT' open(10,file=filename,status='old') 2 read(10,1000,end=899,err=899)ts if(lexc)then c read SCF only to determine total energy for G09: if(ts(1:10).eq.' SCF Done:')escf=exe(ts) if(ts(1:15).eq.' Excited State ')then read(ts(15:18),*)item do 203 i=1,len(ts)-1 203 if(ts(i:i+1).eq.'eV')read(ts(i-10:i-2),*)ev c c record CI coefficients: call readcicoef(10,inow,imax,aci,bci,cci,tci) c identify state best matching the wanted one: if(lfollow)then sp=0.0d0 do 207 i=1,inow do 207 j=1,inowf 207 if(faci(j).eq.aci(i).and.fbci(j).eq.bci(i).and. 1 ftci(j).eq.tci(i)) 1 sp=sp+fcci(j)*cci(i) if(dabs(sp).gt.spmax)then spmax=dabs(sp) ispmax=item endif endif endif if(ts(1:28).eq.' This state for optimization')then ice=item eeau=escf+ev/27.211384205943d0 acie=aci bcie=bci ccie=cci tcie=tci inowe=inow endif c TDDFT excited state energy - for G16 if(ts(1:15).eq.' Total Energy, ')then read(ts(33:len(ts)),*)e write(3,1000)ts loade=1 endif else if(ts(2:29).eq.'ONIOM: extrapolated energy =')then write(3,1000)ts read(ts(30:80),*)e loade=1 endif if(ts(1:8).eq.' Energy='.and.ts(28:32).eq.'NIter')then write(3,1000)ts read(ts(9:27),*)e loade=1 endif if(ts(28:34).eq.'EUMP2 =')then write(3,1000)ts read(ts(35:len(ts)),*)e loade=1 endif if(ts(1:10).eq.' SCF Done:')then write(3,1000)ts e=exe(ts) loade=1 endif endif if(ts(38:59).eq.'Forces (Hartrees/Bohr)')then read(10,*) read(10,*) load=1 ia=0 3 read(10,1000)ts c c transition vector itrv=0 if(ts(1:17).eq.' -2 ')itrv=1 if(ts(2:3).ne.'--'.and.itrv.eq.0)then ia=ia+1 if(ia.gt.nat0)then call bye('Too many atoms') stop endif backspace 10 read(10,*)(gx(3*(ia-1)+ix),ix=1,2),(gx(3*(ia-1)+ix),ix=1,3) goto 3 endif nat=ia endif goto 2 endif 899 close(10) if(load.eq.0)then call bye('Gradient not found') else call msg('Cartesian gradient loaded.') cnorm=0.0d0 cmax=gx(1) do 5 i=1,3*nat cnorm=cnorm+gx(i)**2 5 if(abs(gx(i)).gt.abs(cmax))cmax=gx(i) cnorm=sqrt(cnorm) rms=cnorm/sqrt(dble(3*nat)) write(3,3003)cnorm,cmax,rms 3003 format(' Norm ',f15.8,/,' Max ',f15.8,/,' RMS ',f15.8,/) endif if(lexc)then if(loade.eq.0)then write(3,3010) 3010 format(' No G16 excited state energy') if(ice.ne.0)then e=eeau write(3,3011)e 3011 format(' G09 like excited state energy ',f16.9) else call bye('Excited state energy not found') endif endif c check that we are looking at the desired state: if(lfollow)then write(3,304)ispmax,ice,spmax 304 format(' State desired : ',i4,/, 1 ' State computed: ',i4,/, 1 ' Scalar product: ',f10.4) if(ispmax.ne.ice.and..not.rfix)then write(3,*)'INP.NEW will be redone' call ginproot(ispmax) else call rwfollow(0,inowe,acie,bcie,ccie,tcie) if(ispmax.ne.ice)write(6,*)'state switch ignored' write(3,*)'new state characteristic recorded' endif else call rwfollow(0,inowe,acie,bcie,ccie,tcie) write(3,*)'new state characteristic recorded' endif else if(loade.eq.0)call bye('Energy not found') endif return end c ====================================================================== subroutine projectgr(gx,nat,x) implicit none integer*4 nat real*8 gx(*),x(*) character*4 ts logical*4 lprg lprg=.true. open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a4) if(ts.eq.'gpro')read(10,*)lprg goto 1 999 close(10) if(lprg)call PROJECTG(gx,nat,x) return end c ====================================================================== subroutine loads(NAT0,s,nat,w,nm) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION s(3*NAT0,3*NAT0),w(3*NAT0) character*80 filename,ts logical*4 li0 c filename='F.INP' open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:8).eq.'s-matrix')then read(10,1000)filename goto 999 endif goto 1 999 close(10) c open(10,file=filename) read(10,*,end=9999,err=9999)nm do 2 i=1,nat+1 2 read(10,*) do 3 ia=1,nat do 3 im=1,nm 3 read(10,*)(s((ia-1)*3+ix,im),ix=1,2),(s((ia-1)*3+ix,im),ix=1,3) read(10,*)nm read(10,*)(w(im),im=1,nm) close(10) call msg('S matrix loaded from '//filename//'.') c li0=.true. open(10,file='Q.OPT') 11 read(10,1000,end=699,err=699)ts if(ts(1:11).eq.'ignore zero')read(10,*)li0 goto 11 699 close(10) if(li0)then call msg('Zero modes ignored and deleted') i=0 411 i=i+1 if(abs(w(i)).lt.1.0d0)then do 412 ii=i,nm-1 w(ii)=w(ii+1) do 412 ix=1,3*nat 412 s(ix,ii)=s(ix,ii+1) nm=nm-1 i=i-1 endif if(i.lt.nm)goto 411 endif c return 9999 call bye(filename//'not found') end c ====================================================================== subroutine loadm(NAT0,am,nat) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION am(3*NAT0) character*80 filename,ts c filename='FTRY.INP' open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:6).eq.'masses')then read(10,1000)filename goto 999 endif goto 1 999 close(10) c open(10,file=filename) read(10,*,end=9999,err=9999) read(10,*)nsf do 2 i=1,nsf read(10,*) 2 read(10,*) read(10,444)(qdum,i=1,nat) 444 format(6f12.6) read(10,*) read(10,444)(am((i-1)*3+1),i=1,nat) do 3 i=1,nat am((i-1)*3+1)=sqrt(am((i-1)*3+1)) am((i-1)*3+2)= am((i-1)*3+1) 3 am((i-1)*3+3)= am((i-1)*3+1) close(10) call msg('Atomic masses loaded from '//filename) return 9999 call bye(filename//' not found') end c ====================================================================== subroutine dimensionlesss(NAT0,s,am,nat,nm) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION am(3*NAT0),s(3*NAT0,3*NAT0) c do 1 i=1,3*nat do 1 j=1,nm 1 s(i,j)=s(i,j)*am(i) return end c ====================================================================== subroutine bye(ts) character*(*) ts write(3,6001) write(3,'(a)')ts write(3,6001) 6001 format(80(1h*)) write(3,6000) 6000 format('program stopped') close(3) stop end c ====================================================================== subroutine makes(NAT0,fx,s,nat,w,nm,am,A,TEM,ity,x,fq,ZM,C,EDIAG) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION fx(3*NAT0,3*NAT0),s(3*NAT0,3*NAT0),w(3*NAT0),am(3*NAT0), 1ZM(3*NAT0),A(3*NAT0,3*NAT0),ity(NAT0),x(3*NAT0),fq(3*NAT0,3*NAT0), 1TEM(3*NAT0,3*NAT0),C(3,NAT0),EDIAG(3*NAT0) LOGICAL*4 lsout,li0,lpr CHARACTER*80 ts c lsout=.false. li0=.true. lpr=.true. open(10,file='Q.OPT') 11 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:7).eq.'smatout')read(10,1000)lsout if(ts(1:7).eq.'project')read(10,1000)lpr if(ts(1:11).eq.'ignore zero')read(10,*)li0 goto 11 999 close(10) DO 220 i=1,3*nat 220 ZM(i)=1/am(i) C cc=4.3597482d0/0.529177249d0/0.529177249d0 do 221 i=1,3*nat do 221 j=1,3*nat 221 fq(i,j)=fx(i,j)*cc c if(lpr)CALL PROJECT(fq,NAT0,nat,A,TEM,C,x) C CALL FAST(NAT0,nat,fq,ZM) CALL TRED12(3*NAT0,3*nat,fq,s,w,2,IERR,EDIAG) IF(IERR.NE.0)call bye(' Diagonalization error') nm=3*nat const=sqrt(cc) DO 310 I=1,nm T=w(I) SIGN=1.0d0 IF(T.LT.0.0d0)SIGN=-1.0d0 T=SIGN*SQRT(ABS(T))*5140.48707d0/const 310 w(I)=T C C Calculate the true mass-weighted S-matrix: DO 210 IA=1,3*nat DO 210 IM=1,nm 210 s(IA,IM)=ZM(IA)*s(IA,IM) c if(li0)then call msg('Zero modes ignored and deleted') i=0 411 i=i+1 if(abs(w(i)).lt.1.0d0)then do 412 ii=i,nm-1 w(ii)=w(ii+1) do 412 ix=1,3*nat 412 s(ix,ii)=s(ix,ii+1) nm=nm-1 i=i-1 endif if(i.lt.nm)goto 411 endif if(lsout)call SMATOUT(s,w,NAT0,nat,nm,ity,x) call msg('S matrix made from force field.') return END C SUBROUTINE SMATOUT(s,w,NAT0,nat,nm,ity,x) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION S(3*NAT0,3*NAT0),w(NAT0),ity(NAT0),x(3*NAT0) OPEN(34,FILE='F.INP') WRITE(34,10)nm,3*nat,nat 10 FORMAT(3I5) DO 3 I=1,nat 3 WRITE(34,11)ity(I),(x(3*(I-1)+ix),ix=1,3) 11 FORMAT(I4,3F12.6) WRITE(34,14) 14 FORMAT(' Atom Mode X-disp. Y-disp. Z-disp.') DO 1 I=1,nat DO 1 J=1,nm IU=3*(I-1) 1 WRITE(34,15)I,J,(s(IU+ix,J),ix=1,3) 15 FORMAT(2I5,3F11.6) WRITE(34,11)nm WRITE(34,13)(w(I),I=1,nm) 13 format(6F11.3) CLOSE(34) RETURN END SUBROUTINE TRED12(MX3,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (MX3,N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(MX3,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END C SUBROUTINE WRITEFF(NAT0,nat,fx,fn) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION fx(3*NAT0,3*NAT0) character*(*) fn c open(10,file=fn) N=3*nat N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(10,17)LN,(fx(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) close(10) call msg('Force field written.') RETURN END SUBROUTINE FAST(NAT0,nat,fq,ZM) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION fq(3*NAT0,3*nat),ZM(3*nat) DO 214 I=1,3*nat DO 214 J=I,3*nat fq(I,J)=fq(I,J)*ZM(J)*ZM(I) 214 fq(J,I)=fq(I,J) RETURN END C SUBROUTINE PROJECT(F,NAT0,nat,A,TEM,C,x) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION F(3*NAT0,3*NAT0),A(3*NAT0,3*NAT0),TEM(3*NAT0,3*NAT0), 1C(3,NAT0),x(3*NAT0) C WRITE(3,*)'Projecting transl/rotations from the force field...' CALL DOMA(A,NAT0,nat,TEM,C,x) N=3*nat DO 1 I=1,N DO 1 J=1,N TEM(I,J)=0.0d0 DO 1 II=1,N 1 TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J) DO 2 I=1,N DO 2 J=1,N F(I,J)=0.0d0 DO 2 II=1,N 2 F(I,J)=F(I,J)+TEM(I,II)*A(J,II) RETURN END SUBROUTINE PROJECTG(G,nat,x) IMPLICIT none integer*4 NAT0 parameter (NAT0=900) integer*4 nat,N,I,II real*8 G(*),x(*),A(3*NAT0,3*NAT0),TEM(3*NAT0,3*NAT0), 1C(3,NAT0),t(3*NAT0) C WRITE(3,*)'Projecting transl/rotations from the gradient...' CALL DOMA(A,NAT0,nat,TEM,C,x) N=3*nat DO 1 I=1,N t(I)=0.0d0 DO 1 II=1,N 1 t(I)=t(I)+A(I,II)*G(II) DO 2 I=1,N 2 G(I)=t(I) RETURN END SUBROUTINE DOMA(A,NAT0,nat,TEM,C,x) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION C(3,NAT0),A(3*NAT0,3*nat),TEM(3*NAT0,3*nat), 1x(*) C AMACH=1.0d-14 DO 11 I=1,nat C(1,I)=x(3*(I-1)+1) C(2,I)=x(3*(I-1)+2) 11 C(3,I)=x(3*(I-1)+3) DO 12 I=1,3*nat DO 12 J=1,6 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.0d0 A(3*I-1,2)= 1.0d0 A(3*I ,3)= 1.0d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) N6=6 CALL ORT(A,NAT0,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,3*nat 6 S=S+A(J,ICOL)**2 IF(S.LT.1.0d3*AMACH)THEN DO 7 J=1,3*nat 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE WRITE(3,*)N6,' coordinates projected' DO 2 I=1,3*nat DO 2 J=1,3*nat TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,3*nat DO 4 J=1,3*nat 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END C ============================================================ SUBROUTINE ORT(A,NAT0,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(3*NAT0,N6) AMACH=1.0D-13 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END C ============================================================ subroutine msg(ts) character*(*) ts c do 1 i=min(80,len(ts)),1,-1 1 if(ts(i:i).ne.' ')goto 2 2 write(3,*)ts(1:i) return end C ============================================================ function sp(a,b) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension a(3),b(3) sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end C ============================================================ subroutine getpar(NAT0,x,kind,p,actual,npar0,ipar,ip,pp) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ipar(npar0,4),p(4,3),v1(3),v2(3),v3(3),x(3*NAT0), 1pp(4,3,4,3), 2delta(3,3),v1x(3,4,3),v2x(3,4,3),v3x(3,4,3) data delta/1.0d0,0.0d0,0.0d0,0.0d0,1.0d0,0.0d0,0.0d0,0.0d0,1.0d0/ c actual=0.0d0 fac=180.0d0/3.141592653589d0 if(kind.eq.2)then do 3 ix=1,3 3 v1(ix)=x(3*(ipar(ip,2)-1)+ix)-x(3*(ipar(ip,1)-1)+ix) v11=sp(v1,v1) c c v - derivatives do 1 ix=1,3 do 1 iv=1,3 v1x(iv,1,ix)=-delta(iv,ix) 1 v1x(iv,2,ix)= delta(iv,ix) c actual=dsqrt(v11) c do 2 ia=1,2 do 2 ix=1,3 v1xd1=v1x(1,ia,ix)*v1(1)+v1x(2,ia,ix)*v1(2)+v1x(3,ia,ix)*v1(3) p(ia,ix)=v1xd1/actual do 2 ja=1,2 do 2 jx=1,3 v1yd1=v1x(1,ja,jx)*v1(1)+v1x(2,ja,jx)*v1(2)+v1x(3,ja,jx)*v1(3) v1x1y=v1x(1,ia,ix)*v1x(1,ja,jx)+v1x(2,ia,ix)*v1x(2,ja,jx) 1 +v1x(3,ia,ix)*v1x(3,ja,jx) 2 pp(ia,ix,ja,jx)=-v1xd1*v1yd1/actual**3+v1x1y/actual endif if(kind.eq.3)then do 4 ix=1,3 v1(ix)=x(3*(ipar(ip,2)-1)+ix)-x(3*(ipar(ip,1)-1)+ix) 4 v2(ix)=x(3*(ipar(ip,3)-1)+ix)-x(3*(ipar(ip,2)-1)+ix) v11=sp(v1,v1) v12=sp(v1,v2) v22=sp(v2,v2) c c v - derivatives do 55 ix=1,3 do 55 iv=1,3 v1x(iv,1,ix)=-delta(iv,ix) v2x(iv,1,ix)= 0.0d0 v1x(iv,2,ix)= delta(iv,ix) v2x(iv,2,ix)=-delta(iv,ix) v1x(iv,3,ix)= 0.0d0 55 v2x(iv,3,ix)= delta(iv,ix) c ps=-v12 ah=dsqrt(v11) bh=dsqrt(v22) an=ah*bh o=ps/an actual=fac*acos(o) dp=dsqrt(1.0d0-o**2) c do 56 ia=1,3 do 56 ix=1,3 v1xd1=v1x(1,ia,ix)*v1(1)+v1x(2,ia,ix)*v1(2)+v1x(3,ia,ix)*v1(3) v1xd2=v1x(1,ia,ix)*v2(1)+v1x(2,ia,ix)*v2(2)+v1x(3,ia,ix)*v2(3) v2xd1=v2x(1,ia,ix)*v1(1)+v2x(2,ia,ix)*v1(2)+v2x(3,ia,ix)*v1(3) v2xd2=v2x(1,ia,ix)*v2(1)+v2x(2,ia,ix)*v2(2)+v2x(3,ia,ix)*v2(3) px=-v1xd2-v2xd1 ax=bh*v1xd1/ah+ah*v2xd2/bh c c Ox=Px/A-PAx/A**2 ox=px/an-ps*ax/an**2 p(ia,ix)=-fac*ox/dp do 56 ja=1,3 do 56 jx=1,3 v1yd1=v1x(1,ja,jx)*v1(1)+v1x(2,ja,jx)*v1(2)+v1x(3,ja,jx)*v1(3) v1yd2=v1x(1,ja,jx)*v2(1)+v1x(2,ja,jx)*v2(2)+v1x(3,ja,jx)*v2(3) v2yd1=v2x(1,ja,jx)*v1(1)+v2x(2,ja,jx)*v1(2)+v2x(3,ja,jx)*v1(3) v2yd2=v2x(1,ja,jx)*v2(1)+v2x(2,ja,jx)*v2(2)+v2x(3,ja,jx)*v2(3) v1x1y=v1x(1,ia,ix)*v1x(1,ja,jx)+v1x(2,ia,ix)*v1x(2,ja,jx) 1 +v1x(3,ia,ix)*v1x(3,ja,jx) v1x2y=v1x(1,ia,ix)*v2x(1,ja,jx)+v1x(2,ia,ix)*v2x(2,ja,jx) 1 +v1x(3,ia,ix)*v2x(3,ja,jx) v2x1y=v2x(1,ia,ix)*v1x(1,ja,jx)+v2x(2,ia,ix)*v1x(2,ja,jx) 1 +v2x(3,ia,ix)*v1x(3,ja,jx) v2x2y=v2x(1,ia,ix)*v2x(1,ja,jx)+v2x(2,ia,ix)*v2x(2,ja,jx) 1 +v2x(3,ia,ix)*v2x(3,ja,jx) py=-v1yd2-v2yd1 ay=bh*v1yd1/ah+ah*v2yd2/bh axy=-bh*v1xd1*v1yd1/ah**3+bh*v1x1y/ah+v1xd1*v2yd2/ah/bh 1 + v1yd1*v2xd2/ah/bh-ah*v2xd2*v2yd2/bh**3+ah*v2x2y/bh pxy=-v1x2y-v2x1y c oy=py/an-ps*ay/an**2 c Oxy=Pxy/N-PxNy/N**2-NxPy/N**2-PNxy/N**2+2PNxNy/N**3 oxy=pxy/an-px*ay/an**2- 1 ax*py/an**2-ps*axy/an**2+2.0d0*ps*ax*ay/an**3 56 pp(ia,ix,ja,jx)=-fac*oxy/dp-fac*o*ox*oy/dp**3 endif if(kind.eq.4)then do 5 ix=1,3 v1(ix)=x(3*(ipar(ip,2)-1)+ix)-x(3*(ipar(ip,1)-1)+ix) v2(ix)=x(3*(ipar(ip,3)-1)+ix)-x(3*(ipar(ip,2)-1)+ix) 5 v3(ix)=x(3*(ipar(ip,4)-1)+ix)-x(3*(ipar(ip,3)-1)+ix) v11=sp(v1,v1) v12=sp(v1,v2) v13=sp(v1,v3) v22=sp(v2,v2) v23=sp(v2,v3) v33=sp(v3,v3) c c v - derivatives do 53 ix=1,3 do 53 iv=1,3 v1x(iv,1,ix)=-delta(iv,ix) v2x(iv,1,ix)= 0.0d0 v3x(iv,1,ix)= 0.0d0 v1x(iv,2,ix)= delta(iv,ix) v2x(iv,2,ix)=-delta(iv,ix) v3x(iv,2,ix)= 0.0d0 v1x(iv,3,ix)= 0.0d0 v2x(iv,3,ix)= delta(iv,ix) v3x(iv,3,ix)=-delta(iv,ix) v1x(iv,4,ix)= 0.0d0 v2x(iv,4,ix)= 0.0d0 53 v3x(iv,4,ix)= delta(iv,ix) c c P=-v1.v3 v2**2 + v1*v2 v2*v3 ps=-v13*v22+v12*v23 c c A=aa**.5 bb**.5 aa=v11*v22-v12**2 bb=v22*v33-v23**2 ah=dsqrt(aa) bh=dsqrt(bb) an=ah*bh o=ps/an dp=dsqrt(1.0d0-o**2) c c figure out sign a1=v1(2)*v2(3)-v1(3)*v2(2) a2=v1(3)*v2(1)-v1(1)*v2(3) a3=v1(1)*v2(2)-v1(2)*v2(1) sig=1.0d0 if(a1*v3(1)+a2*v3(2)+a3*v3(3).lt.0.0d0)sig=-1.0d0 actual=fac*acos(o)*sig c c P (first and second) do 54 ia=1,4 do 54 ix=1,3 v1xd1=v1x(1,ia,ix)*v1(1)+v1x(2,ia,ix)*v1(2)+v1x(3,ia,ix)*v1(3) v1xd2=v1x(1,ia,ix)*v2(1)+v1x(2,ia,ix)*v2(2)+v1x(3,ia,ix)*v2(3) v1xd3=v1x(1,ia,ix)*v3(1)+v1x(2,ia,ix)*v3(2)+v1x(3,ia,ix)*v3(3) v2xd1=v2x(1,ia,ix)*v1(1)+v2x(2,ia,ix)*v1(2)+v2x(3,ia,ix)*v1(3) v2xd2=v2x(1,ia,ix)*v2(1)+v2x(2,ia,ix)*v2(2)+v2x(3,ia,ix)*v2(3) v2xd3=v2x(1,ia,ix)*v3(1)+v2x(2,ia,ix)*v3(2)+v2x(3,ia,ix)*v3(3) v3xd1=v3x(1,ia,ix)*v1(1)+v3x(2,ia,ix)*v1(2)+v3x(3,ia,ix)*v1(3) v3xd2=v3x(1,ia,ix)*v2(1)+v3x(2,ia,ix)*v2(2)+v3x(3,ia,ix)*v2(3) v3xd3=v3x(1,ia,ix)*v3(1)+v3x(2,ia,ix)*v3(2)+v3x(3,ia,ix)*v3(3) px=-v1xd3*v22-v3xd1*v22-2.0d0*v13*v2xd2 1 +v1xd2*v23+v2xd1*v23+v12*v2xd3+v12*v3xd2 aax=2.0d0*(v1xd1*v22+v2xd2*v11-v12*v2xd1-v12*v1xd2) bbx=2.0d0*(v2xd2*v33+v3xd3*v22-v23*v3xd2-v23*v2xd3) ax=0.5d0*(aax*bh/ah+bbx*ah/bh) c c Ox=Px/A-PAx/A**2 ox=px/an-ps*ax/an**2 if(abs(dp).lt.1.0d-20)then p(ia,ix)=0.0d0 else p(ia,ix)=-fac*ox/dp*sig endif do 54 ja=1,4 do 54 jx=1,3 v1yd1=v1x(1,ja,jx)*v1(1)+v1x(2,ja,jx)*v1(2)+v1x(3,ja,jx)*v1(3) v1yd2=v1x(1,ja,jx)*v2(1)+v1x(2,ja,jx)*v2(2)+v1x(3,ja,jx)*v2(3) v1yd3=v1x(1,ja,jx)*v3(1)+v1x(2,ja,jx)*v3(2)+v1x(3,ja,jx)*v3(3) v2yd1=v2x(1,ja,jx)*v1(1)+v2x(2,ja,jx)*v1(2)+v2x(3,ja,jx)*v1(3) v2yd2=v2x(1,ja,jx)*v2(1)+v2x(2,ja,jx)*v2(2)+v2x(3,ja,jx)*v2(3) v2yd3=v2x(1,ja,jx)*v3(1)+v2x(2,ja,jx)*v3(2)+v2x(3,ja,jx)*v3(3) v3yd1=v3x(1,ja,jx)*v1(1)+v3x(2,ja,jx)*v1(2)+v3x(3,ja,jx)*v1(3) v3yd2=v3x(1,ja,jx)*v2(1)+v3x(2,ja,jx)*v2(2)+v3x(3,ja,jx)*v2(3) v3yd3=v3x(1,ja,jx)*v3(1)+v3x(2,ja,jx)*v3(2)+v3x(3,ja,jx)*v3(3) v1x1y=v1x(1,ia,ix)*v1x(1,ja,jx)+v1x(2,ia,ix)*v1x(2,ja,jx) 1 +v1x(3,ia,ix)*v1x(3,ja,jx) v1x2y=v1x(1,ia,ix)*v2x(1,ja,jx)+v1x(2,ia,ix)*v2x(2,ja,jx) 1 +v1x(3,ia,ix)*v2x(3,ja,jx) v1x3y=v1x(1,ia,ix)*v3x(1,ja,jx)+v1x(2,ia,ix)*v3x(2,ja,jx) 1 +v1x(3,ia,ix)*v3x(3,ja,jx) v2x1y=v2x(1,ia,ix)*v1x(1,ja,jx)+v2x(2,ia,ix)*v1x(2,ja,jx) 1 +v2x(3,ia,ix)*v1x(3,ja,jx) v2x2y=v2x(1,ia,ix)*v2x(1,ja,jx)+v2x(2,ia,ix)*v2x(2,ja,jx) 1 +v2x(3,ia,ix)*v2x(3,ja,jx) v2x3y=v2x(1,ia,ix)*v3x(1,ja,jx)+v2x(2,ia,ix)*v3x(2,ja,jx) 1 +v2x(3,ia,ix)*v3x(3,ja,jx) v3x1y=v3x(1,ia,ix)*v1x(1,ja,jx)+v3x(2,ia,ix)*v1x(2,ja,jx) 1 +v3x(3,ia,ix)*v1x(3,ja,jx) v3x2y=v3x(1,ia,ix)*v2x(1,ja,jx)+v3x(2,ia,ix)*v2x(2,ja,jx) 1 +v3x(3,ia,ix)*v2x(3,ja,jx) v3x3y=v3x(1,ia,ix)*v3x(1,ja,jx)+v3x(2,ia,ix)*v3x(2,ja,jx) 1 +v3x(3,ia,ix)*v3x(3,ja,jx) aay=2.0d0*(v1yd1*v22+v2yd2*v11-v12*v2yd1-v12*v1yd2) bby=2.0d0*(v2yd2*v33+v3yd3*v22-v23*v3yd2-v23*v2yd3) py=-v1yd3*v22-v3yd1*v22-2.0d0*v13*v2yd2 1 +v1yd2*v23+v2yd1*v23+v12*v2yd3+v12*v3yd2 ay=0.5d0*(aay*bh/ah+bby*ah/bh) aaxy=2.0d0*(v1x1y*v22+2.0d0*v1xd1*v2yd2+v2x2y*v11 1 +2.0d0*v2xd2*v1yd1 1 -v2x1y*v12-v2xd1*v1yd2-v2xd1*v2yd1 2 -v1x2y*v12-v1xd2*v1yd2-v1xd2*v2yd1) bbxy=2.0d0*(v2x2y*v33+2.0d0*v2xd2*v3yd3+v3x3y*v22 1 +2.0d0*v3xd3*v2yd2 1 -v3x2y*v23-v3xd2*v2yd3-v3xd2*v3yd2 2 -v2x3y*v23-v2xd3*v2yd3-v2xd3*v3yd2) axy=0.5d0*ah*bh*(aaxy/aa 1 -0.5d0*(aax*aay/aa**2-(aax*bby+bbx*aay)/aa/bb+bbx*bby/bb**2) 2 +bbxy/bb) pxy=-v1x3y*v22-2.0d0*v1xd3*v2yd2-v3x1y*v22 1 -2.0d0*v3xd1*v2yd2 -2.0d0*v1yd3*v2xd2-2.0d0*v2xd2*v3yd1 2 -2.0d0*v13*v2x2y+v1x2y*v23+v1xd2*v2yd3+v1xd2*v3yd2 3 +v2x1y*v23+v2xd1*v2yd3+v2xd1*v3yd2 4 +v2xd3*v1yd2+v2yd1*v2xd3+v12*v2x3y 5 +v3xd2*v1yd2+v3xd2*v2yd1+v12*v3x2y oy=py/an-ps*ay/an**2 c c Oxy=Pxy/N-PxNy/N**2-NxPy/N**2-PNxy/N**2+2PNxNy/N**3 oxy=pxy/an-px*ay/an**2- 1 ax*py/an**2-ps*axy/an**2+2.0d0*ps*ax*ay/an**3 54 pp(ia,ix,ja,jx)=fac*sig*(-oxy/dp-o*ox*oy/dp**3) endif return end C ============================================================ subroutine loadbar(def,bar,tp) character*80 ts character*5 tp real*8 bar,def c bar=def open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.tp)read(10,*)bar goto 5 999 close(10) return end C ============================================================ subroutine loadpar(npar,ipar,kpar,npar0,bpar,apar,nsame,ijsame, 1bijsame) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ipar(npar0,4),kpar(npar0),bpar(npar0), 1apar(npar0),ijsame(npar0,2),bijsame(npar0,2) character*80 ts c npar=0 nsame=0 open(10,file='Q.OPT') 5 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:5).eq.'fixed')then read(10,*)npar if(npar.gt.npar0)then call msg('npar>npar0') close(10) close(3) stop endif do 17 ii=1,npar 17 read(10,*)kpar(ii),(ipar(ii,jj),jj=1,kpar(ii)),apar(ii),bpar(ii) endif if(ts(1:4).eq.'same')then read(10,*)nsame if(nsame.gt.npar0)then call msg('nsame>npar0') close(10) close(3) stop endif do 1 ii=1,nsame read(10,*)ic backspace 10 if(ic.gt.0)then call msg('Equal parameters') read(10,*)ijsame(ii,1),ijsame(ii,2),bijsame(ii,1) bijsame(ii,2)=0.0d0 else call msg('Constant difference') read(10,*)ijsame(ii,1),ijsame(ii,2),(bijsame(ii,jj),jj=1,2) ijsame(ii,1)=abs(ijsame(ii,1)) endif 1 continue endif goto 5 999 close(10) do 2 ii=1,nsame do 2 i2=1,2 ic=0 do 3 jj=1,npar 3 if(ijsame(ii,i2).eq.jj)ic=ic+1 2 if(ic.eq.0)call bye('Same parameter not defined') return end c ===================================== subroutine fixeg(npar,e,ipar,kpar,npar0,bpar,apar,NAT0,x,gg, 1nat,p,pp,ijsame,bijsame,nsame,p1,p2,pp1,pp2) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ipar(npar0,4),kpar(npar0),bpar(npar0), 1apar(npar0),x(3*NAT0),gg(3*NAT0),ijsame(npar0,2),bijsame(npar0,2), 2p(4,3),pp(4,3,4,3),p1(4,3),p2(4,3),pp1(4,3,4,3),pp2(4,3,4,3) character*10 sk(4) data sk/' unknown ',' distance ',' angle ',' torsion '/ call msg('Adding fixed parameter terms') bohr=0.529177d0 write(3,3001) 3001 format(' parameter barrier desired ', 1 ' current atoms') eold=e do 6 i=1,3*nat 6 gg(i)=0.0d0 c do 1 i=1,npar kind=kpar(i) a0=apar(i) B=bpar(i) call getpar(NAT0,x,kind,p,actual,npar0,ipar,i,pp) c c if over 180, accomodate the sign: if(abs(actual-a0).gt.180.0d0.and.kind.eq.4)then if(abs(actual-a0+360.0d0).lt.abs(actual-a0))then call msg('Changing reference angle -> a0-360') a0=a0-360.0d0 endif if(abs(actual-a0-360.0d0).lt.abs(actual-a0))then call msg('Changing reference angle -> a0+360') a0=a0+360.0d0 endif endif c c E=B.(actual-a0)**2 e=e+B*(actual-a0)**2 c c Ex=2B(actual-a0)actualx pfac=2.0d0*B*(actual-a0)*bohr do 61 ia=1,kind do 61 ix=1,3 ii=3*(ipar(i,ia)-1)+ix c c actually Forces, i.e. negative gradient needed 61 gg(ii)=gg(ii)-pfac*p(ia,ix) c 1 write(3,3000)i,sk(kind),B,a0,actual, 1(ipar(i,j),j=1,kind) 3000 format(i2,a10,f10.5,2(f12.4),4i4) c write(3,30041) 30041 format('Condition par1 par2 Difference Target') do 2 i=1,nsame B=bijsame(i,1) D=bijsame(i,2) i1=ijsame(i,1) i2=ijsame(i,2) kind1=kpar(i1) kind2=kpar(i2) if(kind1.ne.kind2)call msg(' WARNING: KIND1 <> KIND2 !') a01=apar(i1) a02=apar(i2) call getpar(NAT0,x,kind1,p1,actual1,npar0,ipar,i1,pp1) call getpar(NAT0,x,kind2,p2,actual2,npar0,ipar,i2,pp2) c c if over 180, accomodate the sign: if(abs(actual1-a01).gt.180.0d0.and.kind1.eq.4)then if(abs(actual1-a01+360.0d0).lt.abs(actual1-a01))then call msg('Changing reference angle -> a01-360') a01=a01-360.0d0 endif if(abs(actual1-a01-360.0d0).lt.abs(actual1-a01))then call msg('Changing reference angle -> a01+360') a01=a01+360.0d0 endif endif if(abs(actual2-a02).gt.180.0d0.and.kind2.eq.4)then if(abs(actual2-a02+360.0d0).lt.abs(actual2-a02))then call msg('Changing reference angle -> a02-360') a02=a02-360.0d0 endif if(abs(actual2-a02-360.0d0).lt.abs(actual2-a02))then call msg('Changing reference angle -> a02+360') a02=a02+360.0d0 endif endif c c E=B.(actual1-actual2-D)**2 e=e+B*(actual1-actual2-D)**2 c c Ex=2B(actual1-actual2-D)(actual1x-actual2x) c actually Forces, i.e. negative gradient used pfac=2.0d0*B*(actual1-actual2-D)*bohr do 63 ia=1,kind1 do 63 ix=1,3 ii=3*(ipar(i1,ia)-1)+ix 63 gg(ii)=gg(ii)-pfac*p1(ia,ix) do 64 ia=1,kind2 do 64 ix=1,3 ii=3*(ipar(i2,ia)-1)+ix 64 gg(ii)=gg(ii)+pfac*p2(ia,ix) c 2 write(3,3004)i,actual1,actual2,actual1-actual2,D 3004 format(i2,' par1: ',f12.4,' par2: ',3f12.4) c write(3,3002)e,eold 3002 format(' Total and molecular energy: ',2f20.6,/, 1'Parameter Force (Hartree/bohr):') do 62 i=1,nat 62 write(3,3003)i,(gg(3*(i-1)+j),j=1,3) 3003 format(i4,3f12.6) return end c ===================================== subroutine fixegc(e,x,gg,nat,pot,potc,xq,cq,nq) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension x(*),gg(*),pot(*),potc(*),xq(*),cq(*) call msg('Adding charge terms') bohr=0.529177d0 eold=e do 6 i=1,3*nat 6 gg(i)=0.0d0 call loadbar(1.0d0,ALPHA,'ALPHA') c c current potential call gcp(nat,potc,x,xq,cq,nq) c c E=ALPHA*sum(j=1,nat)(potj-pot0j)**2 write(3,3004) 3004 format(' atom current potential desired barrier') do 1 j=1,nat write(3,3005)j,potc(j),pot(j),ALPHA 3005 format(i5,3f15.6) 1 e=e+ALPHA*(potc(j)-pot(j))**2 c c Ex=2*ALPHA*sum(j=1,nat)(potj-potj0)*potjx (potjx approximate) c actually Forces, i.e. negative gradient taken do 3 l=1,nat il=3*(l-1) pfac=2.0d0*ALPHA*(potc(l)-pot(l))*bohr xa=x(il+1) ya=x(il+2) za=x(il+3) do 3 i=1,nq it=3*i-3 d2=(xa-xq(it+1))**2+(ya-xq(it+2))**2+(za-xq(it+3))**2 d=dsqrt(d2) pp=pfac*cq(i)/d/d2 do 3 lx=1,3 ii=il+lx 3 gg(ii)=gg(ii)+pp*(x(ii)-xq(it+lx)) c write(3,3002)e,eold,e-eold 3002 format(' Etot Emol Ees: ',3f17.6,/, 1'Charge Force: (hartree/bohr)') do 62 i=1,nat 62 write(3,3003)i,(gg(3*(i-1)+j),j=1,3) 3003 format(i4,3f12.6) return end c ===================================== subroutine fixff(npar,fx,ipar,kpar,npar0,gg,gx, 1bpar,apar,NAT0,x,nat,p,fxpar,pp,nsame,ijsame,bijsame, 2p1,p2,pp1,pp2) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ipar(npar0,4),kpar(npar0),bpar(npar0), 1apar(npar0),x(3*NAT0),pp(4,3,4,3),gg(3*NAT0),gx(3*NAT0), 2p(4,3),fxpar(3*NAT0,3*NAT0),fx(3*NAT0,3*NAT0), 3ijsame(npar0,2),bijsame(npar0,2), 2p1(4,3),p2(4,3),pp1(4,3,4,3),pp2(4,3,4,3) call msg('Adding fixed parameter terms for FF') do 2 i=1,3*nat do 2 j=1,3*nat 2 fxpar(i,j)=0.0d0 c do 1 i=1,npar kind=kpar(i) a0=apar(i) B=bpar(i) call getpar(NAT0,x,kind,p,actual,npar0,ipar,i,pp) pf=2.0d0*b*(actual-a0) of=2.0d0*b c c Exy=2B(actual-a0) actualxy + 2 B actualx actualy do 1 ia=1,kind do 1 ix=1,3 ii=3*(ipar(i,ia)-1)+ix do 1 ja=1,kind do 1 iy=1,3 jj=3*(ipar(i,ja)-1)+iy 1 fxpar(ii,jj)=fxpar(ii,jj)+pf*pp(ia,ix,ja,iy)+of*p(ia,ix)*p(ja,iy) c c Same constants: do 4 i=1,nsame b=bijsame(i,1) d=bijsame(i,2) i1=ijsame(i,1) i2=ijsame(i,2) kind1=kpar(i1) kind2=kpar(i2) call getpar(NAT0,x,kind1,p1,actual1,npar0,ipar,i1,pp1) call getpar(NAT0,x,kind2,p2,actual2,npar0,ipar,i2,pp2) pf=2.0d0*b*(actual1-actual2-d) of=2.0d0*b c c Exy=2B(actual1-actual2-d)(actual1xy-actual2xy) c + 2 B (actual1x - actual2x)(actual1y - actual2y) do 4 iat=1,nat i11=0 do 41 ia=1,kind1 if(ipar(i1,ia).eq.iat)then i11=ia goto 411 endif 41 continue 411 i12=0 do 42 ja=1,kind2 if(ipar(i2,ja).eq.iat)then i12=ja goto 412 endif 42 continue 412 continue do 4 jat=1,nat i21=0 do 43 ia=1,kind1 if(ipar(i1,ia).eq.jat)then i21=ia goto 413 endif 43 continue 413 i22=0 do 44 ja=1,kind2 if(ipar(i2,ja).eq.jat)then i22=ja goto 414 endif 44 continue 414 continue do 4 ix=1,3 ii=3*(iat-1)+ix do 4 iy=1,3 jj=3*(jat-1)+iy c c 2nd derivative of first parameter pp1s=0.0d0 if(i11.gt.0.and.i21.gt.0)pp1s=pp1(i11,ix,i21,iy) c c 2nd derivative of second parameter pp2s=0.0d0 if(i12.gt.0.and.i22.gt.0)pp2s=pp1(i12,ix,i22,iy) c c 1st derivatives of first parameter p1x=0.0d0 p1y=0.0d0 if(i11.gt.0)p1x=p1(i11,ix) if(i21.gt.0)p1y=p1(i21,iy) c c 1st derivatives of second parameter p2x=0.0d0 p2y=0.0d0 if(i12.gt.0)p2x=p2(i12,ix) if(i22.gt.0)p2y=p2(i22,iy) 4 fxpar(ii,jj)=fxpar(ii,jj)+pf*(pp1s-pp2s)+of*(p1x-p2x)*(p1y-p2y) c call msg('Parameter FF:') c do 4 i=1,3*nat c write(3,*)i c4 write(3,309)(fxpar(i,j),j=1,3*nat) c309 format(8f10.4) bohr=0.529177d0 do 3 i=1,3*nat gx(i)=gx(i)+gg(i) do 3 j=1,3*nat 3 fx(i,j)=fx(i,j)+fxpar(i,j)*bohr**2 call msg('Parameter FF and gradient added.') write(3,3002) 3002 format(' Total gradient (Hartree/bohr):') do 62 i=1,nat 62 write(3,3003)i,(gx(3*(i-1)+j),j=1,3) 3003 format(i4,3f12.6) return end c ===================================== subroutine fixffc(nat,fxpar,NAT0,x,gg,gx,fx,potc,pot,xq,cq,nq) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension x(*),gg(*),gx(*),fxpar(3*NAT0,3*NAT0),fx(3*NAT0,3*NAT0), 1pot(*),potc(*),xq(*),cq(*) c bohr=0.529177d0 call msg('Adding charge terms for FF') do 9 i=1,3*nat do 9 j=1,3*nat 9 fxpar(i,j)=0.0d0 call loadbar(1.0d0,ALPHA,'ALPHA') c c current potential call gcp(nat,potc,x,xq,cq,nq) c do 1 l=1,nat il=3*(l-1) xa=x(il+1) ya=x(il+2) za=x(il+3) s0=0.0d0 sx=0.0d0 sy=0.0d0 sz=0.0d0 sxx=0.0d0 sxy=0.0d0 sxz=0.0d0 syy=0.0d0 syz=0.0d0 szz=0.0d0 do 11 i=1,nq it=3*(i-1) xi=xq(it+1) yi=xq(it+2) zi=xq(it+3) qi=cq(i) d2=(xa-xi)**2+(ya-yi)**2+(za-zi)**2 d=dsqrt(d2) d3=d*d2 d5=d3*d2 s0=s0-qi/d3 sx=sx-qi/d3*(xa-xi) sy=sy-qi/d3*(ya-yi) sz=sz-qi/d3*(za-zi) sxx=sxx+3.0d0*qi/d5*(xa-xi)*(xa-xi) sxy=sxy+3.0d0*qi/d5*(xa-xi)*(ya-yi) sxz=sxz+3.0d0*qi/d5*(xa-xi)*(za-zi) syy=syy+3.0d0*qi/d5*(ya-yi)*(ya-yi) syz=syz+3.0d0*qi/d5*(ya-yi)*(za-zi) 11 szz=szz+3.0d0*qi/d5*(za-zi)*(za-zi) p2=2.0d0*ALPHA*bohr**2 p1=p2*(potc(l)-pot(l)) fxpar(il+1,il+1)=p1*(s0+sxx)+p2*sx*sx fxpar(il+1,il+2)=p1*( +sxy)+p2*sx*sy fxpar(il+1,il+3)=p1*( +sxz)+p2*sx*sz fxpar(il+2,il+1)=p1*( +sxy)+p2*sy*sx fxpar(il+2,il+2)=p1*(s0+syy)+p2*sy*sy fxpar(il+2,il+3)=p1*( +syz)+p2*sy*sz fxpar(il+3,il+1)=p1*( +sxz)+p2*sz*sx fxpar(il+3,il+2)=p1*( +syz)+p2*sz*sy 1 fxpar(il+3,il+3)=p1*(s0+szz)+p2*sz*sz c call msg('Charge FF (hartree/bohr^2):') c do 4 i=1,3*nat c write(3,*)i c4 write(3,309)(fxpar(i,j),j=1,3*nat) c309 format(8f10.4) open(11,file='AG.SCR') open(12,file='AGPOT.SCR') do 31 i=1,nat write(11,1111)(gx(3*(i-1)+j),j=1,3) 31 write(12,1111)(gg(3*(i-1)+j),j=1,3) 1111 format(3f12.6) close(11) close(12) do 3 i=1,3*nat gx(i)=gx(i)+gg(i) do 3 j=1,3*nat 3 fx(i,j)=fx(i,j)+fxpar(i,j) call msg('Charge FF and gradient added.') c call WRITEFF(NAT0,nat,fx,'all.fc') write(3,3002) 3002 format(' Total gradient:') do 62 i=1,nat 62 write(3,3003)i,(gx(3*(i-1)+j),j=1,3) 3003 format(i4,3f12.6) return end c ===================================== subroutine gcp(nat,potc,x,xq,cq,nq) integer*4 i,nat,nq,ic real*8 potc(*),xa,ya,za,x(*),d,xq(*),cq(*) do 2 i=1,nat potc(i)=0.0d0 xa=x(3*(i-1)+1) ya=x(3*(i-1)+2) za=x(3*(i-1)+3) do 2 ic=1,nq d=dsqrt((xa-xq(3*(ic-1)+1))**2 1 +(ya-xq(3*(ic-1)+2))**2 2 +(za-xq(3*(ic-1)+3))**2) 2 potc(i)=potc(i)+cq(ic)/d return end c ===================================== subroutine genchar(NAT,NC,pot,x,xq,cq) c c NC - number of charges, NC.GE.NAT+2 c NAT .. number of atoms c M .. matrix dimension c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) logical lpref parameter (NAT0=900,M0=1000,nc00=1000) dimension pot(*),x(*),xq(*),a(M0,M0),a0(M0,M0),ipiv(M0), 1work(3*M0*M0),b(M0),c(M0),cq(*) c M=NC+NAT+1 if(NC.lt.NAT)call bye('NC < NAT !') if(NC.gt.nc00)call bye('NC > nc00 !') if(NAT.gt.nat0)call bye('nat0 < NAT !') if(M.gt.M0)call bye('M > M0!') c c set up the matrix do 14 i=1,M do 14 j=1,M 14 a(i,j)=0.0d0 do 15 i=1,NC a(i,i)=1.0d0 a(i,M)=1.0d0 15 a(M,i)=1.0d0 do 16 ia=1,NAT xa=x(3*(ia-1)+1) ya=x(3*(ia-1)+2) za=x(3*(ia-1)+3) do 16 ic=1,NC xi=xq(3*(ic-1)+1) yi=xq(3*(ic-1)+2) zi=xq(3*(ic-1)+3) rij=dsqrt((xi-xa)**2+(yi-ya)**2+(zi-za)**2) a(NC+ia,ic)=1.0d0/rij 16 a(ic,NC+ia)=1.0d0/rij do 17 i=1,M do 17 j=1,M 17 a0(i,j)=a(i,j) M1=M M2=M call dgetrf(M1,M2,a,M0,ipiv,info) if(info.eq.0)then write(3,*)' Pivoting OK' else call bye('pivoting failed') endif c lwork=3*M*M call dgetri(M,a,M0,ipiv,work,lwork,info) c if(info.eq.0)write(3,*)'successful inversion' if(info.lt.0)then write(3,*)info call bye('invalid argument') endif c tol=1.0d-9 nmo=M do i=1,nmo do j=1,nmo sum=0.0d0 sum2=0.0d0 do ii=1,nmo sum=sum+a(i,ii)*a0(ii,j) sum2=sum2+a(ii,i)*a0(j,ii) enddo ic=0 if(abs(sum ).gt.tol.and.i.ne.j)ic=1 if(abs(sum2).gt.tol.and.i.ne.j)ic=2 if(abs(sum -1.0d0).gt.tol.and.i.eq.j)ic=3 if(abs(sum2-1.0d0).gt.tol.and.i.eq.j)ic=4 if(ic.ne.0)then write(3,*)'ic = ',ic write(3,*)i,j,sum,sum2 write(3,*)i,j,sum call bye('Inversion error') endif enddo enddo c write(3,*)'Inversion Check OK' c do 4 i=1,M 4 c(i)=0.0d0 inquire(file='OLD.CHA',exist=lpref) if(lpref)then write(3,*)'Prefered charges found in OLD.CHA' open(31,file='OLD.CHA') read(31,*) read(31,*) do 61 i=1,NC 61 read(31,*) c(i),c(i) close(31) else write(3,*)'Minimal charges seeked.' endif do 6 i=1,NAT 6 c(nc+i)=pot(i) c do 18 i=1,M sum=0.0d0 do 19 j=1,M 19 sum=sum+a(i,j)*c(j) 18 b(i)=sum write(3,88) 88 format(' charge old new') do 20 i=1,NC 20 write(3,89)i,cq(i),b(i) 89 format(i6,2f12.4) do 21 i=1,NC 21 cq(i)=b(i) return end c c =========================================== c LAPACK ROUTINES: c =========================================== SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. INTEGER*4 INFO, LDA, LWORK, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * * Purpose * ======= * * DGETRI computes the inverse of a matrix using the LU factorization * computed by DGETRF. * * This method inverts U and then computes inv(A) by solving the system * inv(A)*L = inv(U) for inv(A). * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the factors L and U from the factorization * A = P*L*U as computed by DGETRF. * On exit, if INFO = 0, the inverse of the original matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (input) INTEGER*4 array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) * On exit, if INFO=0, then WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * For optimal performance LWORK >= N*NB, where NB is * the optimal blocksize returned by ILAENV. * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero; the matrix is * singular and its inverse could not be computed. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER*4 I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB, $ NBMIN, NN * .. * .. External Functions .. INTEGER*4 ILAENV EXTERNAL ILAENV * .. * .. External Subroutines .. EXTERNAL DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 ) LWKOPT = N*NB WORK( 1 ) = LWKOPT LQUERY = ( LWORK.EQ.-1 ) IF( N.LT.0 ) THEN INFO = -1 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -3 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRI', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Form inv(U). If INFO > 0 from DTRTRI, then U is singular, * and the inverse is not computed. * CALL DTRTRI( 'Upper', 'N', N, A, LDA, INFO ) IF( INFO.GT.0 ) $ RETURN * NBMIN = 2 LDWORK = N IF( NB.GT.1 .AND. NB.LT.N ) THEN IWS = MAX( LDWORK*NB, 1 ) IF( LWORK.LT.IWS ) THEN NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) ) END IF ELSE IWS = N END IF * * Solve the equation inv(A)*L = inv(U) for inv(A). * IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN * * Use unblocked code. * DO 20 J = N, 1, -1 * * Copy current column of L to WORK and replace with zeros. * DO 10 I = J + 1, N WORK( I ) = A( I, J ) A( I, J ) = ZERO 10 CONTINUE * * Compute current column of inv(A). * IF( J.LT.N ) $ CALL DGEMV( 'N', N, N-J, -ONE, A( 1, J+1 ), $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) 20 CONTINUE ELSE * * Use blocked code. * NN = ( ( N-1 ) / NB )*NB + 1 DO 50 J = NN, 1, -NB JB = MIN( NB, N-J+1 ) * * Copy current block column of L to WORK and replace with * zeros. * DO 40 JJ = J, J + JB - 1 DO 30 I = JJ + 1, N WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) A( I, JJ ) = ZERO 30 CONTINUE 40 CONTINUE * * Compute current block column of inv(A). * IF( J+JB.LE.N ) $ CALL DGEMM( 'N', 'N', N, JB, $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA, $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA ) CALL DTRSM( 'R', 'L', 'N', 'U', N, JB, $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) 50 CONTINUE END IF * * Apply column interchanges. * DO 60 J = N - 1, 1, -1 JP = IPIV( J ) IF( JP.NE.J ) $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 ) 60 CONTINUE * WORK( 1 ) = IWS RETURN * * End of DGETRI * END SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER*4 INFO, LDA, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DTRTI2 computes the inverse of arreal upper or lower triangular * matrix. * * This is the Level 2 BLAS version of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the matrix A is upper or lower triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * DIAG (input) CHARACTER*1 * Specifies whether or not the matrix A is unit triangular. * = 'N': Non-unit triangular * = 'U': Unit triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the triangular matrix A. If UPLO = 'U', the * leading n by n upper triangular part of the array A contains * the upper triangular matrix, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading n by n lower triangular part of the array A contains * the lower triangular matrix, and the strictly upper * triangular part of A is not referenced. If DIAG = 'U', the * diagonal elements of A are also not referenced and are * assumed to be 1. * * On exit, the (triangular) inverse of the original matrix, in * the same storage format. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER*4 J DOUBLE PRECISION AJJ * .. * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. * .. External Subroutines .. EXTERNAL DSCAL, DTRMV, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LBAME( UPLO, 'U' ) NOUNIT = LBAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LBAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LBAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DTRTI2', -INFO ) RETURN END IF * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix. * DO 10 J = 1, N IF( NOUNIT ) THEN A( J, J ) = ONE / A( J, J ) AJJ = -A( J, J ) ELSE AJJ = -ONE END IF * * Compute elements 1:j-1 of j-th column. * CALL DTRMV( 'Upper', 'N', DIAG, J-1, A, LDA, $ A( 1, J ), 1 ) CALL DSCAL( J-1, AJJ, A( 1, J ), 1 ) 10 CONTINUE ELSE * * Compute inverse of lower triangular matrix. * DO 20 J = N, 1, -1 IF( NOUNIT ) THEN A( J, J ) = ONE / A( J, J ) AJJ = -A( J, J ) ELSE AJJ = -ONE END IF IF( J.LT.N ) THEN * * Compute elements j+1:n of j-th column. * CALL DTRMV( 'L', 'N', DIAG, N-J, $ A( J+1, J+1 ), LDA, A( J+1, J ), 1 ) CALL DSCAL( N-J, AJJ, A( J+1, J ), 1 ) END IF 20 CONTINUE END IF * RETURN * * End of DTRTI2 * END SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER*4 INFO, LDA, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DTRTRI computes the inverse of arreal upper or lower triangular * matrix A. * * This is the Level 3 BLAS version of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * = 'U': A is upper triangular; * = 'L': A is lower triangular. * * DIAG (input) CHARACTER*1 * = 'N': A is non-unit triangular; * = 'U': A is unit triangular. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the triangular matrix A. If UPLO = 'U', the * leading N-by-N upper triangular part of the array A contains * the upper triangular matrix, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading N-by-N lower triangular part of the array A contains * the lower triangular matrix, and the strictly upper * triangular part of A is not referenced. If DIAG = 'U', the * diagonal elements of A are also not referenced and are * assumed to be 1. * On exit, the (triangular) inverse of the original matrix, in * the same storage format. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, A(i,i) is exactly zero. The triangular * matrix is singular and its inverse can not be computed. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER*4 J, JB, NB, NN * .. * .. External Functions .. LOGICAL LBAME INTEGER*4 ILAENV EXTERNAL LBAME, ILAENV * .. * .. External Subroutines .. EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LBAME( UPLO, 'U' ) NOUNIT = LBAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LBAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LBAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DTRTRI', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Check for singularity if non-unit. * IF( NOUNIT ) THEN DO 10 INFO = 1, N IF( A( INFO, INFO ).EQ.ZERO ) $ RETURN 10 CONTINUE INFO = 0 END IF * * Determine the block size for this environment. * NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.N ) THEN * * Use unblocked code * CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) ELSE * * Use blocked code * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix * DO 20 J = 1, N, NB JB = MIN( NB, N-J+1 ) * * Compute rows 1:j-1 of current block column * CALL DTRMM( 'L', 'Upper', 'N', DIAG, J-1, $ JB, ONE, A, LDA, A( 1, J ), LDA ) CALL DTRSM( 'R', 'Upper', 'N', DIAG, J-1, $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA ) * * Compute inverse of current diagonal block * CALL DTRTI2( 'U', DIAG, JB, A( J, J ), LDA, INFO ) 20 CONTINUE ELSE * * Compute inverse of lower triangular matrix * NN = ( ( N-1 ) / NB )*NB + 1 DO 30 J = NN, 1, -NB JB = MIN( NB, N-J+1 ) IF( J+JB.LE.N ) THEN * * Compute rows j+jb:n of current block column * CALL DTRMM( 'L', 'L', 'N', DIAG, $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA, $ A( J+JB, J ), LDA ) CALL DTRSM( 'R', 'L', 'N', DIAG, $ N-J-JB+1, JB, -ONE, A( J, J ), LDA, $ A( J+JB, J ), LDA ) END IF * * Compute inverse of current diagonal block * CALL DTRTI2( 'L', DIAG, JB, A( J, J ), LDA, INFO ) 30 CONTINUE END IF END IF * RETURN * * End of DTRTRI * END LOGICAL FUNCTION LBAME( CA, CB ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LBAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER*4 INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LBAME = CA.EQ.CB IF( LBAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LBAME = INTA.EQ.INTB * * RETURN * * End of LBAME * END SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER*4 INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * ===================================================================== * * .. Executable Statements .. * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER*4 M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER*4 I, INFO, J, L, NROWA, NROWB DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LBAME( TRANSA, 'N' ) NOTB = LBAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M c NCOLA = K ELSE NROWA = K c NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LBAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LBAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LBAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of DGEMM . * END SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER*4 INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER*4 I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LBAME( TRANS, 'N' ).AND. $ .NOT.LBAME( TRANS, 'T' ).AND. $ .NOT.LBAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LBAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LBAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * END subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real*8 dx(*),dy(*),dtemp integer*4 i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER*4 M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRMM performs one of the matrix-matrix operations * * B := alpha*op( A )*B, or B := alpha*B*op( A ), * * where alpha is a scalar, B is an m by n matrix, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) multiplies B from * the left or right as follows: * * SIDE = 'L' or 'l' B := alpha*op( A )*B. * * SIDE = 'R' or 'r' B := alpha*B*op( A ). * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the matrix B, and on exit is overwritten by the * transformed matrix. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER*4 I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LBAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LBAME( DIAG , 'N' ) UPPER = LBAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LBAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LBAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LBAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LBAME( DIAG , 'U' ) ).AND. $ ( .NOT.LBAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*A*B. * IF( UPPER )THEN DO 50, J = 1, N DO 40, K = 1, M IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) DO 30, I = 1, K - 1 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 30 CONTINUE IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) B( K, J ) = TEMP END IF 40 CONTINUE 50 CONTINUE ELSE DO 80, J = 1, N DO 70 K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) B( K, J ) = TEMP IF( NOUNIT ) $ B( K, J ) = B( K, J )*A( K, K ) DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * Form B := alpha*A'*B. * IF( UPPER )THEN DO 110, J = 1, N DO 100, I = M, 1, -1 TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE B( I, J ) = ALPHA*TEMP 100 CONTINUE 110 CONTINUE ELSE DO 140, J = 1, N DO 130, I = 1, M TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 120, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 120 CONTINUE B( I, J ) = ALPHA*TEMP 130 CONTINUE 140 CONTINUE END IF END IF ELSE IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*A. * IF( UPPER )THEN DO 180, J = N, 1, -1 TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = 1, M B( I, J ) = TEMP*B( I, J ) 150 CONTINUE DO 170, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 160, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE ELSE DO 220, J = 1, N TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 190, I = 1, M B( I, J ) = TEMP*B( I, J ) 190 CONTINUE DO 210, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 200, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 200 CONTINUE END IF 210 CONTINUE 220 CONTINUE END IF ELSE * * Form B := alpha*B*A'. * IF( UPPER )THEN DO 260, K = 1, N DO 240, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 230, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 230 CONTINUE END IF 240 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 250, I = 1, M B( I, K ) = TEMP*B( I, K ) 250 CONTINUE END IF 260 CONTINUE ELSE DO 300, K = N, 1, -1 DO 280, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 270, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 270 CONTINUE END IF 280 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE END IF 300 CONTINUE END IF END IF END IF * RETURN * * End of DTRMM . * END SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER*4 INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * Purpose * ======= * * DTRMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER*4 I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LBAME( UPLO , 'U' ).AND. $ .NOT.LBAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LBAME( TRANS, 'N' ).AND. $ .NOT.LBAME( TRANS, 'T' ).AND. $ .NOT.LBAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LBAME( DIAG , 'U' ).AND. $ .NOT.LBAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LBAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LBAME( TRANS, 'N' ) )THEN * * Form x := A*x. * IF( LBAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*A( I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*A( I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A'*x. * IF( LBAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = N, 1, -1 TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 110, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 130, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of DTRMV . * END SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER*4 M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER*4 I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LBAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LBAME( DIAG , 'N' ) UPPER = LBAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LBAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LBAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LBAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LBAME( DIAG , 'U' ) ).AND. $ ( .NOT.LBAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of DTRSM . * END c ooooooooooooooo SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. real*8 ALPHA INTEGER*4 INCX, INCY, LDA, M, N * .. Array Arguments .. real*8 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - real*8. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - real*8 array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - real*8 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - real*8 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. real*8 ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. real*8 TEMP INTEGER*4 I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of DGER . * END SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1992 * * .. Scalar Arguments .. INTEGER*4 INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) real*8 A( LDA, * ) * .. * * Purpose * ======= * * DGETF2 computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 2 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) real*8 array, dimension (LDA,N) * On entry, the m by n matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER*4 array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. real*8 ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER*4 J, JP * .. * .. External Functions .. INTEGER*4 IDAMAX EXTERNAL IDAMAX * .. * .. External Subroutines .. EXTERNAL DGER, DSCAL, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETF2', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * DO 10 J = 1, MIN( M, N ) * * Find pivot and test for singularity. * JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Apply the interchange to columns 1:N. * IF( JP.NE.J ) $ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) * ELSE IF( INFO.EQ.0 ) THEN * INFO = J END IF * IF( J.LT.MIN( M, N ) ) THEN * * Update trailing submatrix. * CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, $ A( J+1, J+1 ), LDA ) END IF 10 CONTINUE RETURN * * End of DGETF2 * END SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER*4 INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) real*8 A( LDA, * ) * .. * * Purpose * ======= * * DGETRF computes an LU factorization of a general M-by-N matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 3 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INT * The number of rows of the matrix A. M >= 0. * * N (input) INT * The number of columns of the matrix A. N >= 0. * * A (input/output) real*8 array, dimension (LDA,N) * On entry, the M-by-N matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INT * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER*4 array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. real*8 ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER*4 I, IINFO, J, JB, NB * .. * .. External Subroutines .. EXTERNAL DGEMM, DGETF2, DLASWP, DTRSM, XERBLA * .. * .. External Functions .. INTEGER*4 ILAENV EXTERNAL ILAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRF', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN * * Use unblocked code. * CALL DGETF2( M, N, A, LDA, IPIV, INFO ) ELSE * * Use blocked code. * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Factor diagonal and subdiagonal blocks and test for exact * singularity. * CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) * * Adjust INFO and the pivot indices. * IF( INFO.EQ.0 .AND. IINFO.GT.0 ) $ INFO = IINFO + J - 1 DO 10 I = J, MIN( M, J+JB-1 ) IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE * * Apply interchanges to columns 1:J-1. * CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) * IF( J+JB.LE.N ) THEN * * Apply interchanges to columns J+JB:N. * CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, $ IPIV, 1 ) * * Compute block row of U. * CALL DTRSM( 'L', 'L', 'N', 'U', JB, $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), $ LDA ) IF( J+JB.LE.M ) THEN * * Update trailing submatrix. * CALL DGEMM( 'N', 'N', M-J-JB+1, $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), $ LDA ) END IF END IF 20 CONTINUE END IF RETURN * * End of DGETRF * END SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. INTEGER*4 INCX, K1, K2, LDA, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) real*8 A( LDA, * ) * .. * * Purpose * ======= * * DLASWP performs a series of row interchanges on the matrix A. * One row interchange is initiated for each of rows K1 through K2 of A. * * Arguments * ========= * * N (input) INTEGER * The number of columns of the matrix A. * * A (input/output) real*8 array, dimension (LDA,N) * On entry, the matrix of column dimension N to which the row * interchanges will be applied. * On exit, the permuted matrix. * * LDA (input) INTEGER * The leading dimension of the array A. * * K1 (input) INTEGER * The first element of IPIV for which a row interchange will * be done. * * K2 (input) INTEGER * The last element of IPIV for which a row interchange will * be done. * * IPIV (input) INTEGER*4 array, dimension (M*abs(INCX)) * The vector of pivot indices. Only the elements in positions * K1 through K2 of IPIV are accessed. * IPIV(K) = L implies rows K and L are to be interchanged. * * INCX (input) INTEGER * The increment between successive values of IPIV. If IPIV * is negative, the pivots are applied in reverse order. * * Further Details * =============== * * Modified by * R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA * * ===================================================================== * * .. Local Scalars .. INTEGER*4 I, I1, I2, INC, IP, IX, IX0, J, K, N32 real*8 TEMP * .. * .. Executable Statements .. * * Interchange row I with row IPIV(I) for each of rows K1 through K2. * IF( INCX.GT.0 ) THEN IX0 = K1 I1 = K1 I2 = K2 INC = 1 ELSE IF( INCX.LT.0 ) THEN IX0 = 1 + ( 1-K2 )*INCX I1 = K2 I2 = K1 INC = -1 ELSE RETURN END IF * N32 = ( N / 32 )*32 IF( N32.NE.0 ) THEN DO 30 J = 1, N32, 32 IX = IX0 DO 20 I = I1, I2, INC IP = IPIV( IX ) IF( IP.NE.I ) THEN DO 10 K = J, J + 31 TEMP = A( I, K ) A( I, K ) = A( IP, K ) A( IP, K ) = TEMP 10 CONTINUE END IF IX = IX + INCX 20 CONTINUE 30 CONTINUE END IF IF( N32.NE.N ) THEN N32 = N32 + 1 IX = IX0 DO 50 I = I1, I2, INC IP = IPIV( IX ) IF( IP.NE.I ) THEN DO 40 K = N32, N TEMP = A( I, K ) A( I, K ) = A( IP, K ) A( IP, K ) = TEMP 40 CONTINUE END IF IX = IX + INCX 50 CONTINUE END IF * RETURN * * End of DLASWP * END subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real*8 da,dx(*) integer*4 i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end integer*4 function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real*8 dx(*),dmax integer*4 i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end INTEGER*4 FUNCTION IEEECK( ISPEC, ZERO, ONE ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1998 * * .. Scalar Arguments .. INTEGER*4 ISPEC REAL*8 ONE, ZERO * .. * * Purpose * ======= * * IEEECK is called from the ILAENV to verify that Infinity and * possibly NaN arithmetic is safe (i.e. will not trap). * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies whether to test just for inifinity arithmetic * or whether to test for infinity and NaN arithmetic. * = 0: Verify infinity arithmetic only. * = 1: Verify infinity and NaN arithmetic. * * ZERO (input) REAL * Must contain the value 0.0 * This is passed to prevent the compiler from optimizing * away this code. * * ONE (input) REAL * Must contain the value 1.0 * This is passed to prevent the compiler from optimizing * away this code. * * RETURN VALUE: INTEGER * = 0: Arithmetic failed to produce the correct answers * = 1: Arithmetic produced the correct answers * * .. Local Scalars .. REAL*8 NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGINF, $ NEGZRO, NEWZRO, POSINF * .. * .. Executable Statements .. IEEECK = 1 * POSINF = ONE / ZERO IF( POSINF.LE.ONE ) THEN IEEECK = 0 RETURN END IF * NEGINF = -ONE / ZERO IF( NEGINF.GE.ZERO ) THEN IEEECK = 0 RETURN END IF * NEGZRO = ONE / ( NEGINF+ONE ) IF( NEGZRO.NE.ZERO ) THEN IEEECK = 0 RETURN END IF * NEGINF = ONE / NEGZRO IF( NEGINF.GE.ZERO ) THEN IEEECK = 0 RETURN END IF * NEWZRO = NEGZRO + ZERO IF( NEWZRO.NE.ZERO ) THEN IEEECK = 0 RETURN END IF * POSINF = ONE / NEWZRO IF( POSINF.LE.ONE ) THEN IEEECK = 0 RETURN END IF * NEGINF = NEGINF*POSINF IF( NEGINF.GE.ZERO ) THEN IEEECK = 0 RETURN END IF * POSINF = POSINF*POSINF IF( POSINF.LE.ONE ) THEN IEEECK = 0 RETURN END IF * * * * * Return if we were only asked to check infinity arithmetic * IF( ISPEC.EQ.0 ) $ RETURN * NAN1 = POSINF + NEGINF * NAN2 = POSINF / NEGINF * NAN3 = POSINF / POSINF * NAN4 = POSINF*ZERO * NAN5 = NEGINF*NEGZRO * NAN6 = NAN5*0.0d0 * IF( NAN1.EQ.NAN1 ) THEN IEEECK = 0 RETURN END IF * IF( NAN2.EQ.NAN2 ) THEN IEEECK = 0 RETURN END IF * IF( NAN3.EQ.NAN3 ) THEN IEEECK = 0 RETURN END IF * IF( NAN4.EQ.NAN4 ) THEN IEEECK = 0 RETURN END IF * IF( NAN5.EQ.NAN5 ) THEN IEEECK = 0 RETURN END IF * IF( NAN6.EQ.NAN6 ) THEN IEEECK = 0 RETURN END IF * RETURN END INTEGER*4 FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, $ N4 ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER*4 ISPEC, N1, N2, N3, N4 * .. * * Purpose * ======= * * ILAENV is called from the LAPACK routines to choose problem-dependent * parameters for the local environment. See ISPEC for a description of * the parameters. * * This version provides a set of parameters which should give good, * but not optimal, performance on many of the currently available * computers. Users are encouraged to modify this subroutine to set * the tuning parameters for their particular machine using the option * and problem size information in the arguments. * * This routine will not function correctly if it is converted to all * lower case. Converting it to all upper case is allowed. * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies the parameter to be returned as the value of * ILAENV. * = 1: the optimal blocksize; if this value is 1, an unblocked * algorithm will give the best performance. * = 2: the minimum block size for which the block routine * should be used; if the usable block size is less than * this value, an unblocked routine should be used. * = 3: the crossover point (in a block routine, for N less * than this value, an unblocked routine should be used) * = 4: the number of shifts, used in the nonsymmetric * eigenvalue routines * = 5: the minimum column dimension for blocking to be used; * rectangular blocks must have dimension at least k by m, * where k is given by ILAENV(2,...) and m by ILAENV(5,...) * = 6: the crossover point for the SVD (when reducing an m by n * matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds * this value, a QR factorization is used first to reduce * the matrix to a triangular form.) * = 7: the number of processors * = 8: the crossover point for the multishift QR and QZ methods * for nonsymmetric eigenvalue problems. * = 9: maximum size of the subproblems at the bottom of the * computation tree in the divide-and-conquer algorithm * (used by xGELSD and xGESDD) * =10: ieee NaN arithmetic can be trusted not to trap * =11: infinity arithmetic can be trusted not to trap * * NAME (input) CHARACTER*(*) * The name of the calling subroutine, in either upper case or * lower case. * * OPTS (input) CHARACTER*(*) * The character options to the subroutine NAME, concatenated * into a single character string. For example, UPLO = 'U', * TRANS = 'T', and DIAG = 'N' for a triangular routine would * be specified as OPTS = 'UTN'. * * N1 (input) INTEGER * N2 (input) INTEGER * N3 (input) INTEGER * N4 (input) INTEGER * Problem dimensions for the subroutine NAME; these may not all * be required. * * (ILAENV) (output) INTEGER * >= 0: the value of the parameter specified by ISPEC * < 0: if ILAENV = -k, the k-th argument had an illegal value. * * Further Details * =============== * * The following conventions have been used when calling ILAENV from the * LAPACK routines: * 1) OPTS is a concatenation of all of the character options to * subroutine NAME, in the same order that they appear in the * argument list for NAME, even if they are not used in determining * the value of the parameter specified by ISPEC. * 2) The problem dimensions N1, N2, N3, N4 are specified in the order * that they appear in the argument list for NAME. N1 is used * first, N2 second, and so on, and unused problem dimensions are * passed a value of -1. * 3) The parameter value returned by ILAENV is checked for validity in * the calling subroutine. For example, ILAENV is used to retrieve * the optimal blocksize for STRTRI as follows: * * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) * IF( NB.LE.1 ) NB = MAX( 1, N ) * * ===================================================================== * * .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER*4 I, IC, IZ, NB, NBMIN, NX * .. * .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL * .. * .. External Functions .. INTEGER*4 IEEECK EXTERNAL IEEECK * .. * .. Executable Statements .. * GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000, $ 1100 ) ISPEC * * Invalid value for ISPEC * ILAENV = -1 RETURN * 100 CONTINUE * * Convert NAME to upper case if the first character is lower case. * ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN * * ASCII character set * IF( IC.GE.97 .AND. IC.LE.122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF * ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN * * EBCDIC character set * IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) $ SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF * ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN * * Prime machines: ASCII+128 * IF( IC.GE.225 .AND. IC.LE.250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF * C1 = SUBNAM( 1:1 ) SNAME = C1.EQ.'S' .OR. C1.EQ.'D' CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' IF( .NOT.( CNAME .OR. SNAME ) ) $ RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) * GO TO ( 110, 200, 300 ) ISPEC * 110 CONTINUE * * ISPEC = 1: block size * * In these examples, separate code is provided for setting NB for * real and complex. We assume that NB will take the same value in * single or real*8. * NB = 1 * IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'PO' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NB = 32 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRF' ) THEN NB = 64 ELSE IF( C3.EQ.'TRD' ) THEN NB = 32 ELSE IF( C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( C2.EQ.'GB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'PB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'TR' ) THEN IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN IF( C3.EQ.'EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN * 200 CONTINUE * * ISPEC = 2: minimum block size * NBMIN = 2 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NBMIN = 8 ELSE NBMIN = 8 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN * 300 CONTINUE * * ISPEC = 3: crossover point * NX = 0 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( SNAME .AND. C3.EQ.'TRD' ) THEN NX = 32 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NX = 32 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN * 400 CONTINUE * * ISPEC = 4: number of shifts (used by xHSEQR) * ILAENV = 6 RETURN * 500 CONTINUE * * ISPEC = 5: minimum column dimension (not used) * ILAENV = 2 RETURN * 600 CONTINUE * * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) * ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN * 700 CONTINUE * * ISPEC = 7: number of processors (not used) * ILAENV = 1 RETURN * 800 CONTINUE * * ISPEC = 8: crossover point for multishift (used by xHSEQR) * ILAENV = 50 RETURN * 900 CONTINUE * * ISPEC = 9: maximum size of the subproblems at the bottom of the * computation tree in the divide-and-conquer algorithm * (used by xGELSD and xGESDD) * ILAENV = 25 RETURN * 1000 CONTINUE * * ISPEC = 10: ieee NaN arithmetic can be trusted not to trap * C ILAENV = 0 ILAENV = 1 IF( ILAENV.EQ.1 ) THEN ILAENV = IEEECK( 0, 0.0D0, 1.0D0 ) END IF RETURN * 1100 CONTINUE * * ISPEC = 11: infinity arithmetic can be trusted not to trap * C ILAENV = 0 ILAENV = 1 IF( ILAENV.EQ.1 ) THEN ILAENV = IEEECK( 1, 0.0D0, 1.0D0 ) END IF RETURN * * End of ILAENV * END c ===================================== subroutine filla(nsy0,a2,sl) implicit none integer*4 nsy0,i,j,k real*8 a2(nsy0,3,3) CHARACTER*3 sl(nsy0) sl(1)='c1 ' sl(2)='c2x' sl(3)='c2y' sl(4)='c2z' sl(5)='sxy' sl(6)='syz' sl(7)='szx' sl(8)='c4x' sl(9)='c4y' sl(10)='c4z' sl(11)='m4x' sl(12)='m4y' sl(13)='m4z' do 4 i=1,nsy0 do 4 j=1,3 do 41 k=1,3 41 a2(i,j,k)=0.0d0 4 a2(i,j,j)=1.0d0 c1 ...ok c2x a2(2,2,2)=-1.0d0 a2(2,3,3)=-1.0d0 c2y a2(3,1,1)=-1.0d0 a2(3,3,3)=-1.0d0 c2z a2(4,1,1)=-1.0d0 a2(4,2,2)=-1.0d0 csxy a2(5,3,3)=-1.0d0 csyz a2(6,1,1)=-1.0d0 cszx a2(7,2,2)=-1.0d0 c4x C4x axis, rotation by 90deg a2(8,2,2)= 0.0d0 a2(8,3,3)= 0.0d0 a2(8,2,3)= 1.0d0 a2(8,3,2)=-1.0d0 c4y a2(9,1,1)= 0.0d0 a2(9,3,3)= 0.0d0 a2(9,1,3)= 1.0d0 a2(9,3,1)=-1.0d0 c4z a2(10,1,1)= 0.0d0 a2(10,2,2)= 0.0d0 a2(10,2,1)= 1.0d0 a2(10,1,2)=-1.0d0 cm4x C4x axis, rotation by 90deg, opposite a2(11,2,2)= 0.0d0 a2(11,3,3)= 0.0d0 a2(11,2,3)=-1.0d0 a2(11,3,2)= 1.0d0 cm4y a2(12,1,1)= 0.0d0 a2(12,3,3)= 0.0d0 a2(12,1,3)=-1.0d0 a2(12,3,1)= 1.0d0 cm4z a2(13,1,1)= 0.0d0 a2(13,2,2)= 0.0d0 a2(13,2,1)=-1.0d0 a2(13,1,2)= 1.0d0 return end c =========================================== function ax(is,x,a2,nsy0,ia,ix) implicit none integer*4 is,nsy0,ix,ia real*8 x(*),a2(nsy0,3,3),ax ax=a2(is,ix,1)*x(3*(ia-1)+1) 1 +a2(is,ix,2)*x(3*(ia-1)+2) 2 +a2(is,ix,3)*x(3*(ia-1)+3) return end c ====================================================================== subroutine ginproot(iroot) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) CHARACTER*80 filename,ts c filename='G98.INP' open(10,file='Q.OPT') 1 read(10,1000,end=999,err=999)ts 1000 format(a80) if(ts(1:14).eq.'gaussian input')read(10,1000)filename goto 1 999 close(10) c open(10,file=filename) open(11,file='INP.NEW') 2 read(10,1000,end=899,err=899)ts call rootchek(ts,iroot) write(11,1000)ts goto 2 899 close(10) close(11) call bye('New input written in INP.NEW') end c ====================================================================== subroutine rootchek(ts,iroot) IMPLICIT none CHARACTER*(*) ts CHARACTER*5 s CHARACTER*1 s1 CHARACTER*2 s2 integer*4 iroot,ii,kk,i,ia,j logical lnumber if(iroot.gt.0.and.iroot.lt.10)then write(s1,'(i1)')iroot kk=1 else if(iroot.gt.9.and.iroot.lt.100)then write(s2,'(i2)')iroot kk=2 else write(6,*)iroot call bye('unusual root request') endif endif c do 1 i=1,len(ts)-4 do 2 ii=1,5 ia=ichar(ts(i+ii-1:i+ii-1)) if(ia.ge.65.and.ia.le.90)ia=ia+32 2 s(ii:ii)=char(ia) if(s.eq.'root=')then c where old number ends: do 3 j=i+6,len(ts) 3 if(.not.lnumber(ts(j:j)))goto 4 4 continue c add new number: if(kk.eq.1)ts=ts(1:i+4)//s1//ts(j:len(ts)) if(kk.eq.2)ts=ts(1:i+4)//s2//ts(j:len(ts)) endif 1 continue return end c ====================================================================== function lnumber(c) logical lnumber character*1 c lnumber=c.eq.'0'.or.c.eq.'1'.or.c.eq.'2'.or.c.eq.'3'.or.c.eq.'4' 1 .or.c.eq.'5'.or.c.eq.'6'.or.c.eq.'7'.or.c.eq.'8'.or.c.eq.'9' return end c ====================================================================== subroutine rwfollow(ic,inow,aci,bci,cci,tci) implicit none integer*4 ic,inow,aci(*),bci(*),i real*8 cci(*) character*1 tci(*) open(10,file='FOLLOW.SCR') if(ic.eq.0)then write(10,*)inow do 1 i=1,inow 1 write(10,200)aci(i),bci(i),tci(i),cci(i) 200 format(2i5,a1,f12.5) else read(10,*)inow do 2 i=1,inow 2 read(10,200)aci(i),bci(i),tci(i),cci(i) endif close(10) return end c ====================================================================== subroutine readcicoef(io,inow,imax,a,b,c,t) implicit none integer*4 inow,imax,io,a(*),b(*),oa,ob,j,i1,i2 real*8 c(*),oc logical lnumber character*80 s80 character*1 t(*),aa inow=0 111 read(io,2000,end=299,err=299)s80 2000 format(a80) if(s80(10:11).eq.'->'.or.s80(10:11).eq.'<-')then do 1 i1=12,len(s80) 1 if(lnumber(s80(i1:i1)))goto 2 2 do 3 i2=i1+1,len(s80) 3 if(.not.lnumber(s80(i2:i2)))goto 4 4 i2=i2-1 if(s80(8:8).eq.'A'.or.s80(8:8).eq.'B')then c open shell read(s80( 1: 7),*)oa read(s80(i1:i2),*)ob aa=s80(8:8) else c closed shell read(s80( 1: 8),*,err=993)oa goto 992 993 write(6,2000)s80 call bye('cannot read orbital number') 992 read(s80(i1:i2),*)ob aa='o' endif read(s80(18:30),*)oc if(s80(10:11).eq.'->')then c straight coefficients if(inow.lt.imax)then c if enough space, just add them: inow=inow+1 a(inow)=oa b(inow)=ob c(inow)=oc t(inow)=aa else c if not space, take the biggest: do 209 j=1,imax if(dabs(c(j)).lt.dabs(oc))then a(j)=oa b(j)=ob c(j)=oc t(j)=aa goto 111 endif 209 continue endif else c back coefficients - ignore endif goto 111 endif 299 continue backspace io end