subroutine cctn_tinker1(NAT,atomic,x,y,z,n12,i12,maxval,maxatm) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind,atomic(*),n12(*),i12(maxval,maxatm) PARAMETER (MX=50,MXB=40,IOX=125) c MX maximum number of atoms c MXB maximum number of overlaped atoms for each constant c IOX max # of overlaps common/cctnrest/FF(3*MX,3*MX),R(3,MX),IND(MX),IB(MX,MXB),ND(MX), 1NDS(MX),RS(3,MX),FFS(3*MX,3*MX),IS(MX,MXB),AM(MX,3,3), 2INDBIG(IOX,MX),IBONDSB(4,MX), 3NBO(MX),NBOB(MX),PS(MX,3,3),AS(MX,3,3),PB(MX,3,3),AB(MX,3,3), 4ALPHA(3*MX,3,3),ALPHAS(3*MX,3,3),AA(3*MX,3,3,3), 5AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3), 5AAT(3*MX,3,3,3),GT(3*MX,3,3) LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM logical lxfile CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname DIMENSION NDNAME(IOX,MX),x(*),y(*),z(*) common/cname/asname(IOX,MX,3,3), 1psname(IOX,MX,3,3),aaname(IOX,3*MX,3,3,3),gname(IOX,3*MX,3,3), 2alphaname(IOX,3*MX,3,3),rname(IOX,3,MX),NATNAME(IOX),NAME, 3ffname(IOX,3*MX,3*MX) common/cctnopts/NS,NO,NATS, 1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,CUTOFF,n1,n2,okind, 3nmin,nmax,npoint,hw,temp,wsmin,wsmax,incr,nexp,navmx,nav character*5 s5 c open(3,file='CCT.OUT') nmin=1 nmax=NAT open(20,file='SPEC.OPT') 100 read(20,2000,err=111,end=111)s5 2000 format(a5) if(s5(1:4).eq.'NMIN')read(20,*)nmin if(s5(1:4).eq.'NMAX')read(20,*)nmax goto 100 111 close(20) NATi=nmax-nmin+1 c c NAT ... number of atoms in the tinekr molecule c NATi ... number of atoms in the tinker mol for transfer c NATS ... number of atoms in the small moelcule if(NATi.gt.MX)call report('too many atoms in cctn_tinker1') do 101 i=nmin,nmax ND(i-nmin+1)=atomic(i) R(1,i-nmin+1)=x(i) R(2,i-nmin+1)=y(i) R(3,i-nmin+1)=z(i) NBOB(i-nmin+1)=n12(i) do 101 j=1,NBOB(i-nmin+1) 101 IBONDSB(j,i-nmin+1)=i12(j,i) c Open file with option and overlaps and read it in: c For selected atoms only: OPEN(2,FILE='CCT.INP',STATUS='OLD') CALL readopt(NS,IB,IS,IND,MX,MXB,INDBIG,IOX,NO,NATi, 1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT,LAPT, 2LALF,LRDO,LNAMES,NAME,LPOLYMER,LHALTONERROR,CUTOFF,n1,n2,okind 3) CLOSE(2) c c if tensors constructed from named fragments, load them all: if(LNAMES)then DO 11 I=1,NO NE=80 do 10 NE=67,1,-1 10 if(NAME(I)(NE:NE).ne.' ')goto 12 12 basicname(1:NE)=NAME(I)(1:NE) xname=basicname(1:NE)//'.x' write(6,*)xname(1:NE+2) inquire(file=xname,exist=lxfile) if(.not.lxfile)then write(6,*)xname call report('Coordinates not found') endif open(2,file=xname) CALL READRXS(MX,NATS,RS,NDS,NBO) NATNAME(I)=NATS do 13 ia=1,NATS NDNAME(I,ia)=NDS(ia) do 13 ix=1,3 13 rname(I,ix,ia)=RS(ix,ia) 11 WRITE(6,*)NATS, ' atoms in the fragment' else if(lxfile)then OPEN(2,FILE='SMALL.X',STATUS='OLD') CALL READRXS(MX,NATS,RS,NDS,NBO) else OPEN(2,FILE='SMALL.XYZ',STATUS='OLD') CALL READRS(MX,NATS,RS,NDS,NBO) write(6,*)NATS,' atoms' endif CLOSE(2) endif c names c c check consistency in atomic types: do 22 i=1,NO if(LNAMES)then NATS=NATNAME(i) do 23 ia=1,NATS 23 NDS(ia)=ndname(i,ia) endif DO 22 IAB=1,NATi IF(INDBIG(I,IAB).NE.0)THEN IF(ND(IAB).NE.NDS(ABS(INDBIG(I,IAB))))THEN WRITE(6,3300)I,IAB,ND(IAB),NDS(ABS(INDBIG(I,IAB))) 3300 FORMAT('Overlap ',I2,', atom ',I3,': ',2I4) CALL REPORT(' Different atomic types !') ENDIF IF(IABS(INDBIG(I,IAB)).gt.NATS)then write(3,*)'overlap ',i write(3,*)' NATS ',NATS write(3,*)'IND ',IAB write(3,*)'INDB ',INDBIG(I,IAB) write(6,*)'overlap ',i write(6,*)' NATS ',NATS write(6,*)'IND ',IAB write(6,*)'INDB ',INDBIG(I,IAB) call report('Overflow NATS') ENDIF ENDIF 22 CONTINUE c c c Special section for named fragments: c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn if(lnames)then do 14 i=1,NO NATS=natname(i) NE=80 do 15 NE=80,1,-1 15 if(NAME(I)(NE:NE).ne.' ')goto 16 16 basicname=NAME(I)(1:NE) fname=basicname(1:NE)//'.fc' write(6,*)fname(1:NE+3) INQUIRE(FILE=fname,EXIST=OPT) IF(.NOT.OPT)THEN WRITE(6,*)'Force field of small not found' IF(LDIA.OR.LOFF)CALL REPORT('Cannot do force field') ELSE c c just check here, if it is possible to load it: OPEN(2,FILE=fname,STATUS='OLD') CALL READFF1(MX,NATS,FFS) CLOSE(2) WRITE(6,3002)NATS 3002 FORMAT(' Force field found,',/,I10,' atoms') do 21 ii=1,NATS*3 do 21 jj=1,NATS*3 c write(6,*)i,ii,jj,FFS(ii,jj) 21 ffname(i,ii,jj)=FFS(ii,jj) ENDIF c IF(LABS.OR.LVCD)THEN nname=basicname(1:NE)//'.ten' write(6,*)nname(1:NE+4) OPEN(15,FILE=nname,STATUS='OLD') do 19 ia=1,NATS do 19 ix=1,3 19 RS(ix,ia)=rname(i,ix,ia) CALL READTEN1(MX,PS,AS,LAPT,RS,NATS) CLOSE(15) do 18 ia=1,NATS do 18 ix=1,3 do 18 jx=1,3 asname(i,ia,ix,jx)=AS(ia,ix,jx) 18 psname(i,ia,ix,jx)=PS(ia,ix,jx) ENDIF C IF(LROA.or.LRAM)THEN tname=basicname(1:NE)//'.ttt' write(6,*)tname(1:NE+4) OPEN(15,FILE=tname,STATUS='OLD') CALL READTTT(MX,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO' CALL TTTT1(3*MX,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*MX,ALPHAS,AAS,GS,NATS,1,RS) ENDIF do 20 iax=1,3*NATS do 20 ix=1,3 do 20 jx=1,3 alphaname(i,iax,ix,jx)=ALPHAS(iax,ix,jx) gname(i,iax,ix,jx)=GS(iax,ix,jx) do 20 kx=1,3 20 aaname(i,iax,ix,jx,kx)=AAS(iax,ix,jx,kx) C Now small tensors are in atomic origins ENDIF C 14 continue c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn else INQUIRE(FILE='SMALL.FC',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,*)'Force field of small not found' WRITE(*,*)'Force field of small not found' IF(LDIA.OR.LOFF)CALL REPORT('Cannot do force field') ELSE OPEN(2,FILE='SMALL.FC',STATUS='OLD') CALL READFF1(MX,NATS,FFS) CLOSE(2) WRITE(3,3002)NATS WRITE(6,3002)NATS ENDIF IF(LABS.OR.LVCD)THEN OPEN(15,FILE='SMALL.TEN',STATUS='OLD') CALL READTEN1(MX,PS,AS,LAPT,RS,NATS) CLOSE(15) ENDIF C IF(LROA.OR.LRAM)THEN OPEN(15,FILE='SMALL.TTT',STATUS='OLD') CALL READTTT(MX,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO for small' CALL TTTT1(3*MX,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*MX,ALPHAS,AAS,GS,NATS,1,RS) ENDIF C Now small tensors are in atomic origins C IF(LRDO)THEN WRITE(3,*)' Alpha from SMALL_ALPHA.TTT for small' OPEN(15,FILE='SMALL_ALPHA.TTT',STATUS='OLD') CALL READTTT(MX,ALPHAS,AAT,GT,NATS) CLOSE(15) ENDIF endif ENDIF c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn c lnames c C write(6,*)'CCTN options read in' C c return end subroutine dospectrum(x,y,z,n12,i12,maxatm,maxval,istep) c IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind,n12(*),i12(maxval,maxatm) PARAMETER (MX=50,MXB=40,IOX=125) c MX maximum number of atoms c MXB maximum number of overlaped atoms for each constant c IOX max # of overlaps common/cctnrest/FF(3*MX,3*MX),R(3,MX),IND(MX),IB(MX,MXB),ND(MX), 1NDS(MX),RS(3,MX),FFS(3*MX,3*MX),IS(MX,MXB),AM(MX,3,3), 2INDBIG(IOX,MX),IBONDSB(4,MX), 3NBO(MX),NBOB(MX),PS(MX,3,3),AS(MX,3,3),PB(MX,3,3),AB(MX,3,3), 4ALPHA(3*MX,3,3),ALPHAS(3*MX,3,3),AA(3*MX,3,3,3), 5AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3), 5AAT(3*MX,3,3,3),GT(3*MX,3,3) LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM logical lxfile CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname DIMENSION NDNAME(IOX,MX),x(*),y(*),z(*) common/cctnopts/NS,NO,NATS, 1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,CUTOFF,n1,n2,okind, 3nmin,nmax,npoint,hw,temp,wsmin,wsmax,incr,nexp,navmx,nav c NATi=nmax-nmin+1 c c transcript current geometry: do 101 i=nmin,nmax R(1,i-nmin+1)=x(i) R(2,i-nmin+1)=y(i) 101 R(3,i-nmin+1)=z(i) DO 1 I=1,3*NATi DO 1 J=1,3*NATi 1 FF(I,J)=0.0d0 c IF(LABS.OR.LVCD)THEN DO 2 I=1,NATi DO 2 J=1,3 DO 2 K=1,3 PB(I,J,K)=0.0d0 2 AB(I,J,K)=0.0d0 ENDIF C IF(LROA.OR.LRAM)THEN DO 3 L=1,NATi DO 3 IX=1,3 IIND=3*(L-1)+IX DO 3 I=1,3 DO 3 J=1,3 ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 DO 3 K=1,3 3 AA(IIND,I,J,K)=0.0d0 ENDIF C CALL IMPROVE(FF,FFS,R,RS,NATi,NATS,NS,IND,IB,IS, 1AM,INDBIG,NO,IBONDSB,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA, 2IWG,LROA,LRAM,ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT, 3LNAMES,LHALTONERROR,CUTOFF,n1,n2,okind) c CONST=4.359828d0/0.529177d0/0.529177d0 c DO 4 I=1,3*NATi DO 4 J=1,3*NATi 4 FF(I,J)=FF(I,J)*CONST c c transform back: CALL TRTENVCD(MX,PB,AB,R,NATi) IF(LROA.or.LRAM)THEN CALL TTTT1(3*MX,ALPHA,AA,G,NATi,2,R) ENDIF call newa(NATi,FF,PB,AB,R,ALPHA,AA,G,istep) return END SUBROUTINE IMPROVE(FFB,FFS,RB,RS,NAT,NATS,NS,IND,IB,IS,AM,INDBIG, 1NO,IBONDS,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,IWG,LROA,LRAM, 2ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT,LNAMES, 3LHALTONERROR,CUTOFF,n1,n2,okind) implicit none INTEGER*4 MX,MXB,IOX PARAMETER (MX=50,MXB=40,IOX=125) INTEGER*4 okind,IND(MX),IS(MX,MXB),IB(MX,MXB),INDBIG(IOX,MX), 4IBONDS(4,MX),NBOB(MX),IOL(IOX),ISUM,I,IA,IAS,IEND,IERR,II, 1IIND,IIOII,IIX,IL,INDI,INDIB,INDIS,INDJ,INDS,INDJS,IO, 2IOIJ,IP,IPASS,IPO,ISAT,ISTART,IWG,IX,IXX,J,JA,JJ,IERC, 3JS,JSTART,JX,K,KK,N1,N2,NAT,NATIO,NATS,NE,NF,NO,NS,INDJB,JJX, 4ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,MX),nats_o(IOX) real*8 FFB(3*MX,3*MX),FFS(3*MX,3*MX),RB(3,MX),RS(3,NATS), 1TS(3),TB(3),AM(MX,3,3),B(3,3),TIJ(3),OTIJ(3),PB(MX,3,3), 2AB(MX,3,3),PS(MX,3,3),AS(MX,3,3),WF(IOX),ALPHA(3*MX,3,3), 3ALPHAS(3*MX,3,3),AA(3*MX,3,3,3),AAS(3*MX,3,3,3),G(3*MX,3,3), 4GS(3*MX,3,3),DISTM,SUM,FNEW,PNEW,ANEW,ALPHANEW,GNEW,AANEW, 5DIST,FOLD,POLD,AOLD,AD1,AD2,ALPHAOLD,GOLD,AD0,AAOLD,CTF,CUTOFF, 6RM(IOX,3,3) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LNAMES,LHALTONERROR real*8 A,XS,XB,TOL,RMS,FFNAME integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS real*8 aname,pname,aaname,gname,alphaname,rname integer*4 natname character*80 NAME(IOX) common/cname/aname(IOX,MX,3,3), 1pname(IOX,MX,3,3),aaname(IOX,3*MX,3,3,3),gname(IOX,3*MX,3,3), 2alphaname(IOX,3*MX,3,3),rname(IOX,3,MX),natname(IOX),NAME, 1ffname(IOX,3*MX,3*MX) c IF(NO.EQ.0)GOTO 19 ctf=CUTOFF**2 C C Look at all atom pairs i-j: istart=1 iend=NAT if(n1.gt.0)then istart=n1 iend=n2 endif DO 20 IA= istart,iend DO 20 JA=IA,NAT IF(.NOT.LOFF.AND.IA.NE.JA)GOTO 20 c c look at distance, skip for a large one if(CUTOFF.gt.0.0d0)THEN dist=(RB(1,IA)-RB(1,JA))**2+(RB(2,IA)-RB(2,JA))**2+ 1 (RB(3,IA)-RB(3,JA))**2 if(dist.gt.ctf)then goto 20 endif endif C C Calculate center of i/j radiusvectors: DO 22 IX=1,3 22 TIJ(IX)=0.5d0*(RB(IX,IA)+RB(IX,JA)) C C Look at all overlaps io, find that which fits best the i-j pair: DISTM=1.0d5 NF=0 DO 21 IO=1,NO C C Check, if io contains ij: IF(INDBIG(IO,IA).GT.0.AND.INDBIG(IO,JA).GT.0)THEN NF=NF+1 IOL(NF)=IO C C Calculate center of the overlap and distance from i-j center c i.........TIJ.............j c ----------OTIJ--------------------- DO 23 IX=1,3 23 OTIJ(IX)=0.0D0 NATIO=0 DO 24 II=1,NAT IIOII=INDBIG(IO,II) C If Lstrict, take only transformed atoms IF(LSTRICT.AND.IIOII.LT.0)IIOII=0 IIOII=ABS(IIOII) IF(IIOII.NE.0)THEN NATIO=NATIO+1 DO 25 IX=1,3 25 OTIJ(IX)=OTIJ(IX)+RB(IX,II) ENDIF 24 CONTINUE c c If required, calculate RMS overlap and use this criterium instead c of the distances,remember rotation matrix (RM) if(okind.eq.1)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IO, 1 IA,JA,NBOB,1,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) do 55 IX=1,3 do 55 JX=1,3 55 RM(NF,IX,JX)=A(IX,JX) ipass_o(NF)=IPASS WF(NF)=DSQRT(RMS) ierc_o(NF)=IERC itu_o(NF)=ITU do 58 I=1,ITU 58 ind_o(NF,I)=IND(I) do 59 I=1,3 do 59 J=1,3 59 RM(NF,I,J)=A(I,J) nats_o(NF)=NATS IF(RMS.LT.DISTM)THEN DISTM=RMS IOIJ=IO ENDIF c c Use distances: else DO 26 IX=1,3 26 OTIJ(IX)=OTIJ(IX)/DBLE(NATIO) DIST=(OTIJ(1)-TIJ(1))**2+(OTIJ(2)-TIJ(2))**2+(OTIJ(3)-TIJ(3))**2 IF(DIST.LT.DISTM)THEN DISTM=DIST IOIJ=IO ENDIF WF(NF)=DSQRT(DIST) endif ENDIF 21 CONTINUE c IF(NF.EQ.0)THEN GOTO 20 ELSE C C Calculate the weighting factors: SUM=0.0d0 DO 43 I=1,NF 43 IF(WF(I).GT.1.0d-9)SUM=SUM+1.0d0/WF(I) DO 42 I=1,NF IF(WF(I).LE.1.0d-9)THEN WF(1)=1.0d0 NF=1 GOTO 44 ENDIF WF(I)=1.0d0/WF(I)/SUM 42 CONTINUE 44 CONTINUE C C Loop over possible overlaps: ISUM=0 DO 45 IP=1,NF IPO=IOL(IP) IF(IWG.EQ.0)THEN IF(IPO.NE.IOIJ)GOTO 45 WF(IP)=1.0d0 ENDIF ISUM=ISUM+1 C C Find the atomic groups to be placed on together for IA-JA: c and the rotation matrix if(okind.eq.0)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1 IA,JA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) else IPASS=ipass_o(IP) IERC=ierc_o(IP) ITU=itu_o(IP) do 56 I=1,ITU 56 IND(I)=ind_o(IP,I) do 57 I=1,3 do 57 J=1,3 57 A(I,J)=RM(IP,I,J) NATS=nats_o(IP) endif if(IERC.ne.0)goto 20 C IF(LINV)THEN DO 50 I=1,3 DO 50 J=1,3 50 A(I,J)=-A(I,J) ENDIF C C Transform the force field: IF((LDIA.AND.IA.EQ.JA).OR.(LOFF.AND.IA.NE.JA))THEN if(LNAMES)then c retrieve FF of fragment IPO NATS=NATNAME(IPO) do ix=1,3*NATS do jx=1,3*NATS FFS(ix,jx)=ffname(IPO,ix,jx) enddo enddo ENDIF DO 36 IX=1,3 INDIB=3*IA-3+IX JSTART=1 IF(IA.EQ.JA)JSTART=IX DO 36 JX=JSTART,3 INDJB=3*JA-3+JX FOLD=FFB(INDIB,INDJB) C Zero out at the start of the summation IF(ISUM.EQ.1)FFB(INDIB,INDJB)=0.0d0 FNEW=0.0d0 DO 37 IIX=1,3 INDIS=3*INDBIG(IPO,IA)-3+IIX DO 37 JJX=1,3 INDJS=3*INDBIG(IPO,JA)-3+JJX 37 FNEW=FNEW+A(IX,IIX)*FFS(INDIS,INDJS)*A(JX,JJX) c IF(LWR)WRITE(3,3001)IA,JA, IX,JX,FOLD,FNEW,IIT,TOL,WF(IP) FFB(INDIB,INDJB)=FNEW*WF(IP)+FFB(INDIB,INDJB) FFB(INDJB,INDIB)=FFB(INDIB,INDJB) 36 CONTINUE ENDIF C C Transform the dipole derivatives IF(IA.EQ.JA.AND.(LABS.OR.LVCD))THEN IAS=INDBIG(IPO,IA) DO 40 IX=1,3 DO 40 JX=1,3 POLD=PB(IA,IX,JX) AOLD=AB(IA,IX,JX) C Zero out at the start of the summation IF(ISUM.EQ.1)THEN PB(IA,IX,JX)=0.0d0 AB(IA,IX,JX)=0.0d0 ENDIF PNEW=0.0d0 ANEW=0.0d0 DO 41 II=1,3 AD1=A(IX,II) DO 41 JJ=1,3 AD2=AD1*A(JX,JJ) if(LNAMES)then PNEW=PNEW+pname(IPO,IAS,II,JJ)*AD2 ANEW=ANEW+aname(IPO,IAS,II,JJ)*AD2 else PNEW=PNEW+PS(IAS,II,JJ)*AD2 ANEW=ANEW+AS(IAS,II,JJ)*AD2 endif 41 continue PB(IA,IX,JX)=PNEW*WF(IP)+PB(IA,IX,JX) AB(IA,IX,JX)=ANEW*WF(IP)+AB(IA,IX,JX) IF (.NOT.LVCD)AB(IA,IX,JX)=AOLD IF (.NOT.LVCD)ANEW=AOLD IF (.NOT.LABS)PB(IA,IX,JX)=POLD IF (.NOT.LABS)PNEW=POLD 40 continue ENDIF C C Transform the ROA and Raman tensors IF(IA.EQ.JA.AND.(LROA.or.LRAM))THEN IAS=INDBIG(IPO,IA) DO 46 IX=1,3 IIND=3*(IA-1)+IX DO 46 I=1,3 DO 46 J=1,3 ALPHAOLD=ALPHA(IIND,I,J) GOLD=G(IIND,I,J) C Zero out at the start of the summation IF(ISUM.EQ.1)THEN ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 ENDIF ALPHANEW=0.0d0 GNEW=0.0d0 DO 47 IXX=1,3 AD0=A(IX,IXX) INDS=3*(IAS-1)+IXX DO 47 II=1,3 AD1=AD0*A(I,II) DO 47 JJ=1,3 AD2=AD1*A(J,JJ) if(LNAMES)then ALPHANEW=ALPHANEW+alphaname(IPO,INDS,II,JJ)*AD2 GNEW=GNEW+gname(IPO,INDS,II,JJ)*AD2 else ALPHANEW=ALPHANEW+ALPHAS(INDS,II,JJ)*AD2 GNEW=GNEW+GS(INDS,II,JJ)*AD2 endif 47 continue ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+ALPHANEW*WF(IP) G(IIND,I,J)=G(IIND,I,J)+GNEW*WF(IP) IF (.NOT.LRAM)ALPHA(IIND,I,J)=ALPHAOLD IF (.NOT.LRAM)ALPHANEW =ALPHAOLD IF (.NOT.LROA) G(IIND,I,J)=GOLD IF (.NOT.LROA) GNEW=GOLD 46 continue c6 IF(LWR)WRITE(3,3006)IIND,I,J,ALPHAOLD,ALPHANEW, c 1 GOLD,GNEW,WF(IP) DO 48 IX=1,3 IIND=3*(IA-1)+IX DO 48 I=1,3 DO 48 J=1,3 DO 48 K=1,3 AAOLD=AA(IIND,I,J,K) C Zero out at the start of the summation IF(ISUM.EQ.1)AA(IIND,I,J,K)=0.0d0 AANEW=0.0d0 DO 49 IXX=1,3 AD0=A(IX,IXX) INDS=3*(IAS-1)+IXX DO 49 II=1,3 AD1=AD0*A(I,II) DO 49 JJ=1,3 AD2=AD1*A(J,JJ) DO 49 KK=1,3 if(LNAMES)then AANEW=AANEW+aaname(IPO,INDS,II,JJ,KK)*AD2*A(K,KK) else AANEW=AANEW+AAS(INDS,II,JJ,KK)*AD2*A(K,KK) endif 49 continue AA(IIND,I,J,K)=AANEW*WF(IP)+AA(IIND,I,J,K) IF (.NOT.LROA)AA(IIND,I,J,K)=AAOLD IF (.NOT.LROA)AANEW=AAOLD 48 continue ENDIF 45 CONTINUE ENDIF C 20 CONTINUE c IA,JA C C If only interaction among transfered atoms selcted, the other C force constants will be set to zero in final force field: IF(LSELECT)THEN DO 51 IA=1,NAT DO 52 IO=1,NO 52 IF(INDBIG(IO,IA).GT.0)GOTO 51 c WRITE(3,*)' Atom ',IA,' taken out of the FF' DO 53 IX=1,3 INDIB=3*IA-3+IX DO 53 JA=1,NAT DO 53 JX=1,3 INDJB=3*JA-3+JX FFB(INDIB,INDJB)=0.0d0 53 FFB(INDJB,INDIB)=0.0d0 51 CONTINUE ENDIF RETURN C 19 WRITE(3,3000) 3000 FORMAT(1X,' Old matrix element New element',/, 1 1X,' atoms xi xj F [a.u.] F',5H' #it,/, 2 1X,'-------------------------------------------------') DO 1 IA=1,NS DO 3 I=1,3 TB(I) =0.0d0 3 TS(I)=0.0d0 DO 4 ISAT=1,ABS(IND(IA))+1 DO 4 J=1,3 XS(ISAT,J)=RS(J,IS(IA,ISAT)) XB(ISAT,J)=RB(J,IB(IA,ISAT)) TB(J)=TB(J)+XB(ISAT,J) 4 TS(J)=TS(J)+XS(ISAT,J) DO 5 ISAT=1,ABS(IND(IA))+1 DO 5 J=1,3 XS(ISAT,J)=XS(ISAT,J)-TS(J)/DBLE(ABS(IND(IA))+1.0d0) 5 XB(ISAT,J)=XB(ISAT,J)-TB(J)/DBLE(ABS(IND(IA))+1.0d0) IF(IND(IA).LT.0)THEN WRITE(3,*)' Space inversion applied for vicinity of atom ',IA DO 9 ISAT=1,ABS(IND(IA))+1 DO 9 J=1,3 9 XS(ISAT,J)=-XS(ISAT,J) ENDIF ITU=ABS(IND(IA))+1 CALL DOMATRIX(IERR) IF(IERR.NE.0)CALL REPORT('Erroneous stop') DO 10 I=1,3 DO 10 J=1,3 10 AM(IA,I,J)=A(I,J) DO 1 JA=1,IA DO 11 I=1,3 DO 11 J=1,3 11 B(I,J)=AM(JA,I,J) DO 6 I=1,3 INDI=3*IB(IA,1)-3+I JSTART=1 IF(IA.EQ.JA)JSTART=I DO 6 J=JSTART,3 INDJ=3*IB(JA,1)-3+J FOLD=FFB(INDI,INDJ) FFB(INDI,INDJ)=0.0d0 DO 7 II=1,3 IL=3*IS(IA,1)-3+II DO 7 JJ=1,3 JS=3*IS(IA,1)-3+JJ 7 FFB(INDI,INDJ)=FFB(INDI,INDJ)+A(I,II)*FFS(IL,JS)*B(J,JJ) IF(I.EQ.1.AND.J.EQ.I.AND.JA.EQ.IA)THEN WRITE(3,3001)IA,JA,I,J,FOLD,FFB(INDI,INDJ),IIT ELSE WRITE(3,3001)IA,JA,I,J,FOLD,FFB(INDI,INDJ) ENDIF FFB(INDJ,INDI)=FFB(INDI,INDJ) 6 CONTINUE 3001 FORMAT(4I5,2F11.5,I5,F10.7,F7.4) 1 CONTINUE RETURN END c SUBROUTINE READRS(N0,N,R,ND,NB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),NB(N0) CHARACTER*78 TITDUM READ(2,2000) TITDUM 2000 FORMAT(A78) N=0 1 N=N+1 READ(2,*)ND(N),R(1,N),R(2,N),R(3,N) NB(N)=0 IF(ND(N).EQ.9999)GOTO 2 GOTO 1 2 N=N-1 IF(N.GT.N0)CALL REPORT(' N > MX !') RETURN END c SUBROUTINE READRXS(N0,N,R,ND,NB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),NB(N0) CHARACTER*78 TITDUM READ(2,2000) TITDUM 2000 FORMAT(A78) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > MX !') do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i) 1 NB(i)=0 RETURN END SUBROUTINE READFF1(N0,NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*N0,3*N0) NAT3=3*NA N=NAT3 C Read in the lower triangle of FF, C written in parts as n1,n1 C . . C ln,n1 . ln,ln C . . . C . . . n3,n3 C . . . . C n,n1 . . n,n3 N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(2,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C DO 31 I=1,N DO 31 J=I,N 31 FCAR(I,J)=FCAR(J,I) RETURN END SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=40) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER,RMS IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=1.0d-7 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)then RETURN endif IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) implicit none INTEGER*4 MXB,NAT,ITER,I REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,A,XS,XB,RTOL,RMS,ANG(3),FU PARAMETER (MXB=40) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER,RMS S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 RMS=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 RMS=RMS+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=RMS RETURN END C c SUBROUTINE SREPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(3,3000) WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(3,*)RS WRITE(6,*)RS WRITE(3,3001) WRITE(6,3001) 3001 FORMAT(80(1H*),/,/, 1'PROGRAM CONTINUES BECAUSE LHALTONERROR DISABLED') END C SUBROUTINE READTEN1(N0,P,A,LAPT,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=50) DIMENSION P(N0,3,3),A(N0,3,3),AI(MX,3,3),AJ(MX,3,3),PV(MX,3,3), 1R(3,MX),X(3,NAT) LOGICAL LAPT BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR READ (15,*) NOAT NAT=NOAT WRITE(3,*)NAT,' atoms' WRITE(*,*)NAT,' atoms' C C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) WRITE(6,*)' Dipole derivatives read - in ' WRITE(3,*)' Dipole derivatives read - in ' C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) AJ(L,J,I)=0.0d0 PV(L,J,I)=P(L,J,I) 2203 AI(L,J,I)=SUM WRITE(6,*)' APT part of AAT constructed' WRITE(3,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (AI(L,J,I),I=1,3) DO 230 L=1,NOAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NOAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) WRITE(6,*)' VCD parameters read in ...' WRITE(3,*)' VCD parameters read in ...' C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) WRITE(6,*)'DOG used' WRITE(3,*)'DOG used' C C Total axial tensor = electronic + nuclear: DO 4 L=1,NAT DO 4 I=1,3 DO 4 J=1,3 4 A(L,I,J)=AI(L,I,J)+AJ(L,I,J) C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(3,*)' AAT shifted to atomic origin' WRITE(6,*)' AAT shifted to atomic origin' RETURN END C SUBROUTINE TRTENVCD(N0,P,A,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=50) DIMENSION P(N0,3,3),A(N0,3,3),R(3,MX),X(3,NAT) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM RETURN END C SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA PARAMETER (MX=50) DIMENSION VCD(MX,3,3),VEL(MX,3,3),DIP(MX,3,3), 1C(3,NAT) DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END C SUBROUTINE TTTT1(N,ALPHA,A,G,NAT,ICO,R) C Transforms A and G to local origins (ICO=1) or c common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,MX) C C transfer into atomic coordinates: DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I,J)/0.529177D0 C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID) G(IIND,IA,IB)=G(IIND,IA,IB)+SUM IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0 2 CONTINUE c IF(ICO.EQ.1)WRITE(3,*)' G transformed into atomic origins' c IF(ICO.EQ.2)WRITE(3,*)' G transformed into common origin' c IF(ICO.EQ.3)WRITE(3,*)' G set to zero' C DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM 1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0 3 CONTINUE c IF(ICO.EQ.1)WRITE(3,*)' A transformed into atomic origins' c IF(ICO.EQ.2)WRITE(3,*)' A transformed into common origin' c IF(ICO.EQ.3)WRITE(3,*)' A set to zero' RETURN END c SUBROUTINE readopt(NS,IB,IS,IND,MXT,MXB,INDBIG,IOX,NO, 1NAT,LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT, 2LAPT,LALF,LRDO,LNAMES,NAME,LPOLYMER,LHALTONERROR,CUTOFF, 2n1,n2,okind) implicit none INTEGER*4 NS,IB,IS,IND,MXT,MXB,INDBIG,IOX,NO, 1NAT,n1,n2,okind,I,IA,IAB,IC,IEND,II,IO,ISMALL,IWG,J,JA,NDD REAL*8 CUTOFF DIMENSION IB(MXT,MXB),IS(MXT,MXB),IND(MXT),INDBIG(IOX,MXT), 1NDD(20) CHARACTER*7 STR,STRC CHARACTER*80 NAME(IOX) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT,LAPT, 1LALF,LRDO,LPOLYMER,LNAMES,LHALTONERROR,LRAM c n1=0 n2=0 okind=1 LRDO=.FALSE. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LOFF=.TRUE. LDIA=.TRUE. LHALTONERROR=.TRUE. CUTOFF=-1.0d5 LROA=.FALSE. LRAM=.FALSE. LSTRICT=.TRUE. LINV=.FALSE. LSELECT=.FALSE. LNAMES=.FALSE. IWG=1 9 READ(2,2000)STR 2000 FORMAT(A7) WRITE(3,2000)STR IF(STR(1:4).EQ.'N1N2') READ(2,*)n1,n2 IF(STR(1:3).EQ.'OKI') READ(2,*)okind IF(STR(1:3).EQ.'IWG') READ(2,*)IWG IF(STR(1:6).EQ.'LNAMES') READ(2,*)LNAMES IF(STR.EQ.'LSELECT') READ(2,*)LSELECT IF(STR(1:4).EQ.'LINV') READ(2,*)LINV IF(STR(1:4).EQ.'LAPT') READ(2,*)LAPT IF(STR.EQ.'LSTRICT') READ(2,*)LSTRICT IF(STR(1:3).EQ.'LWR') READ(2,*)LWR IF(STR(1:4).EQ.'LALF') READ(2,*)LALF IF(STR(1:4).EQ.'LRDO') READ(2,*)LRDO IF(STR(1:4).EQ.'LRAM') READ(2,*)LRAM IF(STR(1:4).EQ.'LROA') READ(2,*)LROA IF(STR(1:4).EQ.'LOFF') READ(2,*)LOFF IF(STR(1:4).EQ.'LDIA') READ(2,*)LDIA IF(STR(1:4).EQ.'LABS') READ(2,*)LABS IF(STR(1:4).EQ.'LVCD') READ(2,*)LVCD IF(STR(1:7).EQ.'LHALTON')READ(2,*)LHALTONERROR IF(STR(1:6).EQ.'CUTOFF') READ(2,*)CUTOFF BACKSPACE 2 READ(2,2000)STRC if(STR.ne.STRC)goto 9 NO=0 c c older laborius definition of overlap: IF(STR.NE.'POLYMER')THEN LPOLYMER=.false. REWIND 2 READ(2,*)NS DO 1 I=1,NS READ(2,*)IB(I,1),IS(I,1),IND(I) DO 1 J=1,ABS(IND(I)) 1 READ(2,*)IB(I,J+1),IS(I,J+1) C C NS atoms are assigned C IB(I) atom of the big molecule goes with the IS(I) atom of the smaller. C Each of them has IND(I) satelite atoms to define the reference frame. C Reference atoms are given for big as IB and for small as IS strings. C c automatic definition of overlap: ELSE LPOLYMER=.true. READ(2,*)NO WRITE(3,*)NO,' overlaps found, atoms assigned automatically ...' IF(NO.GT.IOX)CALL REPORT(' NO > IOX ! Change the dimensions.') DO 2 I=1,NO READ(2,*)IC IF(IC.NE.I)THEN WRITE(6,*)IC,I,' nat = ',NAT CALL REPORT('Error in reading ') ENDIF IF(LNAMES)READ(2,2900)NAME(I) IF(LNAMES)write(6,2900)NAME(I) 2900 FORMAT(a80) c write(6,*)NAT READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT) 2001 FORMAT(20I3) DO 5 IAB=1,NAT 5 INDBIG(I,IAB)=0 READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT) c write(6,2001)(INDBIG(I,IAB),IAB=1,NAT) 2 CONTINUE c c Control Output: DO 4 IO=1,NO WRITE(3,2002)IO 2002 FORMAT(' Overlap number ',I5) IA=0 7 IA=IA+1 IEND=IA+19 IF(IEND.GT.NAT)IEND=NAT IF(IA.GT.NAT)GOTO 6 DO 8 II=1,(IEND-IA+1) 8 NDD(II)=II+IA-1 WRITE(3,2004)(NDD(II),II=1,IEND-IA+1) 2004 FORMAT(20I3) IA=IEND GOTO 7 6 WRITE(3,2001)(INDBIG(IO,IAB),IAB=1,NAT) 4 CONTINUE c c check, that atoms are not defined twice in the overlaps: do io=1,no do ia=1,nat ismall=indbig(io,ia) if(ismall.ne.0)then do ja=ia+1,nat if(ismall.eq.indbig(io,ja))then write(3,*)'overlap ',io write(3,*)'big ',ia,ja write(3,*)'small ',ismall call report('umbiqous overlap') endif enddo endif enddo enddo ENDIF RETURN END C subroutine CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1IA,JA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS) c find covalently bond atoms to IA..JA, oevrlap and make the c rotation matrix c IC=0 ... take atoms two bonds away only if directly c bonded do not suffice c IC=1 ... include all by default implicit none integer*4 IOX,MXB,MX PARAMETER (MXB=40) integer*4 IA,IPASS,II,NAT,INDBIG(IOX,MX),IBA,IBONDS(4,MX), 1NBOB(*),JB,IIA,NATS,natname(*),IX,IND(*),iu,IERC,I,IC,IERR, 2IPO,J,JA real*8 RMS,RS(3,MX),RB(3,MX),TB(3),TS(3),rname(IOX,3,MX) logical LNAMES,LHALTONERROR,LSTRICT,LINV real*8 A,XS,XB,TOL integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS IERC=1 IPASS=0 666 ITU=0 IPASS=IPASS+1 C DO 27 II=1,NAT if(INDBIG(IPO,II).NE.0)then C Take IA and JA right away in the list: IF(II.EQ.IA.OR.II.EQ.JA)GOTO 277 C Take also atoms connected to IA or JA: DO 29 IBA=1,NBOB(II) 29 IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA)GOTO 277 C For second pass, take also atoms two bonds away: IF(IPASS.EQ.2.or.IC.eq.1)THEN DO 31 IBA=1,NBOB(II) IIA=IBONDS(IBA,II) DO 31 JB=1,NBOB(IIA) 31 IF(IBONDS(JB,IIA).EQ.IA.OR.IBONDS(JB,IIA).EQ.JA)GOTO 277 ENDIF c c This atom is of no use, try another one: GOTO 27 277 IF(INDBIG(IPO,II).LT.0.AND.LSTRICT)GOTO 27 C Include atom II into the transformation set: ITU=ITU+1 IF(ITU.GT.MXB)CALL REPORT('ITU > MXB !') if(LNAMES)NATS=natname(IPO) DO 28 IX=1,3 if(LNAMES)then XS(ITU,IX)=rname(IPO,IX,ABS(INDBIG(IPO,II))) else XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II))) endif 28 XB(ITU,IX)=RB(IX,II) IND(ITU)=II endif 27 CONTINUE c c IF(ITU.LT.3)WRITE(3,*)' ITU < 3 !' IF(ITU.LT.3.AND.IPASS.EQ.1.and.IC.eq.0)GOTO 666 if(ITU.LT.3)then write(6,*)'Only following atoms found in the big molecule:' write(6,*)(IND(iu) ,iu=1,ITU) if(LHALTONERROR)then CALL REPORT(' ITU < 3 !') else CALL SREPORT(' ITU < 3 !') return endif endif c c C Calculate center of the coordinates that will be rotated: DO 33 I=1,3 TB(I) =0.0d0 33 TS(I)=0.0d0 DO 34 II=1,ITU DO 34 J=1,3 TB(J)=TB(J)+XB(II,J) 34 TS(J)=TS(J)+XS(II,J) DO 38 J=1,3 TB(J)=TB(J)/DBLE(ITU) 38 TS(J)=TS(J)/DBLE(ITU) C Subtract mass-center: DO 35 II=1,ITU DO 35 J=1,3 XS(II,J)=XS(II,J)-TS(J) 35 XB(II,J)=XB(II,J)-TB(J) C IF(LINV)THEN DO 39 II=1,ITU DO 39 J=1,3 39 XS(II,J)=-XS(II,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR) IF(IERR.EQ.1.AND.IPASS.EQ.1)THEN write(6,*)'domat error,try other pass' GOTO 666 ENDIF IF(IERR.EQ.1)then if(LHALTONERROR)then CALL REPORT('Program cannot continue') else CALL SREPORT('Program cannot continue') return endif endif IERC=0 return end subroutine loadnewaoptions(x,y,z,atomic,spalpha) PARAMETER (MX=50,MX3=3*MX,np0=10000) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) logical lex,lopt integer*4 atomic(*) CHARACTER*5 s5 CHARACTER*80 s80 common/opts/lopt(300) common/allopts/kt(MX),r(3,MX),XM(MX),iso(MX),riso(MX), 1wexp(MX),CHAR(MX),SMAT(MX3,MX3),FRQ(MX3), 1BOHR,wmin,wmax,EXCA, 3stable(3*MX,5) real*8 x(*),y(*),z(*) logical LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV, 1LSELECT,LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR integer*4 okind common/cctnopts/NS,NO,NATS, 1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,CUTOFF,n1,n2,okind, 3nmin,nmax,npoint,hw,temp,wsmin,wsmax,incr,nexp,navmx,nav common/spec/spectra(np0,4),espx(np0,4),espy(np0,4),nsp(4), 1iskind(4),xmin(4),xmax(4),avspectra(np0,4) c c stable c index 1 - freq 2-D 3-R 4-Raman 5-ROA/180 C c 1 . calculate frequencies c 2 . generate FTRY.INP c 3 . substitute acidic hydrogens by deuteria c 4 . read in FTRY.INP c 5 . project zero vibrations c 6 . write F.INP c 7 . do FPC c 8 . call new5/abs/vcd c 9 . write DOG.TAB c 10 . write DIP.MOM c 11 . read F.INP c 12 . write dipoles in F.INP c 13 . call new6/Raman/ROA c 14 . polar ROA c 15 . empty c 16 . write ROA.TAB c 17 . large Raman print c 18 . generate abs spectrum c 19 . write abs spectrum do d.prn c 20 . generate vcd spectrum c 21 . write vcd spectrum do d.prn c 22 . generate Raman spectrum c 23 . write Raman spectrum do d.prn c 24 . generate ROA spectrum c 25 . write ROA spectrum do d.prn c 26 . Gaussian c 27 . Gradient from average instead of immediate spectra c c WMIN, WMAX min, max frequency for new5,new6 c WSMIN, WSMAX min, max frequency for theoretical spectra c XMIN, XMAX min, max frequency for exp/comparison c c EXCA excitation wavelength for Raman in A c c hw ... peak width c navmx ... max. No. of averag. spectra c npoint ... No. of spectral points c temp ... temperature c incr write spectra each incr step c NEXP ... # of experimental spectra c np0 ... max. No. of spectral points c spalpha ... energy term from the spectra, E ~ "spalpha*(s-sexp)^2" c BOHR=0.52917715d0 do 3 i=1,300 3 lopt(i)=.false. lopt(1)=.true. lopt(5)=.true. EXCA=4880.0d0 wmin=10.0d0 wmax=1.0d9 wsmin=100.0d0 wsmax=4000.0d0 lopt(26)=.false. npoint=2000 navmx=10 hw=10.0d0 temp=298.15d0 incr=100 nexp=0 nav=0 spalpha=0.1d0 do 9 i=1,4 xmin(i)=0.0d0 xmax(i)=4000.0d0 do 9 j=1,np0 9 avspectra(j,i)=0.0d0 c inquire(file='SPEC.OPT',exist=lex) if(lex)then open(20,file='SPEC.OPT') 1 read(20,200,end=99,err=99)s5 if(s5.eq.'LFREQ')read(20,*)lopt(1) if(s5(1:4).eq.'FTRY')read(20,*)lopt(2) if(s5(1:4).eq.'DEUT')read(20,*)lopt(3) if(s5(1:4).eq.'RFTR')read(20,*)lopt(4) if(s5(1:4).eq.'PROJ')read(20,*)lopt(5) if(s5(1:5).eq.'F.INP')read(20,*)lopt(6) if(s5(1:3).eq.'FPC')read(20,*)lopt(7) if(s5(1:3).eq.'ABS')read(20,*)lopt(8) if(s5(1:3).eq.'DOG')read(20,*)lopt(9) if(s5(1:3).eq.'MOM')read(20,*)lopt(10) if(s5(1:3).eq.'R.F')read(20,*)lopt(11) if(s5(1:3).eq.'W.F')read(20,*)lopt(12) if(s5(1:3).eq.'RAM')read(20,*)lopt(13) if(s5(1:3).eq.'POL')read(20,*)lopt(14) if(s5(1:3).eq.'ROA')read(20,*)lopt(16) if(s5(1:5).eq.'GAUSS')read(20,*)lopt(26) if(s5(1:2).eq.'##')then read(20,*)io read(20,*)lopt(io) endif if(s5(1:4).eq.'WMIN')read(20,*)wmin if(s5(1:4).eq.'WMAX')read(20,*)wmax if(s5(1:4).eq.'WSMI')read(20,*)wsmin if(s5(1:4).eq.'WSMA')read(20,*)wsmax if(s5(1:4).eq.'EXCA')read(20,*)EXCA if(s5(1:2).eq.'HW')read(20,*)hw if(s5(1:4).eq.'TEMP')read(20,*)temp if(s5(1:4).eq.'INCR')read(20,*)incr if(s5(1:5).eq.'ALPHA')read(20,*)spalpha if(s5(1:5).eq.'AVERA')read(20,*)lopt(27) if(s5(1:4).eq.'ISOT')then read(20,*)ISOT read(20,*)(iso(i),i=1,ISOT) read(20,*)(riso(i),i=1,ISOT) endif if(s5(1:5).eq.'NAVMX')read(20,*)navmx if(s5(1:5).eq.'NPOIN')read(20,*)npoint if(s5(1:4).eq.'NEXP')then read(20,*)nexp if(nexp.gt.4)call report('too many spectra') do 5 i=1,nexp read(20,240)s80 write(6,240)s80 240 format(a80) read(20,*)xmin(i),xmax(i) write(6,*)xmin(i),xmax(i) read(20,*)iskind(i) ii=0 open(21,file=s80) 4 read(21,*,end=211,err=211)xx,yy if(xx.ge.xmin(i).and.xx.le.xmax(i))then ii=ii+1 espx(ii,i)=xx espy(ii,i)=yy endif goto 4 211 close(21) write(6,*)ii,' points' 5 nsp(i)=ii endif goto 1 200 format(a5) 99 close(20) endif write(6,*)'newa options loaded' do 2 I=nmin,nmax kt(I-nmin+1)=atomic(I) r(1,I-nmin+1)=x(I) r(2,I-nmin+1)=y(I) 2 r(3,I-nmin+1)=z(I) call INPTRY(XM,MX,nmax-nmin+1,kt,r,iso,riso,ISOT,wexp,nm,CHAR) return end c ==================================================================== subroutine newa(NN,FCAR,PP,II,RI,ALPHA,A,G,istep) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50,MX3=3*MX,np0=10000) REAL*8 II(MX,3,3),PP(MX,3,3),FCAR(MX3,MX3),kk, 2ALPHA(MX3,3,3),A(MX3,3,3,3),G(MX3,3,3),RI(3,MX) logical lex,lopt CHARACTER*5 s5 common/opts/lopt(300) common/allopts/kt(MX),r(3,MX),XM(MX),iso(MX),riso(MX), 1wexp(MX),CHAR(MX),SMAT(MX3,MX3),FRQ(MX3), 1BOHR,wmin,wmax,EXCA, 3stable(3*MX,5) common/spec/spectra(np0,4),espx(np0,4),espy(np0,4),nsp(4), 1iskind(4),xmin(4),xmax(4),avspectra(np0,4) integer*4 NS,NO,NATS,IWG,n1,n2,okind,nmin,nmax,npoint,navmx, 1nav real*8 CUTOFF,hw,temp,wsmin,wsmax logical LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV, 1LSELECT,LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR common/cctnopts/NS,NO,NATS, 1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,CUTOFF,n1,n2,okind, 3nmin,nmax,npoint,hw,temp,wsmin,wsmax,incr,nexp,navmx,nav character*6 s6 real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0) COMMON/MATVARs/s1x,s1y,s2x,s2y,n11,n22,RTOL,ITER,RMS,aa,kk,b, 1nss,jsc include 'sizes.i' include 'atoms.i' c c spectral gradient: logical lexsp real*8 gradsp,xold,yold,zold,ensp,spalpha,sumspold,rshold common /spegrad/gradsp(3*maxatm),xold(maxatm),yold(maxatm), 1zold(maxatm),rshold(maxatm),ensp,lexsp,spalpha,sumspold real*8 rsh(3*maxatm) DO 35 K=1,NN DO 35 I=1,3 35 R(I,K)=RI(I,K)/BOHR do 1 i=1,3*NN do 1 j=1,5 1 stable(i,j)=0.0d0 nm=3*NN call new4(NN,FCAR,CHAR,XM,wexp,nm,SMAT,FRQ,R) c c absorption/vcd spectra if(lopt(8))call new5(PP,II,NN,r,SMAT,FRQ,nm,stable) if(lopt(13)) 1call new6(SMAT,wmin,wmax,ALPHA,G,A,NN,EXCA,nm,FRQ,r,stable) c do i=1,nm c write(6,609)i,(stable(i,j),j=1,5) c09 format(i4,f8.2,4g12.4) c enddo c c generate spectra if(npoint.gt.np0)call report('too many spectral points') do 3 j=1,4 do 3 i=1,npoint 3 spectra(i,j)=0.0d0 dw=(wsmax-wsmin)/dble(npoint-1) if(lopt(26))then dd=hw/2.0d0/dsqrt(log(2.0d0)) else dd=hw/2.0d0 endif do 2 i=1,npoint xx=wsmin+dw*dble(i-1) do 2 i1=1,nm w=stable(i1,1) if(lopt(18))then c abs c write(6,*)stable(i1,2),hw,lopt(26),xx,w, c 1 sf(hw,lopt(26),xx,w) const=108.7d0*stable(i1,2) spectra(i,1)=spectra(i,1)+sf(dd,lopt(26),xx,w)*const endif if(lopt(20))then c vcd const=435.0d0*stable(i1,3) spectra(i,2)=spectra(i,2)+sf(dd,lopt(26),xx,w)*const endif if(lopt(22))then c Raman const=stable(i1,4)/w/(1.0d0-exp(-(1.44d0*w/temp))) spectra(i,3)=spectra(i,3)+sf(dd,lopt(26),xx,w)*const endif if(lopt(24))then c roa const=stable(i1,5)/w/(1.0d0-exp(-(1.44d0*w/temp))) spectra(i,4)=spectra(i,4)+sf(dd,lopt(26),xx,w)*const endif 2 continue c c generate averaged spectra nav=min(navmx,nav+1) wf=1.0d0/dble(nav) if(lopt(18))then do 441 i=1,npoint 441 avspectra(i,1)=(1.0d0-wf)*avspectra(i,1)+wf*spectra(i,1) endif if(lopt(20))then do 442 i=1,npoint 442 avspectra(i,2)=(1.0d0-wf)*avspectra(i,2)+wf*spectra(i,2) endif if(lopt(22))then do 443 i=1,npoint 443 avspectra(i,3)=(1.0d0-wf)*avspectra(i,3)+wf*spectra(i,3) endif if(lopt(24))then do 444 i=1,npoint 444 avspectra(i,4)=(1.0d0-wf)*avspectra(i,4)+wf*spectra(i,4) endif c if((istep/incr)*incr.eq.istep)then write(s6,6000)istep 6000 format(i6) if(lopt(19))then call wrspec2(np0,s6//'d.prn',npoint,wsmin,dw,spectra,1) call wrspec2(np0,s6//'aved.prn',npoint,wsmin,dw,avspectra,1) endif if(lopt(21))then call wrspec2(np0,s6//'r.prn',npoint,wsmin,dw,spectra,2) call wrspec2(np0,s6//'aver.prn',npoint,wsmin,dw,avspectra,2) endif if(lopt(23))then call wrspec2(np0,s6//'ram.prn',npoint,wsmin,dw,spectra,3) call wrspec2(np0,s6//'averam.prn',npoint,wsmin,dw,avspectra,3) endif if(lopt(25))then call wrspec2(np0,s6//'roa.prn',npoint,wsmin,dw,spectra,4) call wrspec2(np0,s6//'averoa.prn',npoint,wsmin,dw,avspectra,4) endif write(6,*)'spectra written' endif nss=2 if((istep/incr)*incr.eq.istep)write(67,6020)istep 6020 format(i8,$) c ------------------------------------------------ sumsp=0.0d0 do 4 i=1,nexp ik=iskind(i) n11=nsp(i) do 5 i1=1,n11 s1x(i1)=espx(i1,i) 5 s1y(i1)=espy(i1,i) n22=npoint if(lopt(27))then do 6 i1=1,n22 s2x(i1)=wsmin+dw*dble(i1-1) 6 s2y(i1)=avspectra(i1,ik) else do 14 i1=1,n22 s2x(i1)=wsmin+dw*dble(i1-1) 14 s2y(i1)=spectra(i1,ik) endif sum=specomp(xmin(i),xmax(i)) sumsp=sumsp+sum if((istep/incr)*incr.eq.istep)then write(6,6019)i,ik,sum 6019 format(i2,' vs ',i2,': ',g14.5) write(67,6021)sum 6021 format(g14.5,$) endif 4 continue c ------------------------------------------------ c c gradient update: c c shift vector do 7 ia=1,n rsh(3*(ia-1)+1)=x(ia)-xold(ia) rsh(3*(ia-1)+2)=y(ia)-yold(ia) 7 rsh(3*(ia-1)+3)=z(ia)-zold(ia) anorm=0.0d0 do 8 i=1,3*n 8 anorm=anorm+rsh(i)**2 if(anorm.gt.1.0d-9)then anorm=1.0d0/dsqrt(anorm) do 9 i=1,3*n 9 rsh(i)=rsh(i)*anorm endif c c projection of the old gradient onto the new shift direction: g0i=0.0d0 do 10 i=1,3*n 10 g0i=g0i+rsh(i)*gradsp(i) c c new projection: gip1=sumsp-sumspold c c product of the old and new shift directions: gamma=0.0d0 do 12 i=1,n 12 gamma=gamma+rsh(i)*rshold(i) c fac=1.0d0-gamma**2 if(dabs(fac).gt.1.0d-9)fac=(gip1-g0i)/fac gnorm=0.0d0 do 13 i=1,n gradsp(i)=gradsp(i)+fac*(rsh(i)-gamma*rshold(i)) 13 gnorm=gnorm+gradsp(i)**2 gnorm=dsqrt(gnorm) if((istep/incr)*incr.eq.istep)write(67,6021)gnorm c do 11 i=1,n 11 rshold(i)=rsh(i) sumspold=sumsp if((istep/incr)*incr.eq.istep)write(67,*) return end c ==================================================================== subroutine wrspec2(np0,sn,nw,wmin,dw,spec,ic) implicit none character*(*) sn integer*4 nw,iw,ic,np0 real*8 wmin,dw,spec(np0,4),w open(8,file=sn) do 1 iw=1,nw w=wmin+dble(iw-1)*dw 1 write(8,805)w,spec(iw,ic) 805 format(f12.2,f25.12) close(8) return end function sf(d,g,x,x0) implicit none logical g real*8 pi,spi,sf,d,x,x0,dd pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.100.0d0)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif return end subroutine INPTRY(XM,MX,N,kt,r,iso,riso,ISOT,wexp,nm,CHAR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MENDELEV=89) DIMENSION XM(MX),CHAR(*),kt(*),r(3,MX),iso(*),riso(*),wexp(*) CHARACTER*4 CODE dimension amas(MENDELEV) data amas/1.007825,4.003, 2 6.941, 9.012, 10.810,12.000,14.003074,15.994915,18.9984,20.179, 3 22.990,24.305, 26.981,28.086,30.9737,31.972,34.968852,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ logical lopt common/opts/lopt(300) C do 1 i=1,N CHAR(i)=0.0d0 1 XM(i)=amas(kt(i)) nm=3*N do 3 i=1,nm 3 wexp(i)=0.0d0 do 2 i=1,ISOT 2 XM(iso(i))=riso(i) c c Substitute all acidic Hs by Ds if(lopt(3))then ND=0 do 8 i=1,N a=r(1,i) y=r(2,i) z=r(3,i) if(kt(i).eq.1)then do 31 j=1,N it=kt(j) dd=sqrt((a-r(1,j))**2+(y-r(2,j))**2+(z-r(3,j))**2) if(it.eq.7.or.it.eq.8.or.it.eq.9.or.it.eq.16.or.it.eq.17.or. 1 it.eq.35.or.it.eq.53)then if(dd.lt.1.23d0)then ND=ND+1 XM(i)=2.014000d0 endif endif 31 continue endif 8 continue write(6,7009)ND 7009 format(i5,' acidic hydrogens have been substituted by deuteria') endif c c generate FTRY.INP if(lopt(2))then CODE='INTY' OPEN(7,FILE='FTRY.INP') WRITE(7,1111)CODE,0,1,0,1,0,1,0 1111 FORMAT(A4,7I4) WRITE(7,333)0 333 FORMAT(20I4) WRITE(7,444)(CHAR(i),I=1,N) 444 format(6f12.6) WRITE(7,333)1 WRITE(7,444)(XM(I),I=1,N) WRITE(7,*)3*N WRITE(7,444)(wexp(i),I=1,3*N) CLOSE(7) write(6,*)'FTRY.INP created' ENDIF c read FTRY.INP if(lopt(4))then OPEN(7,FILE='FTRY.INP') read(7,*) read(7,*) read(7,444)(CHAR(i),I=1,N) read(7,*) read(7,444)(XM(I),I=1,N) read(7,*)nm read(7,444)(wexp(I),I=1,nm) CLOSE(7) write(6,*)'FTRY.INP read' ENDIF return end subroutine new4(NOAT,FCAR,CHAR,ZIM,wexp,nm,SMAT,FRQ,R) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50,MX3=3*MX) DIMENSION ZIM(*),CHAR(*),FCAR(MX3,MX3),wexp(*),FRQ(*), 1SMAT(MX3,MX3),R(3,MX) COMMON/ITERA/ZM(MX3),Q(MX3),ZNU(MX3),AS(3,MX3),NA logical lopt common/opts/lopt(300) C NA=3*NOAT do 130 i=1,3 do 130 j=1,NA 130 AS(i,j)=0.0d0 C CHATOT=0.0d0 DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHAR(I) C BM=0.0d0 DO 102 IA=1,NOAT 102 BM=BM+ZIM(IA) DO 220 I=1,NOAT DO 220 J=1,3 220 ZM(3*I-3+J)=1/SQRT(ZIM(I)) C DO 211 I=1,NA 211 ZNU(I)=0.0d0 DO 212 I=1,nm 212 ZNU(I)=wexp(I) C C IF(lopt(5))CALL PROJECTF4(FCAR,MX3,R) C C ie=0 ADUM=SUMk(0,ie,FCAR,SMAT) C Calculate the true mass-weighted S-matrix: DO 210 IA=1,NA T=ZM(IA) DO 210 IM=1,NA 210 SMAT(IA,IM)=SMAT(IA,IM)*T DO 213 IM=1,NA 213 FRQ(IM)=AS(1,IM) if(lopt(6))CALL SMATOUT(SMAT,AS,NA) if(lopt(7))CALL FPC(Q,CHAR,SMAT,NOAT,3*NOAT,R) return END C SUBROUTINE SMATOUT(S,AS,NAT3) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50,MX3=3*MX) DIMENSION S(MX3,NAT3),AS(3,NAT3) IRANGE=1 JRANGE=NAT3 NI=JRANGE-IRANGE+1 c c look at the phase and make the biggest element always positive: DO 101 J=IRANGE,JRANGE amx=S(1,J) DO 102 I=2,NAT3 102 if(abs(S(I,J)).gt.abs(amx))amx=S(I,J) if(amx.lt.0)then do 103 I=1,NAT3 103 S(I,J)=-S(I,J) endif 101 continue c OPEN(34,FILE='F.INP') WRITE(34,10)NI,NAT3,NAT3/3 10 FORMAT(3I5) OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*)NAT DO 3 I=1,NAT READ(35,*)IAT,X,Y,Z 3 WRITE(34,11)IAT,X,Y,Z 11 FORMAT(I4,3F12.6) CLOSE(35) WRITE(34,14) 14 FORMAT(' Atom Mode X-disp. Y-disp. Z-disp.') DO 1 I=1,NAT3/3 DO 1 J=IRANGE,JRANGE IU=3*(I-1) K=NAT3-J+1 1 WRITE(34,15)I,K,S(IU+1,J),S(IU+2,J),S(IU+3,J) 15 FORMAT(2I5,3F11.6) WRITE(34,11)NI WRITE(34,13)(AS(1,I),I=IRANGE,JRANGE) 13 format(6F11.3) CLOSE(34) WRITE(*,*)' S-matrix written out ...' RETURN END FUNCTION SUMk(iwr,iee,FCAR,SMAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=50,MX3=3*MX) dimension FCAR(MX3,MX3),EWORK(MX3),SMAT(MX3,MX3),FNEW(MX3,MX3) COMMON/ITERA/ 1ZM(MX3), 2Q(MX3),ZNU(MX3), 2AS(3,MX3),NA c wtol=0.01d0 c1=1302.828d0 CALL FASTO(NA,FCAR,FNEW,ZM) IMOD=2 CALL TRED12(MX3,NA,FNEW,SMAT,Q,IMOD,IERR,EWORK) IF(IERR.NE.0)call report('diagonalization error') c RMSW=0.0D0 ie=0 DO 310 I=1,NA T=Q(I) SIGN=1.0d0 C Make imaginary frequencies negative: IF(T.LT.0.0d0)SIGN=-1.0d0 T=SIGN*SQRT(ABS(T))*c1 Q(I)=T TP=ZNU(I) DEL=T-TP if(ABS(TP).le.wtol)then DEL=0.0d0 else ie=ie+1 endif RMSW=RMSW+DEL**2 AS(1,I)=T AS(2,I)=TP 310 AS(3,I)=T-TP iee=ie if(ie.gt.0)then RMSW=DSQRT(RMSW/DFLOAT(ie)) else if(NA.gt.0)RMSW=DSQRT(RMSW/DFLOAT(NA)) endif SO1=0.0D0 SO2=0.0D0 DO 5151 J=1,NA SO1=SO1+AS(1,J)*AS(2,J) 5151 SO2=SO2+AS(1,J)*AS(1,J) AK=SO1/SO2 C if(iwr.gt.0)WRITE(*,5555) 5555 FORMAT(1X,'************ NEW FREQUENCIES ******************') DO 311 J=NA,1,-1 a3=0.0d0 if(AS(2,NA-J+1).gt.wtol)a3=AS(3,NA-J+1) 311 if(iwr.gt.0)WRITE(*,6666)J,NA-J+1,(AS(I,NA-J+1),I=1,2),a3 6666 FORMAT(1X, 1 I5,' (',I4,') ...... ',F8.2,' -',F8.2,'/',F8.2) if(iwr.gt.0)WRITE(*,5551) 5551 FORMAT(1X,'************************************************') if(iwr.gt.0)WRITE(6,6667)ie,RMSW,AK 6667 FORMAT(i6,'experimental frequencies found',/, 1 ' >>>>><< >> ',F8.2,' <<<<<< ',F8.6) c SUMk=RMSW RETURN END SUBROUTINE FPC(V,CHAR,S,NOAT,NA3,X) C PETR BOUR: LAST UPDATE 4-28-92 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50,MX3=3*MX) DIMENSION V(MX3), S(MX3,NA3), X(3,MX), 1 CHAR(MX), DS(3), SCRT(3) DCOR = 388891.4D0 RCOR = 12221.7D0 NA=3*NOAT NVIBS=NA C INCLUDE TRANSLATIONS AND ROTATIONS, TOO OPEN(13,FILE='FPC.TAB') WRITE(13,50) 50 FORMAT(3X,'FREQ',8X,'DIP. STRN.*E39',6X,'ROT. STRN.*E45'/ 1 3X,'CM-1',2X,2(6X,' (FR CM)**2 '),/,'---------------') DO 42 I=1,NVIBS DO 784 IKY=1,3 DS(IKY)=0.0D0 SCRT(IKY)=0.0D0 DO 784 JR=1,NOAT DS(IKY)=DS(IKY)+CHAR(JR)*S(3*(JR-1)+IKY,I) DO 784 IB=1,3 DO 784 IG=1,3 784 SCRT(IKY)=SCRT(IKY)+EPS(IKY,IB,IG)*X(IB,JR)*S(3*(JR-1)+IG,I) 1 *CHAR(JR) DT =DS(1)*DS (1)+DS(2)*DS (2)+DS(3)*DS (3) ROT=DS(1)*SCRT(1)+DS(2)*SCRT(2)+DS(3)*SCRT(3) ROTC =+ROT*RCOR DIPOL=0.0D0 IF (ABS(V(I)).GT.0) DIPOL = DT*DCOR/ABS(V(I)) 42 WRITE(13,85) I,V(I),DIPOL,ROTC,DT,ROT 85 FORMAT(I4,F8.2,' ',F13.6,' ',f13.6,' ',2G14.5) WRITE(13,86) 86 FORMAT('---------------------------------------------') CLOSE(13) RETURN END C SUBROUTINE FASTO(NA,FCAR,FNEW,ZM) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=50,MX3=3*MX) DIMENSION FCAR(MX3,NA),FNEW(MX3,NA),ZM(NA) DO 214 I=1,NA UI=ZM(I) DO 214 J=I,NA FNEW(I,J)=FCAR(I,J)*ZM(J)*UI 214 FNEW(J,I)=FNEW(I,J) RETURN END C C SUBROUTINE PROJECTF4(F,N0,R) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50,MX3=3*MX) DIMENSION F(N0,N0),TEM(MX3,MX3),A(MX3,MX3),R(3,MX) C C Using the common here to save memory for A and TEM COMMON/ITERA/ZM(MX3),Q(MX3),ZNU(MX3),AS(3,MX3),NA N=NA CALL DOMAA(A,R) 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 DOMAA(A,C) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=50,MX3=3*MX) DIMENSION C(3,MX),TEM(MX3,MX3),A(MX3,MX3) COMMON/ITERA/ZM(MX3),Q(MX3),ZNU(MX3),AS(3,MX3),NA N=NA AMACH=0.00000000000001d0 NAT=N/3 DO 12 I=1,N DO 12 J=1,6 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 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 ORT1(A,MX3,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE DO 2 I=1,N DO 2 J=1,N 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,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END C ============================================================ SUBROUTINE ORT1(A,NDIM,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NDIM,N6) AMACH=0.00000000000001d0 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 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 new5(PP,II,NOAT,R,SV,FRQ,nm,stable) implicit none integer*4 MX,A,B,NOAT,I,J,L,III,IR, 1nm PARAMETER (MX=50) REAL*8 PP(MX,3,3),II(MX,3,3),SV(3*MX,3*MX), 1 R(3,MX),FRQ(*), 2 ELD(3),MDT(3), 3 ME(3),MM(3),AI(3,3),TI(3), 4pi,FD,FC,Z0,WLIM,W,WR,ABSW, 5RDO, 6SVL,sumd,stable(3*MX,5) logical lopt common/opts/lopt(300) c c write(6,*)NOAT,nm pi=4.0d0*atan(1.0d0) FD =SQRT(388.9219d0) FC=131.14D-8 Z0=0.0d0 WLIM=0.000001d0 if(lopt(12))then open(17,file='F.INP') 101 read(17,*,end=102,err=102) goto 101 102 continue endif if(lopt(9))OPEN(18,FILE='DOG.TAB') if(lopt(10))OPEN(19,FILE='DIP.MOM') if(lopt(9))WRITE (18,501) if(lopt(10))CALL INERTIA(R,AI,TI,NOAT) if(lopt(10))WRITE (19,5019)(TI(I),I=1,3) 5019 FORMAT( 1' POLARIZATIONS ACCORDING TO THE INERTIA TENSOR AXES',/, 2' ( ',3F15.4, ')',/, 3' MODE ELECTRIC ', 4' MAGNETIC [%] ',/,80(1H-),/) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2/,'-----------------------------------------------------------') C DO 41 I=1,nm W=FRQ(I) ABSW=DABS(W) WR=DSQRT(ABSW) DO 56 B=1,3 ELD(B)=Z0 MDT(B)=Z0 DO 56 L=1,NOAT DO 56 A=1,3 SVL=SV(3*(L-1)+A,I) ELD(B) = ELD(B) + PP(L,A,B)*SVL 56 MDT(B) = MDT(B) + II(L,A,B)*SVL c DO 51 B=1,3 IF(ABSW.GT.WLIM)ELD(B)=ELD(B)/WR*FD MDT(B)=MDT(B)*FC*WR 51 continue C SUMD = ELD(1)*ELD(1)+ELD(2)*ELD(2)+ELD(3)*ELD(3) RDO = MDT(1)*ELD(1)+MDT(2)*ELD(2)+MDT(3)*ELD(3) C DO 39 III=1,3 ME(III)=Z0 MM(III)=Z0 DO 39 J=1,3 MM(III)=MM(III)+AI(J,III)*MDT(J) 39 ME(III)=ME(III)+AI(J,III)*ELD(J) C CALL NORM1(ME) CALL NORM1(MM) C DO 4 III=1,3 MM(III)=MM(III)*MM(III)*100.0d0 4 ME(III)=ME(III)*ME(III)*100.0d0 C if(lopt(10))WRITE(19,1920)I,W,(ME(IR),IR=1,3),(MM(IR),IR=1,3) 1920 FORMAT(I4,F7.1,4(3F10.0,' ')) if(lopt(12))WRITE(17,522)(ELD(IR),IR=1,3),(MDT(IR),IR=1,3) 522 FORMAT(12G13.5) IF(ABSW.LT.WLIM)SUMD=Z0 stable(I,1)=W stable(I,2)=SUMD stable(I,3)=RDO 41 if(lopt(9))WRITE (18,523) I,W,SUMD,RDO,RDO,SQRT(PI+PI)*W*SUMD 523 FORMAT(I4,5G15.6) C if(lopt(9))WRITE(18,1819) if(lopt(10))WRITE(19,1819) 1819 FORMAT('-------------------------------------------------') if(lopt(12))CLOSE(17) if(lopt(9))CLOSE(18) if(lopt(10))CLOSE(19) return END SUBROUTINE INERTIA(R,A,TI,N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (MX=50) DIMENSION R(3,MX),A(3,3),T(3,3),TI(3),CM(3),E(3) DO 2 I=1,3 CM(I)=0.0d0 DO 3 J=1,N 3 CM(I)=CM(I)+R(I,J) 2 CM(I)=CM(I)/N DO 1 I =1,3 DO 1 J =1,3 T(I,J)=0.0d0 DO 1 IS=1,N 1 T(I,J)=T(I,J)+(R(I,IS)-CM(I))*(R(J,IS)-CM(J)) CALL TRED12(3,3,T,A,TI,2,IERR,E) RETURN C DEL(I,J)=SUM(OVER JS) OF A(I,JS)*A(J,JS) C DIAGONAL TD(I,J) = A(IS,I)*T(IS,JS)*A(JS,j) END SUBROUTINE NORM1(A) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3) AN=DSQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3)) IF(ABS(AN).LT.1.0d-30)RETURN A(1)=A(1)/AN A(2)=A(2)/AN A(3)=A(3)/AN RETURN END subroutine new6(S,wmin,wmax,ALPHA,G,A,NAT,EXCA,nm,FRQ,r,stable) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=50,MX3=3*MX) DIMENSION S(MX3,MX3),ALPHA(MX3,3,3),A(MX3,3,3,3),G(MX3,3,3), 3R(3,MX3/3),FRQ(*),stable(3*MX,5) LOGICAL LPOLAR logical lopt common/opts/lopt(300) C LPOLAR=lopt(14) IF(LPOLAR)THEN c set A, G to zero: CALL TrTT(MX3,ALPHA,A,G,NAT,3,R) c make common A, G from alpha: CALL TrTT(MX3,ALPHA,A,G,NAT,2,R) write(6,*)' A, G made from alpha !!!' ENDIF CALL TRTEN(MX3,S,3*NAT,ALPHA,A,G) CALL DORAMAN(MX3,nm,ALPHA,G,A,EXCA,wmin,wmax,FRQ,stable) return END C ========================================== C ============================================================== SUBROUTINE DORAMAN(NS,N,ALPHA,G,A,EXCA,wmin,wmax,E,stable) c c supposedly G/W is supplied c IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) DIMENSION ALPHA(NS,3,3),G(NS,3,3),A(NS,3,3,3),E(*), 1stable(NS,5) logical lopt common/opts/lopt(300) c c if(lopt(16))OPEN(4,FILE='ROA.TAB') if(lopt(16))WRITE(4,3304) 3304 FORMAT( 1'MODE FREQ Ramx (180) Ramy Ramtot CID90', 2' CID-X CID180 DEP ROA180 RamanDCPI',/, 3' cm-1 A^4/AMU x10^4', 4' x10^4 x10^4 A^5/AMUx1e4 au.10^4') if(lopt(16))write(4,3002) 3002 FORMAT(80(1H-)) c DO 1 IQ=1,N ECM=E(IQ) IF(ABS(ECM).GT.0.01d0)THEN if(wmin.eq.0.0d0.or.ECM.gt.wmin)then if(wmax.eq.0.0d0.or.ECM.lt.wmax)then SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 SPAL=0.0d0 DO 3 I=1,3 SPAL=SPAL+ALPHA(IQ,I,I) DO 3 J=1,3 SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J) SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J) SAG0=SAG0+ALPHA(IQ,I,I)* G(IQ,J,J) SAG1=SAG1+ALPHA(IQ,I,J)* G(IQ,I,J) DO 3 K=1,3 DO 3 L=1,3 3 SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0 c c IL+IR, backscattering, DCP_I c SDCP180=3.0d0*SAL1-SAL0 c c IL+/-IR, backscattering S180=7.0d0*SAL1+SAL0 D180=8.0d0*(3.0d0*SAG1-SAG0+SA1) D=D180 IF(S180.gt.1.0d-9)D=D/S180 c c IL+/-IR, 90o, x S90X=7.0d0*SAL1+SAL0 D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1) DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c c IL+/-IR, 90o, z S90Z=2.0d0*(3.0d0*SAL1-SAL0) D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1) DZ=D90Z IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z c c polarized backscattering components (YDY+YDX=S180): YDY=3.0d0*SAL1-SAL0 YDX=4.0d0*SAL1+2.0d0*SAL0 P1=YDY/YDX c c Try to get Gaussian units: AMU=1822.0d0 BOHR=0.529177d0 SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) YDY=6.0d0*YDY*(AMU*BOHR**4) YDX=6.0d0*YDX*(AMU*BOHR**4) c gpisvejc=980.0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA c c c beta(A)^2=Delta2 in Gaussian: betaa=SA1*3.0d0/2.0d0*gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc c c DCP 180 c ram3=24.0d0*beta2 roa3=8.0d0*(12.0d0*betag+4.0d0*betaa) c c unknown: S90X=S90X*(AMU*BOHR**4) D90X=D90x*gpisvejc DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c c ICP(90)=90Z roa2=8.0d0*(3.0d0*betag-betaa) ram2=12.0d0*beta2 CID2=0.0d0 if(ram2.gt.1.0d-9)CID2=roa2/ram2 c c ICP(180): c 8/c*(12*beta(G')^2 + 4*beta(A)^2): roa1=96.0d0*(betag+betaa/3.0d0) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) CID1=0.0d0 if(ram1.gt.1.0d-9)CID1=roa1/ram1 stable(IQ,1)=ECM stable(IQ,4)=ram1 stable(IQ,5)=roa1 if(lopt(16))WRITE(4,3300)IQ,ECM,YDX,YDY, 1 ram1,CID2,DX,CID1,P1,roa1,roa3 3300 FORMAT(I4,f9.2,6f9.2,f6.3,2f11.2) c endif endif ENDIF 1 CONTINUE if(lopt(16))WRITE(4,3002) if(lopt(16))CLOSE(4) RETURN END SUBROUTINE TRTEN(NS,S,N,ALPHA,A,G) IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) PARAMETER (MX=50,MX3=3*MX) DIMENSION S(NS,N),A(NS,3,3,3), 1ALPHA(NS,3,3),G(NS,3,3),V(MX3) logical lopt common/opts/lopt(300) SFAC=0.0234280d0 DO 1 I=1,3 DO 1 J=1,3 DO 3 JQ=1,N 3 V(JQ)=ALPHA(JQ,I,J) DO 1 IQ=1,N ALPHAS=0.0d0 DO 2 IX=1,N 2 ALPHAS=ALPHAS+S(IX,IQ)*V(IX)*SFAC 1 ALPHA(IQ,I,J)=ALPHAS DO 11 I=1,3 DO 11 J=1,3 DO 31 JQ=1,N 31 V(JQ)=G(JQ,I,J) DO 11 IQ=1,N GS=0.0d0 DO 21 IX=1,N 21 GS=GS+S(IX,IQ)*V(IX)*SFAC 11 G(IQ,I,J)=GS DO 12 I=1,3 DO 12 J=1,3 DO 12 K=1,3 DO 32 JQ=1,N 32 V(JQ)=A(JQ,I,J,K) DO 12 IQ=1,N AS=0.0D0 DO 22 IX=1,N 22 AS=AS+S(IX,IQ)*V(IX)*SFAC 12 A(IQ,I,J,K)=AS if(lopt(17))WRITE(6,3001) 3001 FORMAT(1X,'TENSORS HAVE BEEN TRANSFORMED TO Q COORDINATES',/, 11X,'ALPHA DERIVATIVES LISTED:') DO 4 IQ=1,N SP=ALPHA(IQ,1,1)+ALPHA(IQ,2,2)+ALPHA(IQ,3,3) 4 if(lopt(17))WRITE(6,3000)IQ,SP,((ALPHA(IQ,I,J),I=1,3),J=1,3) 3000 FORMAT(1X,I4,' MODE; Trace: ',f15.8,/,3(3F15.8,/),/) if(lopt(17))WRITE(6,3002) 3002 FORMAT(1X, 11X,'G (magn index -->)') DO 5 IQ=1,N SP=G(IQ,1,1)+G(IQ,2,2)+G(IQ,3,3) 5 if(lopt(17))WRITE(6,3000)IQ,SP,((G(IQ,J,I),I=1,3),J=1,3) if(lopt(17))WRITE(6,3003) 3003 FORMAT(1X, 11X,'A DERIVATIVES LISTED:') DO 6 IQ=1,N 6 if(lopt(17))WRITE(6,3004)IQ,(((A(IQ,I,J,K),I=1,3),J=1,3),K=1,3) 3004 FORMAT(1X,I4,' MODE',/,3(9F8.5,/),/) RETURN END c ==================== SUBROUTINE TrTT(N,ALPHA,A,G,NAT,ICO,R) c C Transforms A and G to local origins (ICO=1) or common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in AU C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=600) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,MX) C C transfer into atomic coordinates: DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I,J) C IF(ICO.EQ.1)S=1.0d0 IF(ICO.EQ.2)S=-1.0d0 IF(ICO.EQ.3)S=0.0d0 C DO 1 I=1,3 DO 1 J=1,3 DO 1 L=1,NAT DO 1 IE=1,3 II=3*(L-1)+IE DO 1 K=1,3 DO 1 M=1,3 c G-second index-magnetic 1 G(II,I,J)=G(II,I,J)*S*S+S*0.5d0*EPS(J,K,M)*X(K,L)*ALPHA(II,I,M) IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' C DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 II=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(II,IA,ID) ENDIF A(II,IA,IB,IC)=A(II,IA,IB,IC)+SUM*S 1-1.5d0*S*(X(IB,L)*ALPHA(II,IA,IC)+X(IC,L)*ALPHA(II,IA,IB)) IF(ICO.EQ.3)A(II,IA,IB,IC)=0.0d0 3 CONTINUE c IF(ICO.EQ.1)WRITE(6,*)' A,G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' A,G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A,G set to zero' RETURN END function specomp(xmin,xmax) implicit none integer*4 ns0,np0 parameter (ns0=10,np0=10000) character*80 sn(ns0), s80 character*10 s10 integer*4 is,i,j,js,IERR,ii real*8 as(ns0),ks(ns0),tem,bs(ns0),norm1,normj,so(ns0,ns0), 1vo(ns0),e(ns0,2*ns0),ai(ns0,ns0),TOL,coe,xmin,xmax,specomp integer n1,n2,ITER,ns,jsc real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b COMMON/MATVARs/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,a,k,b, 1ns,jsc real*8 wmin,wmax,aa,kk,bb,al,kin,bin integer*4 npar,iwr,icol,iinv,ino common/optsp/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino c wmin=xmin wmax=xmax call norms(n1,s1x,s1y,norm1,wmin,wmax,ino) call norms(n2,s2x,s2y,normj,wmin,wmax,ino) call DOMATRIXsp specomp=RMS c write(6,6002)js,a,k,b,RMS,ITER c002 format(i4,4f8.3,i5) return end subroutine norms(n,sx,sy,an,wmin,wmax,ino) integer*4 n,i,ino real*8 an,sx(*),sy(*),wmin,wmax an=0.0d0 if(ino.eq.0)then do 1 i=1,n-1 1 if(sx(i).ge.wmin.and.sx(i).le.wmax) 1 an=an+sy(i)**2*abs(sx(i+1)-sx(i)) an=dsqrt(an) else do 11 i=1,n-1 11 if(sx(i).ge.wmin.and.sx(i).le.wmax) 1 an=an+abs(sy(i))*abs(sx(i+1)-sx(i)) endif if(an.gt.1.0d-24)then do 2 i=1,n 2 sy(i)=sy(i)/an endif return end SUBROUTINE DOMATRIXsp C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT none integer*4 npar0,ITMAX,np0,ns0 parameter (npar0=3,ITMAX=500,np0=10000,ns0=10) real*8 P(4,npar0),Y(4),PBAR(npar0),PR(npar0),PRR(npar0), 1FUs,FTOL,YPR,YPRR,ANG(npar0) parameter (FTOL=1.0d-4) integer*4 ILO,INHI,IHI,I,J integer n1,n2,ITER,ns,jsc real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b COMMON/MATVARs/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,a,k,b, 1ns,jsc real*8 wmin,wmax,aa,kk,bb,al,kin,bin integer*4 npar,iwr,icol,iinv,ino common/optsp/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino c C STARTING VALUES FOR THE ITERATION: c P(,1)..a, P(,2)..k, P(,3)..b if(npar.gt.npar0)call report('too many parameters') P(1,1)=1.0d0 P(1,2)=1.0d0 P(1,3)=0.0d0 P(2,1)=1.0d0 P(2,2)=1.1d0 P(2,3)=1.0d1 P(3,1)=1.1d0 P(3,2)=1.0d0 P(3,3)=5.0d0 P(4,1)=1.1d0 P(4,2)=1.1d0 P(4,3)=2.0d0 do 1 i=1,4 ANG(1)=P(i,1) ANG(2)=P(i,2) ANG(3)=P(i,3) 1 Y(i)=FUs(ANG) c ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(6,*)' Iteration has not converged !' RETURN ENDIF ITER=ITER+1 DO 55 I=1,npar 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,npar 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,npar PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FUs(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,npar 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FUs(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,npar 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,npar 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,npar 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,npar 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FUs(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,npar 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,npar PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FUs(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,npar 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FUs(ANG) implicit none INTEGER*4 i,np0,ns0 parameter (np0=10000,ns0=10) REAL*8 FUs,ANG(*),dd,w,kw integer*4 jc,j integer n1,n2,ITER,ns,jsc real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b COMMON/MATVARs/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,a,k,b, 1ns,jsc real*8 wmin,wmax,aa,kk,bb,al,kin,bin integer*4 npar,iwr,icol,iinv,ino common/optsp/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino c a=ANG(1) k=ANG(2) b=ANG(3) RMS=0.0d0 do 1 i=1,n1 w=s1x(i) if(w.ge.wmin.and.w.le.wmax)then kw=k*w if(npar.gt.2)kw=kw+b c do 23 jc=1,n2-1 23 if((s2x(jc).le.kw.and.s2x(jc+1).ge.kw).or. 1 (s2x(jc).ge.kw.and.s2x(jc+1).le.kw))goto 99 goto 1 99 RMS=RMS+dabs(a*s2y(jc)-s1y(i)) c endif 1 continue RMS=RMS+kk*(kin-k)**2+aa*(1.0d0-a)**2 if(npar.gt.2)RMS=RMS+bb*(b-bin)**2 FUs=RMS RETURN END subroutine ioparspecomp character*10 s10 real*8 wmin,wmax,aa,kk,bb,al,kin,bin integer*4 npar,iwr,icol,iinv,ino common/optsp/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino logical lex npar=3 kk=1.0d3 aa=1.0d1 bb=4.0d-3 iwr=0 icol=0 al=0.001d0 iinv=0 kin=1.0d0 bin=0.0d0 ino=0 inquire(file='SP.OPT',exist=lex) if(lex)then open(9,file='SP.OPT') 1 read(9,90,end=99,err=99)s10 90 format(a10) if(s10(1:4).eq.'alph')read(9,*)al if(s10(1:2).eq.'kk')read(9,*)kk if(s10(1:2).eq.'aa')read(9,*)aa if(s10(1:2).eq.'bb')read(9,*)bb if(s10(1:3).eq.'bin')read(9,*)bin if(s10(1:3).eq.'kin')read(9,*)kin if(s10(1:4).eq.'npar')read(9,*)npar if(s10(1:4).eq.'iinv')read(9,*)iinv if(s10(1:3).eq.'ino')read(9,*)ino c c print control: if(s10(1:3).eq.'iwr')read(9,*)iwr c c collective spectra fitting: if(s10(1:4).eq.'icol')read(9,*)icol goto 1 99 close(9) endif write(6,6000)aa,kk,bb,kin,bin,npar 6000 format( 1'parameter restrictions:',/, 2' aa (intenzity) : ',f8.1,/, 3' kk (x-scale multiplier): ',f8.1,/, 4' bb (x-scale shift) : ',f8.4,/,/, 7' ki (syst. x-multiplier): ',f8.1,/, 8' bi (syst. x-shift) : ',f8.1,/, 5' npar : ',i8,/) return end