program dynamiM implicit none integer*4 maxval,nw0,nn0,mm,nmol,npair,nnames,iend, 1istart,K0 parameter (nw0=2000,maxval=8,K0=5) c nw0 ... maximum number of spectral points c maxval maximal valence numbers (number of bonds from one atom) c K0 maximum number of named kinds of molecules integer*4 ia,iw,iprint,n,nw,ix,nat,NAT1,LP,istep,i,iind,ip,ixp,j, 1jp,jx,k,kp,l,NS,n1,n2,okind,NO,NATS,IWG,nsa,iap,ii, 1m1,m2,natp,mp,nib,ibx,iby,ibz,ibxp,ibyp, 1ibzp,mxb integer*4,allocatable:: i12(:,:),iii(:,:),atomic(:),isr(:), 1al(:,:),alp(:,:),IIM(:,:,:),natm(:),natmp(:),imk(:) character*10 snam(K0) character*1 lg logical ltti,ltei,lexp,lper,lpair,ldress,lp27,lp28, 1LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT,lexca, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,lcct,lproject,lani, 1lconv,lp29,lcell,lincr c wsu 1 el _x sin c wsu 2 el _x cos c wsu 3 el _y sin c wsu 4 el _y cos c wsu 5 el _z sin c wsu 6 el _z cos c wsu 7 mag_x sin c wsu 8 mag_x cos c wsu 9 mag_y sin c wsu10 mag_y cos c wsu11 mag_z sin c wsu12 mag_z cos c for all molecules and spectral points c index = i12 + (iw-1)*12 +(mol-1)*nw0*12 c current (x,y,z), previous (xm,ym,zm) and previous c previous(xmm,ymm,zmm) atomic coordinates c rm are the molecular coordinates real*8 dt,spd,smd,POLI(3,3),s3,qx,qy,qz,rdx,rdy,rdz, 5spec(nw0),stem(nw0),t,w,dbb,ramconst,dw,wmax,tp,sol,wmin,res,sc, 6fac,FBOHR,dtau,bohr,pi,x1,y1,z1,x2,y2,z2,kelvin,vx,vy,vz,rd2,rd5, 1dx,dy,dz,cf2,p1(9),p2(9),DM,POLM(K0,9),fac1,dxp,dyp,dzp,sa1, 1sal0,sal1,sag0,sag1,ax,ay,az,ac(3),bc(3),cc(3),s180,d180,ydy,ydx, 1tf,s1,s2,rams(4,nw0),aan,abn,acn,ap,bp,cp,apm,bpm,cpm, 1at1(3,3),at2(3,3),pt1(3,3),pt2(3,3),ut(3,3),rki(3),rd1, 1gs1(3,3,3),as1(3,3,3),au1(3,3,3,3),au2(3,3,3,3),as2(3,3,3), 1gs2(3,3,3),CUTOFF,perx,pery,perz,px2,py2,pz2,rd67 real*8,allocatable::x(:),y(:),z(:),rm(:),xx(:),yy(:),zz(:), 1xm(:),ym(:),zm(:),xmm(:),ymm(:),zmm(:),pchgft(:),pol(:), 2P(:,:,:),A(:,:,:),ALPHAI(:,:,:),AAI(:,:,:,:),GI(:,:,:), 4ALPHA1(:,:,:),AA1(:,:,:,:),G1(:,:,:),ALPHA(:,:,:),po2(:,:,:), 4ALP29(:,:),ALPHM(:,:),ALPHAM(:,:,:),GM(:,:,:),AAM(:,:,:), 5APTM(:,:,:),AM(:,:,:),AA(:,:,:),G(:,:,:),APTI(:,:,:),AI(:,:,:), 6u(:,:,:),Xwork(:,:),uu(:),gt(:,:),at(:,:),aat(:,:), 7wsu(:),wsat(:,:),wcat(:,:),wsgt(:,:),wcgt(:,:), 2wsaat(:,:),wcaat(:,:) 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 logical lnnames logical,allocatable:: luse(:) c cctn variables: integer*4,allocatable::IND(:),IB(:,:),ND(:),NDS(:),IS(:,:), 1INDBIG(:,:),IBONDSB(:,:),IBONDS(:,:),NBO(:),NBOB(:) real*8,allocatable::R(:,:),RS(:,:),FFS(:,:),FF(:,:), 1PS(:,:,:),AS(:,:,:),ALPHAS(:,:,:),AAS(:,:,:,:), 1GS(:,:,:),AAT_c(:,:,:,:),GT_c(:,:,:),AM_c(:,:,:) integer*4 IOX c IOX max # of overlaps parameter (IOX=30) real*8,allocatable::asname(:,:,:,:),psname(:,:,:,:), 1aaname(:,:,:,:,:),gname(:,:,:,:),alphaname(:,:,:,:), 1rname(:,:,:),ffname(:,:,:) integer*4,allocatable::NATNAME(:),NDNAME(:,:) write(6,6009) 6009 format(/,/, 1' DYNAMIM: spectra generation from molecular trajectories',/,/, 1' Petr Bour',/, 1' Academy of Sciences, Prague, 2015',/,/, 1' Input: FILE.X MD trajectories',/, 1' <> tensor files',/, 1' DYNAMIM.OPT options',/, 1' Output: .PRN spectra',//) call readopt3(nmol,lpair) nn0=nw0*nmol open(10,file='FILE.X',status='old') read(10,*) read(10,*)n NAT=n allocate(x(n),y(n),z(n),rm(3*n),xx(n),yy(n),zz(n),xm(n),ym(n), 1zm(n),xmm(n),ymm(n),zmm(n),pchgft(n),pol(n),P(n,3,3),A(n,3,3), 2ALPHAI(3*n,3,3),AAI(3*n,3,3,3),GI(3*n,3,3),ALPHA1(3*n,3,3), 3AA1(3*n,3,3,3),G1(3*n,3,3),ALPHA(9,nmol,3*n),po2(n,3,3), 4ALP29(9,3*n),ALPHM(9,3*n),ALPHAM(9,K0,3*n),GM(9,K0,3*n), 6AAM(27,K0, 3*n),APTM(n,9,K0),AM(n,9,K0),AA(27,nmol,3*n), 7G(9,nmol,3*n),APTI(n,3,3),AI(n,3,3),u(3,3,n),Xwork(3,n), 8i12(maxval,n),iii(4,n),atomic(n),isr(n), 1IIM(K0,4,n),luse(n),natm(nmol),natmp(nmol),imk(nmol), 6uu(6*nmol),gt(9,nmol),at(9,nmol),aat(27,nmol)) call rg(n,x,y,z,atomic,i12,0,maxval) close(10) if(lpair)then npair=nmol*(nmol-1)/2 write(6,6001)npair 6001 format(i4,' pairs') else npair=nmol endif allocate(al(npair,n),alp(npair,n)) call readopt1(K0,n,NSA,iprint,nmol,nnames,isr,natm, 1imk,al, 1dt,kelvin,bohr,pi,dbb,perx,pery,perz,tp,sol,dw,wmin,wmax,res, 1ac,bc,cc, 2lper,lcell,lp27,lp28,lp29,lproject,lnnames,lpair,LABS,LVCD,LRAM, 2LROA,ltti,ldress,ltei,lani,lconv,lg,snam,lincr) c px2=perx/2.0d0 py2=pery/2.0d0 pz2=perz/2.0d0 do 314 ia=1,n pchgft(ia)=0.0d0 314 luse(ia)=.false. if(NSA.ne.0)then do 315 ia=1,NSA 315 luse(isr(ia))=.true. else do 316 ia=1,n 316 luse(ia)=.true. endif write(6,*)'Active atoms:' do i=1,n if(luse(i))then write(6,6012)'X' else write(6,6012)'O' 6012 format(a1,$) endif if(mod(i,60).eq.0)write(6,*) enddo write(6,*) c inquire(file='CHARGE.LST',exist=lexca) if(lexca)then write(6,*)'CHARGE.LST found, partial charges read from it' open(99,file='CHARGE.LST') read(99,*) do 19 ia=1,n 19 read(99,*)pchgft(ia) close(99) endif inquire(file='POLAR.LST',exist=lexp) if(lexp)then write(6,*)'POLAR.LST found, atomic polarizabilities read from it' open(99,file='POLAR.LST') read(99,*) do 18 ia=1,n 18 read(99,*)pol(ia) close(99) if(lpair)then ip=0 do 181 m1=1,nmol do 181 m2=m1+1,nmol ip=ip+1 natp=0 do 182 ia=1,natm(m1) natp=natp+1 182 alp(ip,natp)=al(m1,ia) do 183 ia=1,natm(m2) natp=natp+1 183 alp(ip,natp)=al(m2,ia) 181 natmp(ip)=natp do 184 ip=1,npair natm(ip)=natmp(ip) do 184 ia=1,natm(ip) 184 al(ip,ia)=alp(ip,ia) endif endif if(ltti)then if(lnnames)then do 186 mm=1,nnames do 1851 istart=1,len(snam(mm)) 1851 if(snam(mm)(istart:istart).ne.' ')goto 1852 1852 do 1853 iend=len(snam(mm)),1,-1 1853 if(snam(mm)(iend:iend).ne.' ')goto 1854 1854 open(15, 1 FILE=snam(mm)(istart:iend)//'.int.ttt',status='old') write(6,7006)mm,' Raman',(snam(mm)(i:i),i=istart,iend) 7006 format(' Kind ',i4,', ',a6,' tensors from ',40a1) CALL READTTT(n,ALPHAI,AAI,GI,NAT1) if(NAT1.gt.2)call readatms(15,NAT1,iii,n) close(15) do 18541 L=1,NAT1 do 18542 I=1,4 18542 IIM(mm,I,L)=iii(I,L) do 18541 IX=1,3 IIND=3*(L-1)+IX do 18541 I=1,3 do 18541 J=1,3 ALPHAM(I+(J-1)*3,mm,IIND)=ALPHAI(IIND,I,J) GM(I+(J-1)*3,mm,IIND)=GI(IIND,I,J) do 18541 K=1,3 18541 AAM(I+(J-1)*3+(K-1)*9,mm,IIND)=AAI(IIND,I,J,K) open(15, 1 FILE=snam(mm)(istart:iend)//'.int.pol',status='old') CALL readpol(POLI) c read(15,*) close(15) do 186 IP=1,3 do 186 JP=1,3 186 POLM(mm,IP+(JP-1)*3)=POLI(IP,JP) else write(6,7007) 7007 format(' All molecules: INT.TTT tensors') open(15,FILE='INT.TTT',status='old') CALL READTTT(n,ALPHAI,AAI,GI,NAT1) if(NAT1.gt.2)call readatms(15,NAT1,iii,n) close(15) open(15,FILE='INT.POL',status='old') CALL readpol(POLI) c read(15,*) close(15) endif endif if(ltei)then if(lnnames)then do 189 mm=1,nnames do 1891 istart=1,len(snam(mm)) 1891 if(snam(mm)(istart:istart).ne.' ')goto 1892 1892 do 1893 iend=len(snam(mm)),1,-1 1893 if(snam(mm)(iend:iend).ne.' ')goto 1894 1894 open(15, 1 FILE=snam(mm)(istart:iend)//'.int.ten',status='old') write(6,7006)mm,'dipole',(snam(mm)(i:i),i=istart,iend) CALL READTENI(n,APTI,AI,NAT1) write(6,*)'NAT1',NAT1,' n:',n if(NAT1.gt.2)call readatms(15,NAT1,iii,n) CLOSE(15) do 189 ia=1,NAT1 do 18842 I=1,4 18842 IIM(mm,I,L)=iii(I,L) do 189 I=1,3 do 189 J=1,3 APTM(ia,I+(J-1)*3,mm)=APTI(ia,I,J) 189 AM( ia,I+(J-1)*3,mm)=AI(ia,I,J) else OPEN(15,FILE='INT.TEN',status='old') CALL READTENI(n,APTI,AI,NAT1) if(NAT1.gt.2)call readatms(15,NAT1,iii,n) CLOSE(15) write(6,7008) 7008 format(' All molecules: INT.TEN tensors') endif endif nw=int((wmax-wmin)/dw) if(nw.gt.nw0)then write(6,*)'too many spectral points' stop endif sc=1.0d0 dtau=dt*1.0d-12/2.418D-17 write(6,6013)iprint,dt,dtau,n,NSA,dw,wmin,wmax,nw,res,lper,perx, 1pery,perz,lnnames,lincr,lp27,lp28,lp29,lcell 6013 format(/,' printing interval : ',i10,' steps',/, 1' dt : ',f10.6,' ps (',f10.4,' au)',/, 1' n : ',i10,' atoms',/, 1' NSA : ',i10,' atoms',/, 1' DwA : ',f10.2,' cm-1',/, 1' wmin : ',f10.2,' cm-1',/, 1' wmax : ',f10.2,' cm-1',/, 1' Num. of freqs. : ',i10,/, 1' Convol. broadening: ',f10.2,' cm-1'/, 1' Per. box : ',l1,3f10.2,' A',/, 1' lnnames : ',l1,/, 1' lincr : ',l1,/, 1' lp27 lp28 lp29 lcell: ',4l3,/) c call dc(res,lg,dw,wmin,nw) dw =dw *tp*sol wmax=wmax*tp*sol wmin=wmin*tp*sol c c initiate frequency coefficients: allocate(wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0)) if(lincr)then call ifcr(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) else call ifc(wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) istep=0 endif do 24 ia=1,n xm(ia)=x(ia) ym(ia)=y(ia) zm(ia)=z(ia) xmm(ia)=x(ia) ymm(ia)=y(ia) 24 zmm(ia)=z(ia) inquire(file='CCT.INP',exist=lcct) if(lcct)then mxb=n allocate(R(3,n),IND(n),IB(n,mxb),ND(n),NDS(n),FF(3*n,3*n), 1 FFS(3*n,3*n),IS(n,mxb),INDBIG(IOX,n),IBONDSB(8,n), 2 IBONDS(8,n),NBO(n),NBOB(n),PS(n,3,3),AS(n,3,3),ALPHAS(3*n,3,3), 3 AAS(3*n,3,3,3),GS(3*n,3,3),AAT_c(3*n,3,3,3),GT_c(3*n,3,3), 1 RS(3,n),asname(IOX,n,3,3),psname(IOX,n,3,3), 1 aaname(IOX,3*n,3,3,3),gname(IOX,3*n,3,3),alphaname(IOX,3*n,3,3), 1 rname(IOX,3,n),NATNAME(IOX),ffname(IOX,3*n,3*n),NDNAME(IOX,n), 1 AM_c(n,3,3)) call cctn_tinker2(n,atomic,x,y,z,mxb, 1 R,IND,IB,ND,NDS,RS,FFS,IS,INDBIG,IBONDSB,IBONDS,NBO,NBOB, 2 PS,AS,ALPHAS,AAS,GS,AAT_c,GT_c, 3 asname,psname,aaname,gname,alphaname,rname,NATNAME,ffname, 4 NDNAME,IOX,luse) endif write(6,*)'animated lani ',lani write(6,*)'convolution lconv ',lconv write(6,*)'internal ten ltei ',ltei write(6,*)'internal ttt ltti ',ltti write(6,*)'at polar lexp ',lexp write(6,*)'at charges lexca ',lexca write(6,*)'cct ram LRAM ',LRAM write(6,*)'cct roa LROA ',LROA write(6,*)'cct abs LABS ',LABS write(6,*)'cct vcd LVCD ',LVCD write(6,*)'project pol lproject ',lproject write(6,*)'dress pol ldress ',ldress write(6,*)'conv type lg ',lg c ######################################################## open(10,file='FILE.X') 31111 istep=istep+1 c current time [s] t=dt*dble(istep)*1.0d-12 if(mod(istep,iprint).eq.0)then write(6,6003)istep,int(t*1.0d15) 6003 format(2i20,' fs') endif c read current geometry: read(10,*,end=31112,err=31112) read(10,*,end=31112,err=31112)n call rg(n,x,y,z,atomic,i12,0,maxval) c if(lcct)then call trtena(x,y,z,ALPHA1,G1,AA1,P,A,n,mxb, 1 R,IND,IB,RS,FF,FFS,IS,INDBIG,IBONDSB,NBOB, 2 PS,AS,ALPHAS,AAS,GS,IOX,AM_c,luse, 3 asname,psname,aaname,gname,alphaname,rname,NATNAME,ffname) c project translations and rotations from alpha if(lproject)call protr(x,y,z,ALPHA1,P,n,istep,luse) if(istep.eq.1)write(6,*)'Tensors by transfer' endif if(ltei)then if(istep.eq.1)write(6,*)'DD from internal tensors' do 154 mm=1,nmol c define the transformation matrices do 1561 ia=1,natm(mm) xx(ia)=x(al(mm,ia)) yy(ia)=y(al(mm,ia)) 1561 zz(ia)=z(al(mm,ia)) if(lnnames)then do 158 ia=1,natm(mm) do 158 I=1,4 158 iii(I,ia)=IIM(imk(mm),I,ia) endif call dtm(iii,natm(mm),u,n,xx,yy,zz) c do 154 ia=1,natm(mm) L=al(mm,ia) c c rewrite transformation to faster arrays: DO 14121 I=1,3 DO 14121 J=1,3 14121 ut(I,J)=u(I,J,ia) if(lnnames)then DO 141 I=1,3 DO 141 J=1,3 at1(I,J)= AM( ia,I+(J-1)*3,imk(mm)) 141 pt1(I,J)=APTM(ia,I+(J-1)*3,imk(mm)) else DO 1413 I=1,3 DO 1413 J=1,3 at1(I,J)= AI(ia,I,J) 1413 pt1(I,J)=APTI(ia,I,J) endif c DO 142 I=1,3 DO 142 JP=1,3 at2(I,JP)=0.0d0 pt2(I,JP)=0.0d0 do 142 IP=1,3 at2(I,JP)=at2(I,JP)+ut(IP,I)*at1(IP,JP) 142 pt2(I,JP)=pt2(I,JP)+ut(IP,I)*pt1(IP,JP) c DO 14 I=1,3 DO 14 J=1,3 at1(I,J)=0.0d0 pt1(I,J)=0.0d0 do 14 JP=1,3 at1(I,J)=at1(I,J)+ut(JP,J)*at2(I,JP) 14 pt1(I,J)=pt1(I,J)+ut(JP,J)*pt2(I,JP) c c Transform to common origin: FBOHR=0.25d0/bohr DO 2201 J=1,3 at1(J,1)=at1(J,1)+(y(ia)*pt1(J,3)-z(ia)*pt1(J,2))*FBOHR at1(J,2)=at1(J,2)+(z(ia)*pt1(J,1)-x(ia)*pt1(J,3))*FBOHR 2201 at1(J,3)=at1(J,3)+(x(ia)*pt1(J,2)-y(ia)*pt1(J,1))*FBOHR c c record DO 154 I=1,3 DO 154 J=1,3 A(L,I,J)=at1(I,J) 154 P(L,I,J)=pt1(I,J) endif c if(lexca)then if(istep.eq.1)write(6,*)'DD from charges' do 401 ia=1,n fac=pchgft(ia) do 405 ix=1,3 do 405 jx=1,3 P(ia,ix,jx)=0.0d0 405 A(ia,ix,jx)=0.0d0 P(ia,1,1)=fac*x(ia) P(ia,2,2)=fac*y(ia) P(ia,3,3)=fac*z(ia) A(ia,1,2)= fac*z(ia) A(ia,2,1)=-fac*z(ia) A(ia,2,3)= fac*x(ia) A(ia,3,2)=-fac*x(ia) A(ia,3,1)= fac*y(ia) 401 A(ia,1,3)=-fac*y(ia) endif c if(lexp)then if(istep.eq.1)write(6,*)'Alpha from atomic polarizabilities' c c molecular polarizabilities: DO 17 mm=1,npair DO 179 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 172 IX=1,3 IIND=3*(L-1)+IX DO 172 I=1,3 DO 172 J=1,3 AA1(IIND,I,J,1)=0.0d0 AA1(IIND,I,J,2)=0.0d0 AA1(IIND,I,J,3)=0.0d0 G1(IIND,I,J)=0.0d0 172 ALPHA1(IIND,I,J)=0.0d0 endif 179 continue DO 192 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 178 iap=1,natm(mm) LP=al(mm,iap) if(luse(LP))then rki(1)=x(LP)-x(L) rki(2)=y(LP)-y(L) rki(3)=z(LP)-z(L) if(LP.ne.L)then rd1=dsqrt(rki(1)**2+rki(2)**2+rki(3)**2+1.0d-9) rd67=6.0d0/rd1**7*pol(L)*pol(LP) DO 173 IX=1,3 IIND=3*(L-1)+IX DO 173 I=1,3 DO 173 J=1,3 s1=-5.0d0*rki(I)*rki(J)*rki(IX) if(I.eq.IX)s1=s1+rd1*rd1*rki(J) if(J.eq.IX)s1=s1+rd1*rd1*rki(I) if(I.eq.J )s1=s1+rd1*rd1*rki(IX) 173 ALPHA1(IIND,I,J)=ALPHA1(IIND,I,J)+s1*rd67 endif endif 178 continue endif 192 continue CALL TTTT(ALPHA1,AA1,G1,NAT,2,x,y,z,Xwork) c c transcript into the big arrays: DO 193 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 194 IX=1,3 IIND=3*(L-1)+IX DO 194 I=1,3 DO 194 J=1,3 DO 195 K=1,3 195 AA(I+(J-1)*3+(K-1)*9,mm,IIND)=AA1(IIND,I,J,K) G(I+(J-1)*3,mm,IIND)=G1(IIND,I,J) 194 ALPHA(I+(J-1)*3,mm,IIND)=ALPHA1(IIND,I,J) endif 193 continue 17 continue endif c molecular positions in atomic units (x y z are in A): DO 197 ii=1,3*nmol 197 rm(ii)=0.0d0 DO 196 mm=1,nmol DM=dble(natm(mm)) DO 198 ia=1,natm(mm) L=al(mm,ia) rm(1+(mm-1)*3)=rm(1+(mm-1)*3)+x(L) rm(2+(mm-1)*3)=rm(2+(mm-1)*3)+y(L) 198 rm(3+(mm-1)*3)=rm(3+(mm-1)*3)+z(L) rm(1+(mm-1)*3)=rm(1+(mm-1)*3)/DM/bohr rm(2+(mm-1)*3)=rm(2+(mm-1)*3)/DM/bohr 196 rm(3+(mm-1)*3)=rm(3+(mm-1)*3)/DM/bohr c if(ltti)then if(istep.eq.1)write(6,*)'Alpha from internal tensors' if(istep.eq.1)write(6,*)'Check the atomic order!' c propagate the internal tensors to all molecules, c order of atoms must be the same do 1515 mm=1,nmol c define the transformation matrices do 156 ia=1,natm(mm) xx(ia)=x(al(mm,ia)) yy(ia)=y(al(mm,ia)) 156 zz(ia)=z(al(mm,ia)) if(lnnames)then do 159 ia=1,natm(mm) do 159 I=1,4 159 iii(I,ia)=IIM(imk(mm),I,ia) endif c if not enough bonds: if(natm(mm).lt.2)then DO 1551 ia=1,natm(mm) DO 1551 I=1,3 DO 1552 J=1,3 1552 u(I,J,ia)=0.0d0 1551 u(I,I,ia)=1.0d0 else call dtm(iii,natm(mm),u,n,xx,yy,zz) endif do 15 ia=1,natm(mm) L=al(mm,ia) c c rewrite transformation to a faster array: DO 1412 I=1,3 DO 1412 J=1,3 if(u(I,J,ia).lt.-1.0d0.or.u(I,J,ia).gt.1.0d0)then write(6,*)ia,L,mm,I,J,u(I,J,ia) stop endif 1412 ut(I,J)=u(I,J,ia) c c polarizability according to the system of first atom: if(ia.eq.1)then DO 20 I=1,3 DO 20 J=1,3 20 po2(mm,I,J)=0.0d0 if(lnames)then DO 21 I=1,3 DO 21 J=1,3 DO 21 IP=1,3 DO 21 JP=1,3 21 po2(mm,I,J)=po2(mm,I,J) 1 +ut(IP,I)*POLM(imk(mm),IP+(JP-1)*3)*ut(JP,J) else DO 155 I=1,3 DO 155 J=1,3 DO 155 IP=1,3 DO 155 JP=1,3 155 po2(mm,I,J)=po2(mm,I,J)+ut(IP,I)*POLI(IP,JP)*ut(JP,J) endif endif c c rewrite ALPHA G tensors to faster arrays: if(lnnames)then DO 151 IXP=1,3 DO 151 IP=1,3 DO 151 JP=1,3 as1(IXP,IP,JP)=ALPHAM(IP+(JP-1)*3,imk(mm),3*(ia-1)+IXP) 151 gs1(IXP,IP,JP)=GM(IP+(JP-1)*3,imk(mm),3*(ia-1)+IXP) else DO 1512 IXP=1,3 DO 1512 IP=1,3 DO 1512 JP=1,3 as1(IXP,IP,JP)=ALPHAI(3*(ia-1)+IXP,IP,JP) 1512 gs1(IXP,IP,JP)=GI(3*(ia-1)+IXP,IP,JP) endif c c first index trafo DO 157 IX=1,3 DO 157 IP=1,3 DO 157 JP=1,3 as2(IX,IP,JP)=0.0d0 gs2(IX,IP,JP)=0.0d0 DO 157 IXP=1,3 as2(IX,IP,JP)=as2(IX,IP,JP)+ut(IXP,IX)*as1(IXP,IP,JP) 157 gs2(IX,IP,JP)=gs2(IX,IP,JP)+ut(IXP,IX)*gs1(IXP,IP,JP) c second index trafo DO 152 IX=1,3 DO 152 I=1,3 DO 152 JP=1,3 as1(IX,I,JP)=0.0d0 gs1(IX,I,JP)=0.0d0 DO 152 IP=1,3 as1(IX,I,JP)=as1(IX,I,JP)+ut(IP,I)*as2(IX,IP,JP) 152 gs1(IX,I,JP)=gs1(IX,I,JP)+ut(IP,I)*gs2(IX,IP,JP) c last index trafo DO 153 IX=1,3 IIND=3*(L-1)+IX DO 153 I=1,3 DO 153 J=1,3 ALPHA1(IIND,I,J)=0.0d0 G1(IIND,I,J)=0.0d0 DO 153 JP=1,3 ALPHA1(IIND,I,J)=ALPHA1(IIND,I,J)+ut(JP,J)*as1(IX,I,JP) 153 G1(IIND,I,J) =G1(IIND,I,J) +ut(JP,J)*gs1(IX,I,JP) c A-similarly as above alpha G: if(lnnames)then DO 16 IXP=1,3 IIND=3*(ia-1)+IXP DO 16 IP=1,3 DO 16 JP=1,3 DO 16 KP=1,3 16 au1(IXP,IP,JP,KP)=AAM(IP+(JP-1)*3+(KP-1)*3,imk(mm),IIND) else DO 165 IXP=1,3 IIND=3*(ia-1)+IXP DO 165 IP=1,3 DO 165 JP=1,3 DO 165 KP=1,3 165 au1(IXP,IP,JP,KP)=AAI(IIND,IP,JP,KP) endif c first DO 161 IX=1,3 DO 161 IP=1,3 DO 161 JP=1,3 DO 161 KP=1,3 au2(IX,IP,JP,KP)=0.0d0 DO 161 IXP=1,3 161 au2(IX,IP,JP,KP)=au2(IX,IP,JP,KP)+ut(IXP,IX)*au1(IXP,IP,JP,KP) c second DO 162 IX=1,3 DO 162 I =1,3 DO 162 JP=1,3 DO 162 KP=1,3 au1(IX,I,JP,KP)=0.0d0 DO 162 IP=1,3 162 au1(IX,I,JP,KP)=au1(IX,I,JP,KP)+ut(IP,I)*au2(IX,IP,JP,KP) c third DO 163 IX=1,3 DO 163 I =1,3 DO 163 J =1,3 DO 163 KP=1,3 au2(IX,I,J,KP)=0.0d0 DO 163 JP=1,3 163 au2(IX,I,J,KP)=au2(IX,I,J,KP)+ut(JP,J)*au1(IX,I,JP,KP) c fourth DO 15 IX=1,3 IIND=3*(L-1)+IX DO 15 I =1,3 DO 15 J =1,3 DO 15 K =1,3 AA1(IIND,I,J,K)=0.0d0 DO 15 KP=1,3 15 AA1(IIND,I,J,K)=AA1(IIND,I,J,K)+ut(KP,K)*au2(IX,I,J,KP) c transform to common origin: CALL TTTT(ALPHA1,AA1,G1,NAT,2,x,y,z,Xwork) do 1515 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 300 IX=1,3 IIND=3*(L-1)+IX DO 300 I=1,3 DO 300 J=1,3 DO 306 K=1,3 306 AA(I+(J-1)*3+(K-1)*9,mm,IIND)=AA1(IIND,I,J,K) G (I+(J-1)*3,mm,IIND)=G1 (IIND,I,J) 300 ALPHA(I+(J-1)*3,mm,IIND)=ALPHA1(IIND,I,J) endif 1515 continue c c dress the polarizabilities to account for the influence of the c other molecules if(ldress)then c calculate cube center qx=0.0d0 qy=0.0d0 qz=0.0d0 do 1003 mm=1,nmol qx=qx+rm(1+3*(mm-1)) qy=qy+rm(2+3*(mm-1)) 1003 qz=qz+rm(3+3*(mm-1)) qx=qx/dble(nmol) qy=qy/dble(nmol) qz=qz/dble(nmol) c **************************************** IF(lp29)THEN c true superbox c c molecules in the box1: do 11517 mm=1,nmol c save original polarizability for each molecule do 406 ia=1,natm(mm) L=al(mm,ia) DO 406 IX=1,3 IIND=3*(L-1)+IX DO 406 I=1,3 DO 406 J=1,3 406 ALP29(I+(J-1)*3,IIND)=ALPHA(I+(J-1)*3,mm,IIND) c temporary polarizability where all is collected do 410 ia=1,natm(mm) L=al(mm,ia) DO 410 IX=1,3 IIND=3*(L-1)+IX DO 410 I=1,3 DO 410 J=1,3 410 ALPHM(I+(J-1)*3,IIND)=0.0d0 c save static polarizability to faster array: do 409 I=1,3 do 409 J=1,3 409 p1(I+3*(J-1))=po2(mm,I,J) c box1 do 11516 ibx=1,3 dx=dble(ibx-2)*perx do 11516 iby=1,3 dy=dble(iby-2)*pery do 11516 ibz=1,3 dz=dble(ibz-2)*perz x1=rm(1+3*(mm-1))+dx y1=rm(2+3*(mm-1))+dy z1=rm(3+3*(mm-1))+dz c damp factor towards box walls: fac1=1.0d0 fac1=fac1*(1.0d0+cos(pi*(qx-x1)/px2/3.0d0))/2.0d0 fac1=fac1*(1.0d0+cos(pi*(qy-y1)/py2/3.0d0))/2.0d0 fac1=fac1*(1.0d0+cos(pi*(qz-z1)/pz2/3.0d0))/2.0d0 c damp original polarizability, add all to molecule mm: do 407 ia=1,natm(mm) L=al(mm,ia) DO 407 IX=1,3 IIND=3*(L-1)+IX DO 407 I=1,3 DO 407 J=1,3 407 ALPHM(I+(J-1)*3,IIND)=ALPHM(I+(J-1)*3,IIND) 1 +ALP29(I+(J-1)*3,IIND)*fac1 c molecules in the box2: do 11516 mp=1,nmol c save polarizability to faster array: do 408 I=1,3 do 408 J=1,3 408 p2(I+3*(J-1))=po2(mp,I,J) c box2 do 11516 ibxp=1,3 dxp=dble(ibxp-2)*perx do 11516 ibyp=1,3 dyp=dble(ibyp-2)*pery do 11516 ibzp=1,3 dzp=dble(ibzp-2)*perz x2=rm(1+3*(mp-1))+dxp y2=rm(2+3*(mp-1))+dyp z2=rm(3+3*(mp-1))+dzp rdx=x1-x2 rdy=y1-y2 rdz=z1-z2 rd2=rdx**2+rdy**2+rdz**2 rd5=rd2**2*dsqrt(rd2) c damp pair contributions towards box walls: fac=fac1 fac=fac*(1.0d0+cos(pi*(qx-x2)/px2/3.0d0))/2.0d0 fac=fac*(1.0d0+cos(pi*(qy-y2)/py2/3.0d0))/2.0d0 fac=fac*(1.0d0+cos(pi*(qz-z2)/pz2/3.0d0))/2.0d0 if(rd2.gt.1.0d-1)then c add all interactions to molecule mm: do 15162 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 402 IX=1,3 IIND=3*(L-1)+IX DO 402 I=1,3 DO 402 J=1,3 s1=p2(J )*ALP29(1+(I-1)*3,IIND) 1 + p2(J+3)*ALP29(2+(I-1)*3,IIND) 2 + p2(J+6)*ALP29(3+(I-1)*3,IIND) s2=p2(J)*rdx+p2(J+3)*rdy+p2(J+6)*rdz s3=ALP29(1+(I-1)*3,IIND)*rdx 1 + ALP29(2+(I-1)*3,IIND)*rdy 2 + ALP29(3+(I-1)*3,IIND)*rdz 402 ALPHM(I+(J-1)*3,IIND)=ALPHM(I+(J-1)*3,IIND)+ 1 2.0d0*fac*(3.0d0*s2*s3-s1*rd2)/rd5 c c corrections on molecular positions: DO 404 I=1,3 DO 404 J=1,3 s1=p2(I)*p1(J)+p2(I+3)*p1(J+3)+p2(I+6)*p1(J+6) s2=p2(I)*rdx+p2(I+3)*rdy+p2(I+6)*rdz s3=p1(J)*rdx+p1(J+3)*rdy+p1(J+6)*rdz smd=s3*p2(I) spd=s2*p1(J) ALPHM(I+(J-1)*3,3*(L-1)+1)=ALPHM(I+(J-1)*3,3*(L-1)+1) 1 +6.0d0*fac*(smd+spd+s1*rdx-5.0d0*s2*s3*rdx/rd2) 1 /DM/rd5 smd=s3*p2(I+3) spd=s2*p1(J+3) ALPHM(I+(J-1)*3,3*(L-1)+2)=ALPHM(I+(J-1)*3,3*(L-1)+2) 1 +6.0d0*fac*(smd+spd+s1*rdy-5.0d0*s2*s3*rdy/rd2) 1 /DM/rd5 smd=s3*p2(I+6) spd=s2*p1(J+6) 404 ALPHM(I+(J-1)*3,3*(L-1)+3)=ALPHM(I+(J-1)*3,3*(L-1)+3) 1 +6.0d0*fac*(smd+spd+s1*rdz-5.0d0*s2*s3*rdz/rd2) 1 /DM/rd5 endif 15162 continue endif 11516 continue c record resultant derivatives: do 11517 ia=1,natm(mm) L=al(mm,ia) DO 11517 IX=1,3 IIND=3*(L-1)+IX DO 11517 I=1,3 DO 11517 J=1,3 11517 ALPHA(I+(J-1)*3,mm,IIND)=ALPHM(I+(J-1)*3,IIND) c **************************************** ELSE c c set superbox if required if(lp27)then nib=3 else nib=1 endif c c molecules in the box: do 1517 mm=1,nmol DM=dble(natm(mm)) c save original polarizability for each molecule do 411 ia=1,natm(mm) L=al(mm,ia) DO 411 IX=1,3 IIND=3*(L-1)+IX DO 411 I=1,3 DO 411 J=1,3 411 ALP29(I+(J-1)*3,IIND)=ALPHA(I+(J-1)*3,mm,IIND) c temporary polarizability where all is collected do 412 ia=1,natm(mm) L=al(mm,ia) DO 412 IX=1,3 IIND=3*(L-1)+IX DO 412 I=1,3 DO 412 J=1,3 412 ALPHM(I+(J-1)*3,IIND)=0.0d0 c save static polarizability to faster array: do 413 I=1,3 do 413 J=1,3 413 p1(I+3*(J-1))=po2(mm,I,J) x1=rm(1+3*(mm-1)) y1=rm(2+3*(mm-1)) z1=rm(3+3*(mm-1)) c molecules in the box and the boxes around (if lp27): do 1516 ibx=1,nib dx=dble(ibx-nib+1)*perx do 1516 iby=1,nib dy=dble(iby-nib+1)*pery do 1516 ibz=1,nib dz=dble(ibz-nib+1)*perz do 1516 mp=1,nmol if(mm.ne.mp.or. 1 (lp27.and.(ibx.ne.2.or.iby.ne.2.or.ibz.ne.2)))then c save polarizability to faster array: do 414 I=1,3 do 414 J=1,3 414 p2(I+3*(J-1))=po2(mp,I,J) x2=rm(1+3*(mp-1)) y2=rm(2+3*(mp-1)) z2=rm(3+3*(mp-1)) if(nib.eq.3)then x2=x2+dx y2=y2+dy z2=z2+dz endif rdx=x1-x2 rdy=y1-y2 rdz=z1-z2 rd2=rdx**2+rdy**2+rdz**2 rd5=rd2**2*dsqrt(rd2) fac=1.0d0 if(lp27)then c if(ibx.eq.1)fac=fac*(1.0d0+cos(pi*(qx-px2-x2)/perx))/2.0d0 c if(ibx.eq.3)fac=fac*(1.0d0+cos(pi*(qx+px2-x2)/perx))/2.0d0 c if(iby.eq.1)fac=fac*(1.0d0+cos(pi*(qy-py2-y2)/pery))/2.0d0 c if(iby.eq.3)fac=fac*(1.0d0+cos(pi*(qy+py2-y2)/pery))/2.0d0 c if(ibz.eq.1)fac=fac*(1.0d0+cos(pi*(qz-pz2-z2)/perz))/2.0d0 c if(ibz.eq.3)fac=fac*(1.0d0+cos(pi*(qz+pz2-z2)/perz))/2.0d0 endif c damp towards box walls: if(lp28)then fac=fac*(1.0d0+cos(pi*(qx-x2)/px2))/2.0d0 fac=fac*(1.0d0+cos(pi*(qy-y2)/py2))/2.0d0 fac=fac*(1.0d0+cos(pi*(qz-z2)/pz2))/2.0d0 endif if(rd5.eq.0.0d0)then write(6,*)mm,mp call report(' zero distance between molecules!') endif do 15161 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 400 IX=1,3 IIND=3*(L-1)+IX DO 400 I=1,3 DO 400 J=1,3 s1=p2(J )*ALP29(1+(I-1)*3,IIND) 1 + p2(J+3)*ALP29(2+(I-1)*3,IIND) 2 + p2(J+6)*ALP29(3+(I-1)*3,IIND) s2=p2(J)*rdx+p2(J+3)*rdy+p2(J+6)*rdz s3=ALP29(1+(I-1)*3,IIND)*rdx 1 + ALP29(2+(I-1)*3,IIND)*rdy 2 + ALP29(3+(I-1)*3,IIND)*rdz 400 ALPHM(I+(J-1)*3,IIND)=ALPHM(I+(J-1)*3,IIND)+ 1 2.0d0*fac*(3.0d0*s2*s3-s1*rd2)/rd5 c c corrections on molecular positions: DO 403 I=1,3 DO 403 J=1,3 s1=p2(I)*p1(J)+p2(I+3)*p1(J+3)+p2(I+6)*p1(J+6) s2=p2(I)*rdx+p2(I+3)*rdy+p2(I+6)*rdz s3=p1(J)*rdx+p1(J+3)*rdy+p1(J+6)*rdz ALPHM(I+(J-1)*3,3*(L-1)+1)=ALPHM(I+(J-1)*3,3*(L-1)+1) 1 +6.0d0*fac*(s3*p2(I )+s2*p1(J )+(s1-5.0d0*s2*s3/rd2)*rdx) 1 /DM/rd5 ALPHM(I+(J-1)*3,3*(L-1)+2)=ALPHM(I+(J-1)*3,3*(L-1)+2) 1 +6.0d0*fac*(s3*p2(I+3)+s2*p1(J+3)+(s1-5.0d0*s2*s3/rd2)*rdy) 1 /DM/rd5 403 ALPHM(I+(J-1)*3,3*(L-1)+3)=ALPHM(I+(J-1)*3,3*(L-1)+3) 1 +6.0d0*fac*(s3*p2(I+6)+s2*p1(J+6)+(s1-5.0d0*s2*s3/rd2)*rdz) 1 /DM/rd5 endif 15161 continue endif 1516 continue c colect back to slow array: do 1517 ia=1,natm(mm) L=al(mm,ia) DO 1517 IX=1,3 IIND=3*(L-1)+IX DO 1517 I=1,3 DO 1517 J=1,3 1517 ALPHA(I+(J-1)*3,mm,IIND)=ALP29(I+(J-1)*3,IIND)+ 1 ALPHM(I+(J-1)*3,IIND) ENDIF c **************************************** if(istep.eq.1)write(6,*)'Alpha dressed' endif endif if(lcct)then if(istep.eq.1)write(6,*)'Alpha from transferred ' DO 201 mm=1,nmol DO 199 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 200 IX=1,3 IIND=3*(L-1)+IX DO 200 I=1,3 DO 200 J=1,3 DO 206 K=1,3 206 AA(I+(J-1)*3+(K-1)*9,mm,IIND)=AA1(IIND,I,J,K) G (I+(J-1)*3,mm,IIND)=G1 (IIND,I,J) 200 ALPHA(I+(J-1)*3,mm,IIND)=ALPHA1(IIND,I,J) endif 199 continue 201 continue endif c c correct for periodic jumps: if(lper)then if(lcell)then aan=dsqrt(ac(1)**2+ac(2)**2+ac(3)**2) abn=dsqrt(bc(1)**2+bc(2)**2+bc(3)**2) acn=dsqrt(cc(1)**2+cc(2)**2+cc(3)**2) do 185 ia=1,n ap=(x(ia)*ac(1)+y(ia)*ac(2)+z(ia)*ac(3))/aan bp=(x(ia)*bc(1)+y(ia)*bc(2)+z(ia)*bc(3))/abn cp=(x(ia)*cc(1)+y(ia)*cc(2)+z(ia)*cc(3))/acn apm=(xm(ia)*ac(1)+ym(ia)*ac(2)+zm(ia)*ac(3))/aan bpm=(xm(ia)*bc(1)+ym(ia)*bc(2)+zm(ia)*bc(3))/abn cpm=(xm(ia)*cc(1)+ym(ia)*cc(2)+zm(ia)*cc(3))/acn if(ap-apm.gt.aan)then call trabc(ia,xm ,ym ,zm ,ac,1.0d0) call trabc(ia,xmm,ymm,zmm,ac,1.0d0) endif if(apm-ap.gt.aan)then call trabc(ia,xm ,ym ,zm ,ac,-1.0d0) call trabc(ia,xmm,ymm,zmm,ac,-1.0d0) endif if(bp-bpm.gt.abn)then call trabc(ia,xm , ym, zm,bc,1.0d0) call trabc(ia,xmm,ymm,zmm,bc,1.0d0) endif if(bpm-bp.gt.abn)then call trabc(ia, xm, ym, zm,bc,-1.0d0) call trabc(ia,xmm,ymm,zmm,bc,-1.0d0) endif if(cp-cpm.gt.acn)then call trabc(ia, xm, ym, zm,cc,1.0d0) call trabc(ia,xmm,ymm,zmm,cc,1.0d0) endif if(cpm-cp.gt.acn)then call trabc(ia, xm, ym, zm,cc,-1.0d0) call trabc(ia,xmm,ymm,zmm,cc,-1.0d0) endif 185 continue else do 180 ia=1,n if(x(ia)-xm(ia).gt.px2)then xm(ia)=xm(ia)+perx xmm(ia)=xmm(ia)+perx endif if(xm(ia)-x(ia).gt.px2)then xm(ia)=xm(ia)-perx xmm(ia)=xmm(ia)-perx endif if(y(ia)-ym(ia).gt.py2)then ym(ia)=ym(ia)+pery ymm(ia)=ymm(ia)+pery endif if(ym(ia)-y(ia).gt.py2)then ym(ia)=ym(ia)-pery ymm(ia)=ymm(ia)-pery endif if(z(ia)-zm(ia).gt.pz2)then zm(ia)=zm(ia)+perz zmm(ia)=zmm(ia)+perz endif if(zm(ia)-z(ia).gt.pz2)then zm(ia)=zm(ia)-perz zmm(ia)=zmm(ia)-perz endif 180 continue endif endif c c update dipoles: if(ltei.or.lexca.or.LABS.or.LVCD)then c dipole time derivatives du/dt=du/dx . du/dt: c constants changing vx,ax into atomic units fac=1.0d0/bohr/dtau cf2=1.0d0/bohr/dtau/dtau do 209 mm=1,nmol do 2010 ia=1,6 2010 uu((mm-1)*6+ia)=0.0d0 do 209 ia=1,natm(mm) L=al(mm,ia) vx=(x(L)-xm(L))*fac vy=(y(L)-ym(L))*fac vz=(z(L)-zm(L))*fac ax=(x(L)+xmm(L)-xm(L)*2.0d0)*cf2 ay=(y(L)+ymm(L)-ym(L)*2.0d0)*cf2 az=(z(L)+zmm(L)-zm(L)*2.0d0)*cf2 c uu-index: 1 - electric moment x c 2 - electric moment y c 3 - electric moment z c 4 - magnetic moment x c 5 - magnetic moment y c 6 - magnetic moment z uu((mm-1)*6+1)=uu((mm-1)*6+1)+vx*P(L,1,1)+vy*P(L,2,1)+vz*P(L,3,1) uu((mm-1)*6+2)=uu((mm-1)*6+2)+vx*P(L,1,2)+vy*P(L,2,2)+vz*P(L,3,2) uu((mm-1)*6+3)=uu((mm-1)*6+3)+vx*P(L,1,3)+vy*P(L,2,3)+vz*P(L,3,3) uu((mm-1)*6+4)=uu((mm-1)*6+4)+ax*A(L,1,1)+ay*A(L,2,1)+az*A(L,3,1) uu((mm-1)*6+5)=uu((mm-1)*6+5)+ax*A(L,1,2)+ay*A(L,2,2)+az*A(L,3,2) 209 uu((mm-1)*6+6)=uu((mm-1)*6+6)+ax*A(L,1,3)+ay*A(L,2,3)+az*A(L,3,3) if(lconv)then call update2c(wmin,dw,nw,wsu,uu,t,sc,nmol) else call update2(wmin,dw,nw,wsu,uu,t,sc,nw0,nmol) endif endif c c update polarizabilities: if(ltti.or.lexp.or.LROA.or.LRAM)then c polarizability time derivatives da/dt=da/dx . dx/dt: fac=1.0d0/dtau/bohr do 2002 mm=1,nmol do 2001 ix=1,9 at(ix,mm)=0.0d0 gt(ix,mm)=0.0d0 aat(ix,mm)=0.0d0 aat(ix+9,mm)=0.0d0 2001 aat(ix+18,mm)=0.0d0 do 2002 ia=1,natm(mm) L=al(mm,ia) vx=(x(L)-xm(L))*fac vy=(y(L)-ym(L))*fac vz=(z(L)-zm(L))*fac do 2002 ii=1,9 gt(ii,mm)=gt(ii,mm) +G(ii,mm,3*(L-1)+1)*vx 1 + G(ii,mm,3*(L-1)+2)*vy 2 + G(ii,mm,3*(L-1)+3)*vz at(ii,mm)=at(ii,mm) +ALPHA(ii,mm,3*(L-1)+1)*vx 1 + ALPHA(ii,mm,3*(L-1)+2)*vy 2 + ALPHA(ii,mm,3*(L-1)+3)*vz aat(ii ,mm)=aat(ii ,mm)+AA(ii ,mm,3*(L-1)+1)*vx 1 + AA(ii ,mm,3*(L-1)+2)*vy 2 + AA(ii ,mm,3*(L-1)+3)*vz aat(ii+ 9,mm)=aat(ii+ 9,mm)+AA(ii+ 9,mm,3*(L-1)+1)*vx 1 + AA(ii+ 9,mm,3*(L-1)+2)*vy 2 + AA(ii+ 9,mm,3*(L-1)+3)*vz 2002 aat(ii+18,mm)=aat(ii+18,mm)+AA(ii+18,mm,3*(L-1)+1)*vx 1 + AA(ii+18,mm,3*(L-1)+2)*vy 2 + AA(ii+18,mm,3*(L-1)+3)*vz c if(lconv)then call update3c(wmin,dw,nw,wsat, 1 wcat,wsgt,wcgt,wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol) else call update3(wmin,dw,nw,wsat, 1 wcat,wsgt,wcgt,wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol,nw0) endif endif c c remember last point do 1 ia=1,n xmm(ia)=xm(ia) ymm(ia)=ym(ia) zmm(ia)=zm(ia) xm(ia)=x(ia) ym(ia)=y(ia) 1 zm(ia)=z(ia) c if(mod(istep,iprint).eq.0)then c write the normalized power spectra if(ltei.or.lexca.or.LABS.or.LVCD)then c abs/vcd: do 711 iw=1,nw rams(1,iw)=0.0d0 rams(2,iw)=0.0d0 do 711 mm=1,nmol w=wmin+dble(iw-1)*dw fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol 1 /dble(istep) ii=(iw-1)*12+(mm-1)*12*nw0 rams(1,iw)=rams(1,iw)+fac* 1 (wsu(1+ii)**2+wsu(3+ii)**2+wsu(5+ii)**2 2 +wsu(2+ii)**2+wsu(4+ii)**2+wsu(6+ii)**2) 711 rams(2,iw)=rams(2,iw)+fac* 1 (wsu(1+ii)*wsu(8+ii)+wsu(3+ii)*wsu(10+ii)+wsu(5+ii)*wsu(12+ii) 2 -wsu(2+ii)*wsu(7+ii)+wsu(9+ii)*wsu(10+ii)+wsu(6+ii)*wsu(11+ii)) c do 1238 iw=1,nw 1238 spec(iw)=rams(1,iw) if(lconv)then call wrspec(istep,lani,'WD.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WDraw.PRN',nw,wmin,dw,spec,tp,sol) c convolute it with a smoothing function: call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WD.PRN',nw,wmin,dw,spec,tp,sol) endif c do 1239 iw=1,nw 1239 spec(iw)=rams(2,iw) if(lconv)then call wrspec(istep,lani,'WR.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WRraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WR.PRN',nw,wmin,dw,spec,tp,sol) endif endif c Raman: if(ltti.or.lexp.or.LROA.or.LRAM)then do 1233 iw=1,nw w=wmin+dble(iw-1)*dw rams(1,iw)=0.0d0 rams(2,iw)=0.0d0 rams(3,iw)=0.0d0 rams(4,iw)=0.0d0 do 1233 mm=1,nmol c temperature factor:*1d5(arbitrary factor) if(dabs(w).lt.1.0d-40)then tf=0.0d0 else tf=1.0d0/w/(1.0d0-exp(-w*1.44d0/kelvin))*1.0d5 endif if(lexp)then tf=tf*1.0d-10 else tf=tf*1.0d6 endif sal0=0.0d0 sal1=0.0d0 sag0=0.0d0 sag1=0.0d0 sa1=0.0d0 do 121 jx=1,3 ii=iw+(mm-1)*nw0 sa1=sa1+wsat(1+(jx-1)*3,ii)* 1 (wsaat(2+(3-1)*3+(jx-1)*9,ii)-wsaat(3+(2-1)*3+(jx-1)*9,ii)) 2 +wcat(1+(jx-1)*3,ii)* 2 (wcaat(2+(3-1)*3+(jx-1)*9,ii)-wcaat(3+(2-1)*3+(jx-1)*9,ii)) 3 +wsat(2+(jx-1)*3,ii)* 3 (wsaat(3+(1-1)*3+(jx-1)*9,ii)-wsaat(1+(3-1)*3+(jx-1)*9,ii)) 4 +wcat(2+(jx-1)*3,ii)* 4 (wcaat(3+(1-1)*3+(jx-1)*9,ii)-wcaat(1+(3-1)*3+(jx-1)*9,ii)) 5 +wsat(3+(jx-1)*3,ii)* 5 (wsaat(1+(2-1)*3+(jx-1)*9,ii)-wsaat(2+(1-1)*3+(jx-1)*9,ii)) 6 +wcat(3+(jx-1)*3,ii)* 6 (wcaat(1+(2-1)*3+(jx-1)*9,ii)-wcaat(2+(1-1)*3+(jx-1)*9,ii)) do 121 ix=1,3 sal0=sal0+wsat(ix+(ix-1)*3,ii)*wsat(jx+(jx-1)*3,ii) 1 +wcat(ix+(ix-1)*3,ii)*wcat(jx+(jx-1)*3,ii) sal1=sal1+wsat(ix+(jx-1)*3,ii)*wsat(ix+(jx-1)*3,ii) 1 +wcat(ix+(jx-1)*3,ii)*wcat(ix+(jx-1)*3,ii) c coordinates sag0=sag0+wsat(ix+(ix-1)*3,ii)*wsgt(jx+(jx-1)*3,ii) 1 +wcat(ix+(ix-1)*3,ii)*wcgt(jx+(jx-1)*3,ii) sag1=sag1+wsat(ix+(jx-1)*3,ii)*wsgt(ix+(jx-1)*3,ii) 1 +wcat(ix+(jx-1)*3,ii)*wcgt(ix+(jx-1)*3,ii) 121 continue sa1=sa1/3.0d0 c IL+IR Backscattering: s180=7.0d0*sal1+sal0 c IL+IR Backscattering-xpol: ydy=3.0d0*sal1-sal0 c IL+IR Backscattering-ypol: ydx=4.0d0*sal1+2.0d0*sal0 c IR-IL Backscattering: d180=8.0d0*(3.0d0*sag1-sag0+sa1) ramconst=(w/sol/1000.0d0)**2/dble(istep) rams(1,iw)=rams(1,iw)+s180*tf*ramconst rams(2,iw)=rams(2,iw)+ydy *tf*ramconst rams(3,iw)=rams(3,iw)+ydx *tf*ramconst 1233 rams(4,iw)=rams(4,iw)+d180*tf*ramconst do 1234 iw=1,nw 1234 spec(iw)=rams(1,iw) if(lconv)then call wrspec(istep,lani,'WRam180.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WRam180raw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WRam180.PRN',nw,wmin,dw,spec,tp,sol) endif c do 1235 iw=1,nw c1235 spec(iw)=rams(2,iw) c call wrspec(istep,lani,'WRamyraw.PRN',nw,wmin,dw,spec,tp,sol) c call conv(res,nw,stem,spec,wmin,dw,tp,sol) c call wrspec(istep,lani,'WRamy.PRN',nw,wmin,dw,spec,tp,sol) c do 1236 iw=1,nw c1236 spec(iw)=rams(3,iw) c call wrspec(istep,lani,'WRamxraw.PRN',nw,wmin,dw,spec,tp,sol) c call conv(res,nw,stem,spec,wmin,dw,tp,sol) c call wrspec(istep,lani,'WRamx.PRN',nw,wmin,dw,spec,tp,sol) c do 1237 iw=1,nw c1237 spec(iw)=rams(4,iw) do 1237 iw=1,nw 1237 spec(iw)=rams(4,iw) if(lconv)then call wrspec(istep,lani,'WRoa180.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WRoa180raw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WRoa180.PRN',nw,wmin,dw,spec,tp,sol) endif endif endif c ######################################################## if(lincr.and.mod(istep,iprint).eq.0)call ifcw(istep, 1wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) goto 31111 31112 close(10) write(6,6002)istep 6002 format(i10,' geometries in FILE.X') stop end subroutine conv(res,nw,stem,spec,wmin,dw,tp,sol) implicit none real*8 res,stem(*),spec(*),wmin,dw,fac,tp,sol,w,wp integer*4 nw,iw,iwp if(res.gt.0)then do 72 iw=1,nw 72 stem(iw)=spec(iw) do 73 iw=1,nw w=(wmin+dble(iw-1)*dw)/tp/sol spec(iw)=0.0d0 do 73 iwp=1,nw wp=(wmin+dble(iwp-1)*dw)/tp/sol fac=((wp-w)/res)**2 73 if(fac.lt.20.0d0)spec(iw)=spec(iw)+stem(iwp)*exp(-fac) endif return end subroutine wrspec(istep,lani,sn,nw,wmin,dw,spec,tp,sol) implicit none character*(*) sn real*8 wmin,dw,spec(*),w,tp,sol integer*4 nw,iw,is,istep logical*4 lani character*20 ts if(lani)then write(ts,90)istep 90 format(i20) do 2 is=1,20 2 if(ts(is:is).eq.' ')ts(is:is)='0' open(8,file=ts(1:20)//sn) else open(8,file=sn) endif do 1 iw=1,nw w=wmin+dble(iw-1)*dw 1 write(8,805)w/tp/sol,spec(iw) 805 format(f12.2,' ',g15.6) close(8) return end subroutine update2(wmin,dw,nw,wsu,u,t,sc,nw0,nmol) implicit none integer*4 nw,iw,nmol,mm,ii,nw0,ix real*8 wmin,w,t,ss,cc,sc,dw,tp,sol,wsu(*),wau,u(*) c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 1 mm=1,nmol ii=(iw-1)*12+(mm-1)*12*nw0 do 1 ix=1,6 wsu( 2*ix-1+ii)=wsu(2*ix-1+ii)*sc+u((mm-1)*6+ix)*ss 1 wsu( 2*ix +ii)=wsu(2*ix +ii)*sc+u((mm-1)*6+ix)*cc return end subroutine update2c(wmin,dw,nw,wsu,u,t,sc,nmol) c use convolution implicit none integer*4 nw,iw,nmol,mm,ii,nw0,ix,nc,iwp parameter (nw0=2000) real*8 wmin,w,t,ss,cc,sc,dw,tp,sol,wsu(*),wau,u(*),fc,ff common/convc/ff(2*nw0+1),nc c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 1 iwp=max(1,iw-nc+1),min(iw+nc-1,nw) fc=ff(nc+iwp-iw) do 1 mm=1,nmol ii=(iwp-1)*12+(mm-1)*12*nw0 do 1 ix=1,6 wsu(2*ix-1+ii)=wsu(2*ix-1+ii)*sc+u((mm-1)*6+ix)*ss*fc 1 wsu(2*ix +ii)=wsu(2*ix +ii)*sc+u((mm-1)*6+ix)*cc*fc return end subroutine update3c(wmin,dw,nw,wsat,wcat,wsgt,wcgt, 1wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol) c alpha,A,G -> Raman, ROA c use convolution implicit none integer*4 nw,iw,ix,nn0,nmol,mm,ii,nw0 parameter (nw0=2000) real*8 wmin,w,at(9,nmol),gt(9,nmol),aat(27,nmol), 1t,ss,cc,sc,dw,tp,sol, 2wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0),wcgt(9,nn0), 2wsaat(27,nn0),wcaat(27,nn0),wau integer*4 nc,iwp real*8 fc,ff common/convc/ff(2*nw0+1),nc c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw c do iwp=1,2*nc-1 c write(6,*)iwp,ff(iwp) c enddo c stop do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 1 iwp=max(1,iw-nc+1),min(iw+nc-1,nw) fc=ff(nc+iwp-iw) do 1 mm=1,nmol ii=iwp+(mm-1)*nw0 do 1 ix=1,9 wsat(ix,ii)=wsat(ix,ii)*sc+at(ix,mm)*ss*fc wcat(ix,ii)=wcat(ix,ii)*sc+at(ix,mm)*cc*fc wsgt(ix,ii)=wsgt(ix,ii)*sc+gt(ix,mm)*ss*fc wcgt(ix,ii)=wcgt(ix,ii)*sc+gt(ix,mm)*cc*fc wsaat(ix,ii)=wsaat(ix,ii)*sc+aat(ix,mm)*ss*fc wcaat(ix,ii)=wcaat(ix,ii)*sc+aat(ix,mm)*cc*fc wsaat(ix+9,ii)=wsaat(ix+9,ii)*sc+aat(ix+9,mm)*ss*fc wcaat(ix+9,ii)=wcaat(ix+9,ii)*sc+aat(ix+9,mm)*cc*fc wsaat(ix+18,ii)=wsaat(ix+18,ii)*sc+aat(ix+18,mm)*ss*fc 1 wcaat(ix+18,ii)=wcaat(ix+18,ii)*sc+aat(ix+18,mm)*cc*fc return end subroutine update3(wmin,dw,nw,wsat,wcat,wsgt,wcgt, 1wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol,nw0) c alpha,A,G -> Raman, ROA implicit none integer*4 nw,iw,ix,nn0,nmol,mm,ii,nw0 real*8 wmin,w,at(9,nmol),gt(9,nmol),aat(27,nmol), 1t,ss,cc,sc,dw,tp,sol, 2wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0),wcgt(9,nn0), 2wsaat(27,nn0),wcaat(27,nn0),wau c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 1 mm=1,nmol ii=iw+(mm-1)*nw0 do 1 ix=1,9 wsat(ix,ii)=wsat(ix,ii)*sc+at(ix,mm)*ss wcat(ix,ii)=wcat(ix,ii)*sc+at(ix,mm)*cc wsgt(ix,ii)=wsgt(ix,ii)*sc+gt(ix,mm)*ss wcgt(ix,ii)=wcgt(ix,ii)*sc+gt(ix,mm)*cc wsaat(ix,ii)=wsaat(ix,ii)*sc+aat(ix,mm)*ss wcaat(ix,ii)=wcaat(ix,ii)*sc+aat(ix,mm)*cc wsaat(ix+9,ii)=wsaat(ix+9,ii)*sc+aat(ix+9,mm)*ss wcaat(ix+9,ii)=wcaat(ix+9,ii)*sc+aat(ix+9,mm)*cc wsaat(ix+18,ii)=wsaat(ix+18,ii)*sc+aat(ix+18,mm)*ss 1 wcaat(ix+18,ii)=wcaat(ix+18,ii)*sc+aat(ix+18,mm)*cc return end subroutine ifc(wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) implicit none integer*4 ix,nn0,ii real*8 wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0) do 2 ii=1,12*nn0 2 wsu(ii)=0.0d0 do 1 ii=1,nn0 do 1 ix=1,9 wsat(ix,ii)=0.0d0 wcat(ix,ii)=0.0d0 wsgt(ix,ii)=0.0d0 wcgt(ix,ii)=0.0d0 wcaat(ix,ii)=0.0d0 wsaat(ix,ii)=0.0d0 wcaat(ix+9,ii)=0.0d0 wsaat(ix+9,ii)=0.0d0 wcaat(ix+18,ii)=0.0d0 1 wsaat(ix+18,ii)=0.0d0 write(6,*)'coefficients initiated' return end subroutine ifcr(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) implicit none integer*4 ix,nn0,ii,istep real*8 wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0) logical lex inquire(file='RESTART.SCR',exist=lex) if(.not.lex)then istep=0 call ifc(wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) write(6,*)'RESTART.SCR not found' else open(9,file='RESTART.SCR',form='unformatted') read(9)istep read(9)(wsu(ii),ii=1,12*nn0) read(9)((wsat(ix,ii),ix=1,9),ii=1,nn0) read(9)((wcat(ix,ii),ix=1,9),ii=1,nn0) read(9)((wsgt(ix,ii),ix=1,9),ii=1,nn0) read(9)((wcgt(ix,ii),ix=1,9),ii=1,nn0) read(9)((wsaat(ix,ii),ix=1,27),ii=1,nn0) read(9)((wcaat(ix,ii),ix=1,27),ii=1,nn0) close(9) write(6,*)'RESTART.SCR read' endif return end subroutine ifcw(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) implicit none integer*4 ix,nn0,ii,istep real*8 wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0) open(9,file='RESTART.SCR',form='unformatted') write(9)istep write(9)(wsu(ii),ii=1,12*nn0) write(9)((wsat(ix,ii),ix=1,9),ii=1,nn0) write(9)((wcat(ix,ii),ix=1,9),ii=1,nn0) write(9)((wsgt(ix,ii),ix=1,9),ii=1,nn0) write(9)((wcgt(ix,ii),ix=1,9),ii=1,nn0) write(9)((wsaat(ix,ii),ix=1,27),ii=1,nn0) write(9)((wcaat(ix,ii),ix=1,27),ii=1,nn0) close(9) write(6,*)'RESTART.SCR updated' return end SUBROUTINE READTTT(N,ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*N,3,3),G(3*N,3,3),A(3*N,3,3,3) READ(15,*) READ(15,*)NAT READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) RETURN END SUBROUTINE READTENI(N,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N,3,3),A(N,3,3) READ (15,*)NAT DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) DO 230 L=1,NAT DO 230 J=1,3 230 READ (15,*) DO 100 L=1,NAT DO 100 J=1,3 100 READ (15,*) RETURN END subroutine readatms(io,NAT,iii,n) integer*4 n,NAT,iii(4,n),ia,io do 1 ia=1,NAT 1 read(io,*)(iii(ii,ia),ii=1,1),(iii(ii,ia),ii=1,4) return end subroutine dtm(ii,NAT,u,n,x,y,z) implicit none integer*4 n,NAT,ii(4,n),ia,i1,i2,j1,j2,ix,jx,ixx real*8 u(3,3,n),x(*),v1(3),v2(3),v3(3),y(*),z(*),dxx c if(NAT.eq.1)then do 2 ix=1,3 do 2 jx=1,3 2 u(jx,ix,1)=0.0d0 do 3 jx=1,3 3 u(jx,jx,1)=1.0d0 else if(NAT.eq.2)then do 4 ia=1,NAT i1=ii(1,ia) i2=ii(2,ia) v1(1)=x(i2)-x(i1) v1(2)=y(i2)-y(i1) v1(3)=z(i2)-z(i1) v2(1)=0.0d0 v2(2)=0.0d0 v2(3)=0.0d0 dxx=dabs(v1(1)) ixx=1 if(dabs(v1(2)).lt.dxx)then ixx=2 dxx=dabs(v1(2)) endif if(dabs(v1(3)).lt.dxx)ixx=3 v2(ixx)=1.0d0 call make3v(v1,v2,v3) do 4 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 4 u(3,ix,ia)=v3(ix) else do 1 ia=1,NAT i1=ii(1,ia) i2=ii(2,ia) j1=ii(3,ia) j2=ii(4,ia) if(i1.gt.0.and.i2.gt.0.and.j1.gt.0.and.j2.gt.0)then v1(1)=x(i2)-x(i1) v2(1)=x(j2)-x(j1) v1(2)=y(i2)-y(i1) v2(2)=y(j2)-y(j1) v1(3)=z(i2)-z(i1) v2(3)=z(j2)-z(j1) call make3v(v1,v2,v3) do 13 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 13 u(3,ix,ia)=v3(ix) else do 14 ix=1,3 do 14 jx=1,3 14 u(jx,ix,ia)=0.0d0 u(1,1,ia)=1.0d0 u(2,2,ia)=1.0d0 u(3,3,ia)=1.0d0 endif 1 continue endif endif return end SUBROUTINE TTTT(ALPHA,A,G,NAT,ICO,xx,y,z,X) C Transforms A and G to local origins (ICO=1) c or 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) DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3),xx(*), 1y(*),z(*),X(3,NAT) C C transfer into atomic coordinates: DO 5 J=1,NAT X(1,J)=xx(J)/0.529177D0 X(2,J)= y(J)/0.529177D0 5 X(3,J)= z(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(6,*)' G transformed into atomic origins' c IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' c 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 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(6,*)' A transformed into atomic origins' c IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' c IF(ICO.EQ.3)WRITE(6,*)' A set to zero' RETURN END subroutine report(s) character*(*) s write(6,*)s stop end FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF(I.EQ.1)then if(J.EQ.2.AND.K.EQ.3)then EPS= 1.0d0 else IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 endif else IF(I.EQ.2)then IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 else IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 endif endif RETURN END subroutine norm(v) implicit none real*8 v(*),n n=dsqrt(v(1)**2+v(2)**2+v(3)**2) if(n.gt.0.0d0)then v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n endif return end subroutine rg(n,x,y,z,zim,i12,ic,maxval) implicit none integer*4 n,ia,maxval,ic,ii real*8 x(*),y(*),z(*) integer*4 zim(*) integer*4 i12(maxval,n) if(ic.eq.1)then do 2 ia=1,n 2 read(10,*)zim(ia),x(ia),y(ia),z(ia),(i12(ii,ia),ii=1,7) else do 1 ia=1,n 1 read(10,*)zim(ia),x(ia),y(ia),z(ia) endif return end subroutine cctn_tinker2(NAT,atomic,x,y,z,mxb, 1R,IND,IB,ND,NDS,RS,FFS,IS,INDBIG,IBONDSB,IBONDS,NBO,NBOB, 2PS,AS,ALPHAS,AAS,GS,AAT,GT, 3asname,psname,aaname,gname,alphaname,rname,NATNAME,ffname, 4NDNAME,IOX,luse) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind,atomic(*) dimension R(3,NAT),IND(NAT),IB(NAT,mxb), 1ND(NAT),NDS(NAT),FFS(3*NAT,3*NAT),IS(NAT,mxb), 1INDBIG(IOX,NAT),IBONDSB(8,NAT),IBONDS(8,NAT),NBO(NAT), 1NBOB(NAT),PS(NAT,3,3),AS(NAT,3,3),ALPHAS(3*NAT,3,3), 3AAS(3*NAT,3,3,3),GS(3*NAT,3,3),AAT(3*NAT,3,3,3),GT(3*NAT,3,3), 1RS(3,NAT), 1asname(IOX,NAT,3,3),psname(IOX,NAT,3,3), 1aaname(IOX,3*NAT,3,3,3),gname(IOX,3*NAT,3,3), 1alphaname(IOX,3*NAT,3,3),rname(IOX,3,NAT),NATNAME(IOX), 1ffname(IOX,3*NAT,3*NAT),NDNAME(IOX,NAT) LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM,lpol logical lxfile CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname DIMENSION x(*),y(*),z(*),POLS(3,3) 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 logical luse(*) c open(3,file='CCT.OUT') NATi=NAT 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 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,NAT,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, 3lpol) 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) read(2,*) read(2,*)NATS if(NATS.gt.NAT)call report('NATS>NAT') rewind 2 CALL READRXS(NAT,NATS,RS,NDS,NBO,IBONDS) 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 OPEN(2,FILE='SMALL.X',STATUS='OLD') read(2,*) read(2,*)NATS if(NATS.gt.NAT)call report('NATS>NAT') rewind 2 CALL READRXS(NAT,NATS,RS,NDS,NBO,IBONDS) write(6,*)NATS,' atoms in SMALL.X' CLOSE(2) endif c names c i12 bond table for the small molecule,rewrite to big for all overlaps c BIG: i -----------ib NBOB c SMALL: INDBIG(1,i)-----------iss NBO do 2011 i=1,NAT 2011 NBOB(i)=0 do 201 i=1,NAT if(luse(i))then do 2012 io=1,NO if(NBOB(i).eq.0)then if(INDBIG(io,i).ne.0)then do 202 j=1,NBO(INDBIG(io,i)) iss=IBONDS(j,INDBIG(io,i)) ibb=0 do 203 ii=1,NAT if(INDBIG(io,ii).eq.iss)then ibb=ii goto 204 endif 203 continue 204 if(ibb.gt.0)then NBOB(i)=NBOB(i)+1 IBONDSB(NBOB(i),i)=ibb endif 202 continue endif c write(6,*)i,(IBONDSB(ii,i),ii=1,NBOB(i)) endif 2012 continue endif 201 continue do 101 i=1,NAT ND(i)=atomic(i) R(1,i)=x(i) R(2,i)=y(i) 101 R(3,i)=z(i) 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) c 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(NAT,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 do 19 ia=1,NATS do 19 ix=1,3 19 RS(ix,ia)=rname(i,ix,ia) c IF(LABS.OR.LVCD)THEN nname=basicname(1:NE)//'.ten' write(6,*)nname(1:NE+4) OPEN(15,FILE=nname,STATUS='OLD') CALL READTEN1(NAT,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(NAT,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(6,*)' A and G obtained from alpha as DO' CALL TTTT1(3*NAT,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*NAT,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(NAT,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(NAT,PS,AS,LAPT,RS,NATS) write(6,*)'SMALL.TEN read' CLOSE(15) ENDIF C IF(LROA.OR.LRAM)THEN OPEN(15,FILE='SMALL.TTT',STATUS='OLD') CALL READTTT(NAT,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO for small' CALL TTTT1(3*NAT,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*NAT,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(NAT,ALPHAS,AAT,GT,NATS) CLOSE(15) ENDIF write(6,*)'Raman tensors read' endif if(lpol)then open(15,file='POL.TTT') call readpol(POLS) write(6,*)'Polarizability read' endif close(15) ENDIF c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn c lnames c write(6,*)'CCTN options read in' close(3) return end subroutine trtena(x,y,z,ALPHA,G,AA,PB,AB,n,mxb, 1R,IND,IB,RS,FF,FFS,IS,INDBIG,IBONDSB,NBOB, 2 PS,AS,ALPHAS,AAS,GS,IOX,AM_c,luse, 3 asname,psname,aaname,gname,alphaname,rname,NATNAME,ffname) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind dimension R(3,n),IND(n),IB(n,mxb),FF(3*n,3*n), 1RS(3,n),FFS(3*n,3*n),IS(n,mxb),INDBIG(IOX,n), 2IBONDSB(8,n),NBOB(n),PS(n,3,3),AS(n,3,3), 4ALPHAS(3*n,3,3),AAS(3*n,3,3,3),GS(3*n,3,3), 5PB(n,3,3),AB(n,3,3),ALPHA(3*n,3,3),AA(3*n,3,3,3), 1G(3*n,3,3),AM_c(n,3,3), 1asname(IOX,n,3,3),psname(IOX,n,3,3), 1aaname(IOX,3*n,3,3,3),gname(IOX,3*n,3,3), 1alphaname(IOX,3*n,3,3),rname(IOX,3,n),NATNAME(IOX), 1ffname(IOX,3*n,3*n) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM,luse(*) DIMENSION 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 c NATi=n c c transcript current geometry: do 101 i=1,NATi R(1,i)=x(i) R(2,i)=y(i) 101 R(3,i)=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,AM_c,INDBIG, 1NO,IBONDSB,NBOB,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,IWG,LROA,LRAM, 1ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT,LNAMES, 3LHALTONERROR,CUTOFF,n1,n2,okind,IERR,LWR,luse,IOX,mxb, 1asname,psname,aaname,gname,alphaname,rname,NATNAME, 1ffname) c if(IERR.ne.0)then write(6,*)'MD point skipped' return endif 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(n,PB,AB,R,NATi) IF(LROA.or.LRAM)THEN CALL TTTT1(3*n,ALPHA,AA,G,NATi,2,R) ENDIF return END SUBROUTINE IMPROVE(FFB,FFS,RB,RS,NAT,NATS,NS,IND,IB,IS,AM,INDBIG, 1NO,IBONDS ,NBOB,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,IERR,LWR,luse,IOX,mxb, 4aname,pname,aaname,gname,alphaname,rname,NATNAME, 1ffname) implicit none INTEGER*4 mxb,IOX INTEGER*4 okind,IND(NAT),IS(NAT,mxb),IB(NAT,mxb),INDBIG(IOX,NAT), 4IBONDS(8,NAT),NBOB(*),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,NF,NO,NS,INDJB,JJX, 4ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,NAT),nats_o(IOX) real*8 FFB(3*NAT,3*NAT),FFS(3*NAT,3*NAT),RB(3,NAT),RS(3,NATS), 1TS(3),TB(3),AM(NAT,3,3),B(3,3),TIJ(3),OTIJ(3),PB(NAT,3,3), 2AB(NAT,3,3),PS(NAT,3,3),AS(NAT,3,3),WF(IOX),ALPHA(3*NAT,3,3), 3ALPHAS(3*NAT,3,3),AA(3*NAT,3,3,3),AAS(3*NAT,3,3,3),G(3*NAT,3,3), 4GS(3*NAT,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),aname(IOX,NAT,3,3),pname(IOX,NAT,3,3), 1aaname(IOX,3*NAT,3,3,3),gname(IOX,3*NAT,3,3), 2alphaname(IOX,3*NAT,3,3),rname(IOX,3,NAT), 1ffname(IOX,3*NAT,3*NAT),A(3,3),RMS LOGICAL LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LNAMES,LHALTONERROR,LWR real*8,allocatable:: XS(:,:),XB(:,:) integer*4 ITU integer*4 natname(IOX) logical luse(*) common/rmssom/RMS c allocate(XS(mxb,3),XB(mxb,3)) IPASS=0 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 if(luse(IA))then DO 201 JA=IA,NAT if(luse(JA))then IF(.NOT.LOFF.AND.IA.NE.JA)GOTO 201 if(LWR)write(6,*)' Pair ',IA,JA 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 201 endif endif IF(okind.ne.2.or.(IA.eq.istart.and.JA.eq.istart))THEN 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.gt.0)then call CRMS(IPASS,NAT,INDBIG,IOX,IO, 1 IA,JA,NBOB,okind,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS,IERR, 1A,XS,XB,ITU,mxb) if(IERR.ne.0)return 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 else c Use distances: DO 26 IX=1,3 26 OTIJ(IX)=OTIJ(IX)/DBLE(NATIO) DIST=(OTIJ(1)-TIJ(1))**2+(OTIJ(2)-TIJ(2))**2 1 +(OTIJ(3)-TIJ(3))**2 IF(DIST.LT.DISTM)THEN DISTM=DIST IOIJ=IO ENDIF WF(NF)=DSQRT(DIST) endif ENDIF 21 CONTINUE if(LWR)write(6,3003)NF,IOIJ 3003 format(i5,' overlaps found; best is # ',i4) 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 cccccccccccccccccccccccccccccccccc ENDIF C C Loop over possible overlaps: ISUM=0 DO 202 IP=1,NF IPO=IOL(IP) IF(IWG.EQ.0)THEN IF(IPO.NE.IOIJ)GOTO 202 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,IPO, 1 IA,JA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS,IERR, 1A,XS,XB,ITU,mxb) if(IERR.ne.0)return 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 202 if(LWR)then write(6,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU) 3004 format('Overlap ',I2,': weighting factor: ',F7.4,';', 1 I4,' atom bed:',/, 2 ' Big :',20I3) write(6,3005)(INDBIG(IPO,IND(I)),I=1,ITU) 3005 format( ' Small:',20I3,/) endif 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) 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 202 CONTINUE endif 201 CONTINUE endif 20 CONTINUE c IA,JA,NF 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,A,XB,XS,mxb,ITU) 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) 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 READRXS(N0,N,R,ND,NB,IBONDS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),NB(N0) integer*4 BT(7),IBONDS(8,N0) READ(2,*) read(2,*)N do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(BT(k),k=1,7) NB(i)=0 do 1 k=1,7 if(BT(k).ne.0)then NB(i)=NB(i)+1 IBONDS(NB(i),i)=BT(k) endif 1 continue 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,A,XB,XS,mxb,ITU) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3),A(3,3), 1XS(mxb,3),XB(mxb,3) do i=1,3 do j=1,3 A(i,j)=0.0d0 enddo enddo IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.1d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG,ITU,A,XS,XB) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG,ITU,A,XS,XB) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG,ITU,A,XS,XB) 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,ITU,A,XS,XB) FTOL=1.0d-7 ITMAX=5000 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(6,*)'mxb',mxb write(6,*)'ITU',ITU write(6,*)1,XS(1,1),XS(1,2),XS(1,3),XB(1,1),XB(1,2),XB(1,3) write(6,*)mxb,XS(mxb,1),XS(mxb,2),XS(mxb,3), 1 XB(mxb,1),XB(mxb,2),XB(mxb,3) write(6,*)'A',A WRITE(6,*)' Rotation has not converged !',ITER 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,ITU,A,XS,XB) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR,ITU,A,XS,XB) 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,ITU,A,XS,XB) 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,ITU,A,XS,XB) 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,n,A,XS,XB) implicit none INTEGER*4 I,n REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,ANG(3),FU, 1A(3,3),XS(n,3),XB(n,3),RMS common/rmssom/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,n 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) DIMENSION P(N0,3,3),A(N0,3,3),X(3,NAT) real*8,allocatable::PV(:,:,:),R(:,:),AI(:,:,:),AJ(:,:,:) LOGICAL LAPT allocate(PV(NAT,3,3),R(3,NAT),AI(NAT,3,3),AJ(NAT,3,3)) 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) DIMENSION P(N0,3,3),A(N0,3,3),X(3,NAT) real*8,allocatable::R(:,:) allocate(R(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 DIMENSION VCD(NAT,3,3),VEL(NAT,3,3),DIP(NAT,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) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT) real*8,allocatable::X(:,:) allocate(X(3,NAT)) 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,n,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,lpol) implicit none INTEGER*4 NS,IB,IS,IND,n,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(n,mxb),IS(n,mxb),IND(n),INDBIG(IOX,n) 1,NDD(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,lpol,loldformat c n1=0 n2=0 okind=1 LRDO=.FALSE. Lpol=.FALSE. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LOFF=.TRUE. LDIA=.TRUE. LHALTONERROR=.TRUE. loldformat=.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.'LPOL') READ(2,*)LPOL 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:7).EQ.'OLDFORM')READ(2,*)loldformat 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 if(loldformat)then 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) else READ(2,*)(INDBIG(I,IAB),IAB=1,NAT) DO 51 IAB=1,NAT 51 INDBIG(I,IAB)=0 READ(2,*)(INDBIG(I,IAB),IAB=1,NAT) endif 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,IPO, 1IA,JA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS,IERR, 1A,XS,XB,ITU,mxb) 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 two bonds away by default c IC=2 ... include all transfered atoms implicit none integer*4 IA,IPASS,II,NAT,INDBIG(IOX,NAT),IBA,IBONDS(8,NAT), 1NBOB(*),JB,IIA,NATS,natname(*),IX,IND(*),iu,IERC,I,IC,IERR, 2IPO,J,JA,ITU,IOX,mxb real*8 RS(3,NAT),RB(3,NAT),TB(3),TS(3),rname(IOX,3,NAT), 1A(3,3),XS(mxb,3),XB(mxb,3) logical LNAMES,LHALTONERROR,LSTRICT,LINV IERC=1 IPASS=0 666 ITU=0 IPASS=IPASS+1 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.or.IC.eq.2)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) write(6,*)'Pair',IA,JA write(6,*)'pass ic',IPASS,IC 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,A,XB,XS,mxb,ITU) 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 c ==================================================================== subroutine protr(x,y,z,ALPHA,PB,N,ic,luse) c project translations/rotations from the electric tensors IMPLICIT none integer*4 I,J,K,IA,IX,IIA,ic,IIX,N,II,nc,jx,kx real*8 x(*),y(*),z(*),PB(N,3,3),ALPHA(3*N,3,3),r(3,N), 1t(3*N),PT(N,3,3),AT(3*N,3,3) real*8,allocatable::A(:,:) logical luse(*) allocate(A(3*N,3*N)) c c N is number of atoms here nc=0 do 1 i=1,N if(luse(i))then nc=nc+1 r(1,nc)=x(i) r(2,nc)=y(i) r(3,nc)=z(i) do 2 ix=1,3 do 2 jx=1,3 PT(nc,ix,jx)=PB(i,ix,jx) do 2 kx=1,3 2 AT(ix+(nc-1)*3,jx,kx)=ALPHA(ix+(i-1)*3,jx,kx) endif 1 continue CALL DOMA(r,A,nc) if(ic.eq.1)call rules(AT,nc,N) c DO 6 J=1,3 DO 11 I=1,3*nc t(I)=0.0d0 DO 11 IIA=1,nc DO 11 IIX=1,3 11 t(I)=t(I)+PT(IIA,IIX,J)*A(I,IIX+3*(IIA-1)) do 6 IA=1,nc do 6 IX=1,3 6 PT(IA,IX,J)=t(IX+3*(IA-1)) c DO 3 J=1,3 DO 3 K=1,3 DO 33 I=1,3*nc t(I)=0.0d0 DO 33 II=1,3*nc 33 t(I)=t(I)+AT(II,J,K)*A(I,II) do 3 I=1,3*nc 3 AT(I,J,K)=t(I) if(ic.eq.1)write(6,*)'projection ok' if(ic.eq.1)call rules(AT,nc,N) nc=0 do 4 i=1,N if(luse(i))then nc=nc+1 do 5 ix=1,3 do 5 jx=1,3 PB(i,ix,jx)=PB(nc,ix,jx) do 5 kx=1,3 5 ALPHA(ix+(i-1)*3,jx,kx)=AT(ix+(nc-1)*3,jx,kx) endif 4 continue RETURN END C ============================================================ SUBROUTINE DOMA(r,A,nat) implicit none integer*4 n,nat,N6,ICOL,I,J,II real*8 r(3,nat),A(3*nat,3*nat),S real*8,allocatable::TEM(:,:) allocate(TEM(3*nat,3*nat)) n=3*nat DO 12 I=1,n DO 12 J=1,n 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)=-r(2,I) A(3*I-1,4)= r(1,I) A(3*I-1,5)=-r(3,I) A(3*I ,5)= r(2,I) A(3*I-2,6)= r(3,I) 1 A(3*I ,6)=-r(1,I) N6=6 CALL ORT(A,n,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.1.0d-10)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 ORT(A,NDIM,N6,NAT) implicit none integer*4 NDIM,N6,NAT,NAT3,I,J,IC,ICP REAL*8 A(NDIM,NDIM),C NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.0.0d0)C=1.0d0/DSQRT(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 C=0.0d0 DO 5 J=1,NAT3 5 C=C+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-C*A(J,ICP) C Normalize : C=0.0d0 DO 7 J=1,NAT3 7 C=C+A(J,IC)**2 IF(C.GT.0.0d0)C=1.0d0/DSQRT(C) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*C RETURN END SUBROUTINE rules(ALPHA,NAT,N0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*N0,3,3),ap(3,3) character*2 xyz(3) data xyz/' X',' Y',' Z'/ c write(6,6077)' Sum rules-alpha:' 6077 format(a17) do 11 ii=1,3 do 21 ix=1,3 do 21 jx=1,3 ap(ix,jx)=0.0d0 do 21 ia=2,NAT 21 ap(ix,jx)=ap(ix,jx)+ALPHA((ia-1)*3+ii,ix,jx) write(6,6017)xyz(ii) 6017 format(a2) write(6,600)(( -ap(ix,jx),jx=1,3),ix=1,3) 600 format(9F8.3) 11 write(6,600)((ALPHA(ii,ix,jx),jx=1,3),ix=1,3) return end subroutine readpol(a) real*8 a(3,3) read(15,*) read(15,*) read(15,*)((a(i,j),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(j,i)=a(i,j) return end c subroutine dc(d,lg,dw,wmin,nw) c predefine convolution function implicit none real*8 dd,ff,w,dw,wmin,pi,spi,d character*1 lg integer*4 nw0,i,nc,nw parameter (nw0=2000) common/convc/ff(2*nw0+1),nc pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) if(nw.gt.nw0)call report('nw too big in dc') w=wmin-dw if(lg.eq.'L')then do 1 nc=1,nw w=w+dw dd=((w-wmin)/d)**2 ff(nc)=1.0d0/d/(dd+1.0d0)/pi 1 if(dd.gt.1.0d3)goto 2 else do 3 nc=1,nw w=w+dw dd=((w-wmin)/d)**2 ff(nc)=exp(-dd)/d/spi 3 if(dd.gt.20.0d0)goto 2 endif 2 do 4 i=1,nc-1 4 ff(nc+i)=ff(i+1) ff(nc)=ff(1) do 5 i=1,nc-1 5 ff(nc-i)=ff(nc+i) return end subroutine trabc(ia,xm,ym,zm,a,ss) integer*4 ia real*8 xm(*),ym(*),zm(*),a(*),ss xm(ia)=xm(ia)+a(1)*ss ym(ia)=ym(ia)+a(2)*ss zm(ia)=zm(ia)+a(3)*ss return end subroutine make3v(v1,v2,v3) c makes system of orthonormal vectors from v1 and v2 implicit none real*8 v1(*),v2(*),v3(*) v3(1)=v1(2)*v2(3)-v1(3)*v2(2) v3(2)=v1(3)*v2(1)-v1(1)*v2(3) v3(3)=v1(1)*v2(2)-v1(2)*v2(1) call norm(v1) call norm(v3) v2(1)=v3(2)*v1(3)-v3(3)*v1(2) v2(2)=v3(3)*v1(1)-v3(1)*v1(3) v2(3)=v3(1)*v1(2)-v3(2)*v1(1) call norm(v2) return end subroutine readopt1(K0,n,NSA,iprint,nmol,nnames,isr,natm, 1imk,al, 1dt,kelvin,bohr,pi,dbb,perx,pery,perz,tp,sol,dw,wmin,wmax,res, 1ac,bc,cc, 2lper,lcell,lp27,lp28,lp29,lproject,lnnames,lpair,LABS,LVCD,LRAM, 2LROA,ltti,ldress,ltei,lani,lconv,lg,snam,lincr) implicit none integer*4 K0,n,NSA,iprint,nmol,nnames,ia,isr(*),natm(*), 1imk(*),al(nmol,n),i,mm real*8 dt,kelvin,bohr,pi,dbb,perx,pery,perz,tp,sol,dw,wmin,wmax, 1res,ac(*),bc(*),cc(*) logical lper,lcell,lp27,lp28,lp29,lproject,lnnames,lpair,LABS, 1LVCD,LRAM,LROA,ltti,ldress,ltei,lani,lconv,ldyn,lincr character*1 lg character*10 s10,snam(*) NSA=0 iprint=100 nmol=0 dt=0.00009672d0 kelvin=298.15d0 bohr=0.529177d0 pi=4.0d0*datan(1.0d0) c au^2 -> debye^2 dbb=2.541765d0**2 perx=0.0d0 pery=0.0d0 perz=0.0d0 tp=2.0d0*pi sol=3.0d10 c dw,wmin,wmax step and frequency limits in cm-1 dw=1.0d0 wmin=0.0d0 wmax=4000.0d0 res=5.0d0 lper=.false. lcell=.false. lp27=.false. lp28=.false. lp29=.false. lproject=.false. lnnames=.false. lpair=.false. LABS=.false. LVCD=.false. LRAM=.false. LROA=.false. ltti=.false. ldress=.false. ltei=.false. lani=.false. lconv=.false. lincr=.false. lg='G' inquire(file='DYNAMIM.OPT',exist=ldyn) if(.not.ldyn)call report('DYNAMIM.OPT not found') open(10,file='DYNAMIM.OPT') 311 read(10,71,end=3111,err=3111)s10 71 format(a10) if(s10(1:4).eq.'NRES')then read(10,*)NSA read(10,*)(isr(ia),ia=1,NSA) endif if(s10(1:4).eq.'TSTE')read(10,*)dt if(s10(1:4).eq.'IPRI')read(10,*)iprint if(s10(1:2).eq.'DW')read(10,*)dw if(s10(1:3).eq.'RES')read(10,*)res if(s10(1:4).eq.'WMIN')read(10,*)wmin if(s10(1:4).eq.'WMAX')read(10,*)wmax if(s10(1:4).eq.'TTTI')read(10,*)ltti if(s10(1:4).eq.'TEMK')read(10,*)kelvin if(s10(1:4).eq.'DRES')read(10,*)ldress if(s10(1:4).eq.'TENI')read(10,*)ltei if(s10(1:4).eq.'ANIM')read(10,*)lani if(s10(1:4).eq.'CONV')read(10,*)lconv if(s10(1:4).eq.'INCR')read(10,*)lincr if(s10(1:4).eq.'CONT')read(10,*)lg if(s10(1:4).eq.'PROJ')read(10,*)lproject if(s10(1:4).eq.'PERI')read(10,*)lper,perx,pery,perz if(s10(1:4).eq.'CELL')then read(10,*)lcell do 72 i=1,3 72 read(10,*)ac(i),bc(i),cc(i) endif if(s10(1:3).eq.'P27')read(10,*)lp27 if(s10(1:3).eq.'P28')read(10,*)lp28 if(s10(1:3).eq.'P29')read(10,*)lp29 if(s10(1:4).eq.'PAIR')read(10,*)lpair if(s10(1:4).eq.'NMOL')then read(10,*)nmol do 191 mm=1,nmol 191 read(10,*)natm(mm),(al(mm,ia),ia=1,natm(mm)) endif c c named tensor files: names and molecular assignment imk must be given: if(s10(1:6).eq.'NNAMES')then read(10,*)lnnames if(lnnames)then read(10,*)nnames if(nnames.gt.K0)call report('too many molecular kinds') do 1911 i=1,nnames 1911 read(10,71)snam(i) read(10,*)nmol do 160 mm=1,nmol 160 read(10,*)imk(mm) endif endif goto 311 c NSA - number of atoms for spectral response (followed by their list) c TSTE - time step in picoseconds c lper ... correct for periodic conditions given by perx pery perz (A) c lcell ... non-orthogonal crystall cell c lp27 ... calculate with the 3^3 large periodic box c lp28 ... damp towards walls of the smaller box c lp29 ... real superbox c lincr ... restart accumulation of spectra from the disk 3111 close(10) return end subroutine readopt3(nmol,lpair) integer*4 nmol logical ldyn,lpair character*10 s10 lpair=.false. nmol=0 inquire(file='DYNAMIM.OPT',exist=ldyn) if(.not.ldyn)call report('DYNAMIM.OPT not found') open(10,file='DYNAMIM.OPT') 311 read(10,71,end=3111,err=3111)s10 71 format(a10) if(s10(1:4).eq.'NMOL')read(10,*)nmol if(s10(1:4).eq.'PAIR')read(10,*)lpair goto 311 3111 close(10) return end