program dynamiX c "dynamicX" read MD trajectory from FILE.X and computes spectra c time correlation functions, spectra from FT transformations: implicit none integer*4 maxval,nat0,nw0 parameter (nat0=300,nw0=10000,maxval=8) integer*4 ia,iprintft,iw,iprint,iwrite,n,nw,ix,nat,NAT1, 1LP,i12(maxval,nat0),istep,iii(4,nat0),i,iind,ip,ixp,j,jp,jx,k, 2kp,kx,l,atomic(nat0),NS,n1,n2,okind,NO,NATS,IWG,isr(nat0),nsa character*10 s10 logical lexc,lres,lqftry,ltti,ltei,lexp,ldyn,lexx,lper, 1LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,lcct,luse(nat0), 1lproject common/usec/luse real*8 x(nat0),y(nat0),z(nat0),xdum,dt, 1xm(nat0),ym(nat0),zm(nat0),xmm(nat0),ymm(nat0),zmm(nat0), 2tsux(nw0),tsuy(nw0),tsuz(nw0),tsmx(nw0),tsmy(nw0),tsmz(nw0), 3tcux(nw0),tcuy(nw0),tcuz(nw0),tcmx(nw0),tcmy(nw0),tcmz(nw0), 2wsux(nw0),wsuy(nw0),wsuz(nw0),wsmx(nw0),wsmy(nw0),wsmz(nw0), 3wcux(nw0),wcuy(nw0),wcuz(nw0),wcmx(nw0),wcmy(nw0),wcmz(nw0), 5spec(nw0),stem(nw0),t,w,dbb,ramconst,dw,wmax,tp,sol,wmin,res,sc, 6fac,FBOHR,dtau,bohr,pi,pchgft(nat0),pol(nat0), 2kelvin,P(nat0,3,3),A(nat0,3,3),vx,vy,vz,ux, 4ALPHAI(3*nat0,3,3),AAI(3*nat0,3,3,3),GI(3*nat0,3,3), 4ALPHA(3*nat0,3,3), AA(3*nat0,3,3,3), G(3*nat0,3,3), 5APTI(nat0,3,3),AI(nat0,3,3),u(3,3,nat0),Xwork(3,nat0), 2wsat(3,3,nw0),wcat(3,3,nw0),wsgt(3,3,nw0),wcgt(3,3,nw0), 2wsaat(3,3,3,nw0),wcaat(3,3,3,nw0),sa1,sal0,sal1,sag0,sag1, 1s180,d180,ydy,ydx,tf,at(3,3),gt(3,3),aat(3,3,3),s1,s2, 1at1(3,3),at2(3,3),pt1(3,3),pt2(3,3),ut(3,3),rki(3),rd, 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,perx2,pery2,perz2 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 NSA=0 dt=0.00009672d0 iprint=100 kelvin=273.15d0 bohr=0.529177d0 pi=4.0d0*datan(1.0d0) c au^2 -> debye^2 dbb=2.541765d0**2 lres=.false. lper=.false. perx=0.0d0 pery=0.0d0 perz=0.0d0 lqftry=.true. tp=2.0d0*pi sol=3.0d10 iprintft=1 c dw,wmin,wmax step and frequency limits in cm-1 dw=1.0d0 wmin=0.0d0 wmax=4000.0d0 res=5.0d0 ltti=.false. lproject=.false. ltei=.false. inquire(file='DYNAMIX.OPT',exist=ldyn) if(.not.ldyn)call report('DYNAMIX.OPT not found') open(10,file='DYNAMIX.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.'LQFT')read(10,*)lqftry if(s10(1:4).eq.'TTTI')read(10,*)ltti if(s10(1:4).eq.'TENI')read(10,*)ltei if(s10(1:2).eq.'IP')read(10,*)iprintft if(s10(1:4).eq.'CONT')read(10,*)lres if(s10(1:4).eq.'PROJ')read(10,*)lproject if(s10(1:4).eq.'PERI')then read(10,*)lper read(10,*)perx,pery,perz 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 prey perz (A) 3111 close(10) perx2=perx/2.0d0 pery2=pery/2.0d0 perz2=perz/2.0d0 if(lper)write(6,6000)perx,pery,perz 6000 format(' Periodic box: ',3f10.2) inquire(file='BOND.X',exist=lexx) if(lexx)then open(10,file='BOND.X') read(10,*) read(10,*)n if(n.gt.nat0)call report('too many atoms for FT') call rg(n,x,y,z,atomic,i12,1,nat0,maxval) close(10) write(6,*)'Bonds read from BOND.X' endif inquire(file='FILE.X',exist=lexx) if(.not.lexx)call report('FILE.X not found') open(10,file='FILE.X') read(10,*) read(10,*)n if(n.gt.nat0)call report('too many atoms for FT') call rg(n,x,y,z,atomic,i12,0,nat0,maxval) close(10) NAT=n 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,6013)iprint,dt,n,NSA 6013 format('DYNAMIX.OPT:',/,' printing interval : ',i5,'steps',/, 1' dt : ',g15.6,' ps',/, 1' n : ',i5,' atoms',/, 1' NSA : ',i5,' atoms',/) 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,80).eq.0)write(6,*) enddo write(6,*) iwrite=iprint c inquire(file='CHARGE.LST',exist=lexc) if(lexc)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) endif if(ltti)then open(15,FILE='INT.TTT',status='old') CALL READTTT(nat0,ALPHAI,AAI,GI,NAT1) call readatms(15,NAT1,iii,nat0) close(15) if(NAT1.ne.n)call report('Strange number of atoms in INT.TTT') endif if(ltei)then OPEN(15,FILE='INT.TEN',status='old') CALL READTENI(nat0,APTI,AI,N) call readatms(15,N,iii,nat0) CLOSE(15) endif inquire(file='FTRY.INP',exist=lexc) if(lexc)then open(99,file='FTRY.INP') read(99,*) read(99,*) if(lqftry)then read(99,6009)(pchgft(ia),ia=1,n) 6009 format(6f12.6) write(6,*)'partial charges read from FTRY.INP' else read(99,6009)(xdum,ia=1,n) endif read(99,*) read(99,6009)(xdum,ia=1,n) close(99) endif dw=dw*tp*sol wmax=wmax*tp*sol wmin=wmin*tp*sol 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 if(iprintft.ge.1)then write(6,*)nw,' frequencies' write(6,6700)sc 6700 format(' Scaling factor ',f9.6) write(6,6701)dt,dtau 6701 format(f10.5,' ps step = ',f10.4,' au') c time step in atomic units: write(6,6702)res 6702 format(f10.5,' cm-1 resolution via convolution') endif c c initiate frequency coefficients c zero out first to avoid compiler messages: call ifc(nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux, 1wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux, 3tsuy,tcuy,tsuz,tcuz,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0) if(lres)then call rwfc(1,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux, 1 wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux, 3 tsuy,tcuy,tsuz,tcuz,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0) write(6,*)'Restart from disk' else write(6,*)'Coefficients zeroed-out' 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)call cctn_tinker2(n,atomic,x,y,z) write(6,*)'ltei ',ltei write(6,*)'ltti ',ltti write(6,*)'lexp ',lexp write(6,*)'LRAM ',LRAM write(6,*)'LROA ',LROA write(6,*)'LABS ',LABS write(6,*)'LVCD ',LVCD write(6,*)'lproject ',lproject c ######################################################## open(10,file='FILE.X') istep=0 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,nat0,maxval) if(lcct)then call trtena(x,y,z,ALPHA,G,AA,P,A,n) c project translations and rotations from alpha if(lproject)call protr(x,y,z,ALPHA,G,AA,P,A,n) if(istep.eq.1)write(6,*)'Tensors by transfer' else c internal coordinate tensors, define the transformation c matrices: if(ltei.or.ltti)call dtm(iii,n,u,nat0,x,y,z) if(ltei)then DO 14 L=1,NAT DO 1411 IP=1,3 DO 1411 J=1,3 1411 ut(IP,J)=u(IP,J,L) DO 141 IP=1,3 DO 141 JP=1,3 at1(IP,JP)=AI(L,IP,JP) 141 pt1(IP,JP)=APTI(L,IP,JP) 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) DO 14 I=1,3 DO 14 J=1,3 s1=0.0d0 s2=0.0d0 do 143 JP=1,3 s1=s1+ut(JP,J)*at2(I,JP) 143 s2=s2+ut(JP,J)*pt2(I,JP) A(L,I,J)=s1 14 P(L,I,J)=s2 c c Transform to common origin: FBOHR=0.25d0/bohr DO 2201 L=1,NAT DO 2201 J=1,3 A(L,J,1)=A(L,J,1)+(y(L)*P(L,J,3)-z(L)*P(L,J,2))*FBOHR A(L,J,2)=A(L,J,2)+(z(L)*P(L,J,1)-x(L)*P(L,J,3))*FBOHR 2201 A(L,J,3)=A(L,J,3)+(x(L)*P(L,J,2)-y(L)*P(L,J,1))*FBOHR if(istep.eq.1)write(6,*)'DD from internals' else if(istep.eq.1)write(6,*)'DD from charges' do 401 ia=1,n fac=pchgft(ia) do 402 ix=1,3 do 402 jx=1,3 P(ia,ix,jx)=0.0d0 402 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 c polarizability from atomic polarizabilities DO 179 L=1,n if(luse(L))then DO 172 IX=1,3 IIND=3*(L-1)+IX DO 172 I=1,3 DO 172 J=1,3 DO 177 K=1,3 177 AA(IIND,I,J,K)=0.0d0 G(IIND,I,J)=0.0d0 172 ALPHA(IIND,I,J)=0.0d0 endif 179 continue DO 17 L=1,n if(luse(L))then DO 178 LP=1,n 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 rd=dsqrt(rki(1)**2+rki(2)**2+rki(3)**2+1.0d-9) 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-rd*rd*rki(J) if(J.eq.IX)s1=s1-rd*rd*rki(I) if(I.eq.J)s1=s1-rd*rd*rki(IX) 173 ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+6.0d0*pol(L)*pol(LP)*s1/rd**7 endif endif 178 continue endif 17 continue CALL TTTT(3*nat0,ALPHA,AA,G,NAT,2,x,y,z,Xwork,nat0) if(istep.eq.1)write(6,*)'Alpha from atomic polarizabilities' else c lesp if(ltti)then DO 15 L=1,n c rewrite DO 1412 IP=1,3 DO 1412 J=1,3 1412 ut(IP,J)=u(IP,J,L) DO 151 IXP=1,3 IIND=3*(L-1)+IXP DO 151 IP=1,3 DO 151 JP=1,3 as1(IXP,IP,JP)=ALPHAI(IIND,IP,JP) 151 gs1(IXP,IP,JP)=G(IIND,IP,JP) c first index trafo DO 154 IX=1,3 DO 154 IP=1,3 DO 154 JP=1,3 as2(IX,IP,JP)=0.0d0 gs2(IX,IP,JP)=0.0d0 DO 154 IXP=1,3 as2(IX,IP,JP)=as2(IX,IP,JP)+ut(IXP,IX)*as1(IXP,IP,JP) 154 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 ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 DO 153 JP=1,3 ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+ut(JP,J)*as1(IX,I,JP) 153 G(IIND,I,J) =G(IIND,I,J) +ut(JP,J)*gs1(IX,I,JP) c A-rewrite DO 16 IXP=1,3 IIND=3*(L-1)+IXP DO 16 IP=1,3 DO 16 JP=1,3 DO 16 KP=1,3 16 au1(IXP,IP,JP,KP)=AAI(IIND,IP,JP,KP) 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 AA(IIND,I,J,K)=0.0d0 DO 15 KP=1,3 15 AA(IIND,I,J,K)=AA(IIND,I,J,K)+ut(KP,K)*au2(IX,I,J,KP) CALL TTTT(3*nat0,ALPHA,AA,G,NAT,2,x,y,z,Xwork,nat0) if(istep.eq.1)write(6,*)'Alpha from internal tensors' else c ltti call report('Alpha not defined') endif c ltti endif c lesp endif c if lcct c correct for periodic jumps: if(lper)then do 180 ia=1,n if(x(ia)-xm(ia).gt.perx2)then xm(ia)=xm(ia)+perx xmm(ia)=xmm(ia)+perx endif if(xm(ia)-x(ia).gt.perx2)then xm(ia)=xm(ia)-perx xmm(ia)=xmm(ia)-perx endif if(y(ia)-ym(ia).gt.pery2)then ym(ia)=ym(ia)+pery ymm(ia)=ymm(ia)+pery endif if(ym(ia)-y(ia).gt.pery2)then ym(ia)=ym(ia)-pery ymm(ia)=ymm(ia)-pery endif if(z(ia)-zm(ia).gt.perz2)then zm(ia)=zm(ia)+perz zmm(ia)=zmm(ia)+perz endif if(zm(ia)-z(ia).gt.perz2)then zm(ia)=zm(ia)-perz zmm(ia)=zmm(ia)-perz endif 180 continue endif c time-dependent polarizability derivatives: do 2001 ix=1,3 do 2001 jx=1,3 at(ix,jx)=0.0d0 gt(ix,jx)=0.0d0 do 2001 kx=1,3 2001 aat(ix,jx,kx)=0.0d0 fac=1.0d0/dtau/bohr do 2002 ia=1,n if(luse(ia))then vx=(x(ia)-xm(ia))*fac vy=(y(ia)-ym(ia))*fac vz=(z(ia)-zm(ia))*fac do 20021 ix=1,3 do 20021 jx=1,3 gt(ix,jx)=gt(ix,jx) +G(3*(ia-1)+1,ix,jx)*vx 1 + G(3*(ia-1)+2,ix,jx)*vy+G(3*(ia-1)+3,ix,jx)*vz at(ix,jx)=at(ix,jx) +ALPHA(3*(ia-1)+1,ix,jx)*vx 1 +ALPHA(3*(ia-1)+2,ix,jx)*vy+ALPHA(3*(ia-1)+3,ix,jx)*vz do 20021 kx=1,3 20021 aat(ix,jx,kx)=aat(ix,jx,kx) +AA(3*(ia-1)+1,ix,jx,kx)*vx 1 +AA(3*(ia-1)+2,ix,jx,kx)*vy+AA(3*(ia-1)+3,ix,jx,kx)*vz endif 2002 continue c A) selected atoms c update dipoles: call updt(1,n,wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz, 1wcmz,wsux,wcux,wsuy,wcuy,wsuz,wcuz,t,sc,xm,ym,zm,xmm,ymm,zmm, 2bohr,dtau,x,y,z,A,P,nat0,ux,luse) c update poalrizabilities: if(ltti.or.lexp.or.LROA.or.LRAM)call update3(wmin,dw,nw,wsat, 1wcat,wsgt,wcgt,wsaat,wcaat,at,gt,aat,t,sc,3*nat0) c c B) all atoms c update dipoles: call updt(1,n,wmin,dw,nw,tsmx,tcmx,tsmy,tcmy,tsmz, 1tcmz,tsux,tcux,tsuy,tcuy,tsuz,tcuz,t,sc,xm,ym,zm,xmm,ymm,zmm, 2bohr,dtau,x,y,z,A,P,nat0,ux,luse) 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 c write Fourier coefficients for restart: if(mod(istep,iwrite).eq.0)then call rwfc(0,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux, 1 wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux, 3 tsuy,tcuy,tsuz,tcuz, 4 wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0) c c selected atoms c write the normalized power spectrum do 711 iw=1,nw w=wmin+dble(iw-1)*dw fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol 711 spec(iw)=fac* 1(wsux(iw)*wsux(iw)+wsuy(iw)*wsuy(iw)+wsuz(iw)*wsuz(iw) 2+wcux(iw)*wcux(iw)+wcuy(iw)*wcuy(iw)+wcuz(iw)*wcuz(iw)) c c convolute it with a smoothing function and write on disk call wrspec('WDraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WD.PRN',nw,wmin,dw,spec,tp,sol) do 8 iw=1,nw w=wmin+dble(iw-1)*dw fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol 8 spec(iw)=fac* 1(wsux(iw)*wcmx(iw)+wsuy(iw)*wcmy(iw)+wsuz(iw)*wcmz(iw) 2-wcux(iw)*wsmx(iw)-wcuy(iw)*wsmy(iw)-wcuz(iw)*wsmz(iw)) c call wrspec('WRraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WR.PRN',nw,wmin,dw,spec,tp,sol) c whole system response: do 9 iw=1,nw w=wmin+dble(iw-1)*dw fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol 9 spec(iw)=fac* 1(tsux(iw)*tsux(iw)+tsuy(iw)*tsuy(iw)+tsuz(iw)*tsuz(iw) 2+tcux(iw)*tcux(iw)+tcuy(iw)*tcuy(iw)+tcuz(iw)*tcuz(iw)) call wrspec('WDTraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WDT.PRN',nw,wmin,dw,spec,tp,sol) do 11 iw=1,nw w=wmin+dble(iw-1)*dw fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol 11 spec(iw)=fac* 1(tsux(iw)*tcmx(iw)+tsuy(iw)*tcmy(iw)+tsuz(iw)*tcmz(iw) 2-tcux(iw)*tsmx(iw)-tcuy(iw)*tsmy(iw)-tcuz(iw)*tsmz(iw)) call wrspec('WRTraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WRT.PRN',nw,wmin,dw,spec,tp,sol) c c Raman-selected atoms only: if(ltti.or.lexp.or.LROA.or.LRAM)then open(51,file='RAM.SCR',form='unformatted') do 1233 iw=1,nw w=wmin+dble(iw-1)*dw 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/298.0d0))*1.0d5 endif if(lexp)then tf=tf*w*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 sa1=sa1+wsat(1,jx,iw)*(wsaat(2,3,jx,iw)-wsaat(3,2,jx,iw)) 1 +wcat(1,jx,iw)*(wcaat(2,3,jx,iw)-wcaat(3,2,jx,iw)) 2 +wsat(2,jx,iw)*(wsaat(3,1,jx,iw)-wsaat(1,3,jx,iw)) 3 +wcat(2,jx,iw)*(wcaat(3,1,jx,iw)-wcaat(1,3,jx,iw)) 4 +wsat(3,jx,iw)*(wsaat(1,2,jx,iw)-wsaat(2,1,jx,iw)) 5 +wcat(3,jx,iw)*(wcaat(1,2,jx,iw)-wcaat(2,1,jx,iw)) do 121 ix=1,3 sal0=sal0+wsat(ix,ix,iw)*wsat(jx,jx,iw) 1 +wcat(ix,ix,iw)*wcat(jx,jx,iw) sal1=sal1+wsat(ix,jx,iw)*wsat(ix,jx,iw) 1 +wcat(ix,jx,iw)*wcat(ix,jx,iw) c coordinates sag0=sag0+wsat(ix,ix,iw)*wsgt(jx,jx,iw) 1 +wcat(ix,ix,iw)*wcgt(jx,jx,iw) sag1=sag1+wsat(ix,jx,iw)*wsgt(ix,jx,iw) 1 +wcat(ix,jx,iw)*wcgt(ix,jx,iw) c momenta: c sag0=sag0+wsat(ix,ix,iw)*wcgt(jx,jx,iw) c 1 -wcat(ix,ix,iw)*wsgt(jx,jx,iw) c sag1=sag1+wsat(ix,jx,iw)*wcgt(ix,jx,iw) c 1 -wcat(ix,jx,iw)*wsgt(ix,jx,iw) 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) 1233 write(51)s180*tf,ydy*tf,ydx*tf,d180*tf rewind 51 ramconst=(w/sol/1000.0d0)**2/dble(istep)**2 do 1234 iw=1,nw read(51)s180,ydy,ydx,d180 1234 spec(iw)=ramconst*s180 call wrspec('WRam180raw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WRam180.PRN',nw,wmin,dw,spec,tp,sol) rewind 51 do 1235 iw=1,nw read(51)s180,ydy,ydx,d180 1235 spec(iw)=ramconst*ydy call wrspec('WRamyraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WRamy.PRN',nw,wmin,dw,spec,tp,sol) rewind 51 do 1236 iw=1,nw read(51)s180,ydy,ydx,d180 1236 spec(iw)=ramconst*ydx call wrspec('WRamxraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WRamx.PRN',nw,wmin,dw,spec,tp,sol) rewind 51 do 1237 iw=1,nw read(51)s180,ydy,ydx,d180 1237 spec(iw)=ramconst*d180 call wrspec('WRoa180raw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec('WRoa180.PRN',nw,wmin,dw,spec,tp,sol) close(51) endif endif c ######################################################## goto 31111 31112 close(10) write(6,6002)istep 6002 format(i10,' geometries in FILE.X') 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(sn,nw,wmin,dw,spec,tp,sol) implicit none character*(*) sn real*8 wmin,dw,spec(*),w,tp,sol integer*4 nw,iw open(8,file=sn) 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,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz, 1wsux,wcux,wsuy,wcuy,wsuz,wcuz,mx,my,mz,ux,uy,uz,t,sc) implicit none integer*4 nw,iw real*8 wmin,w,mx,my,mz,ux,uy,uz,t,ss,cc,sc,dw,tp,sol, 2wsux(*),wsuy(*),wsuz(*),wsmx(*),wsmy(*),wsmz(*), 3wcux(*),wcuy(*),wcuz(*),wcmx(*),wcmy(*),wcmz(*),wau w=wmin-dw tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 if(abs(w).gt.1.0d-10)then ss=-cos(w*t)/wau cc= sin(w*t)/wau wsmx(iw)=wsmx(iw)*sc+mx*ss wcmx(iw)=wcmx(iw)*sc+mx*cc wsmy(iw)=wsmy(iw)*sc+my*ss wcmy(iw)=wcmy(iw)*sc+my*cc wsmz(iw)=wsmz(iw)*sc+mz*ss wcmz(iw)=wcmz(iw)*sc+mz*cc wsux(iw)=wsux(iw)*sc+ux*ss wcux(iw)=wcux(iw)*sc+ux*cc wsuy(iw)=wsuy(iw)*sc+uy*ss wcuy(iw)=wcuy(iw)*sc+uy*cc wsuz(iw)=wsuz(iw)*sc+uz*ss wcuz(iw)=wcuz(iw)*sc+uz*cc endif 1 continue return end subroutine update3(wmin,dw,nw,wsat,wcat,wsgt,wcgt, 1wsaat,wcaat,at,gt,aat,t,sc,MX3) c alpha,A,G -> Raman, ROA implicit none integer*4 nw,iw,ix,jx,kx,MX3 real*8 wmin,w,at(3,3),gt(3,3),aat(3,3,3),t,ss,cc,sc,dw,tp,sol, 2wsat(3,3,MX3),wcat(3,3,MX3),wsgt(3,3,MX3),wcgt(3,3,MX3), 2wsaat(3,3,3,MX3),wcaat(3,3,3,MX3),wau 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 if(abs(w).gt.1.0d-10)then ss=-cos(w*t)/wau cc= sin(w*t)/wau do 2 ix=1,3 do 2 jx=1,3 wsat(ix,jx,iw)=wsat(ix,jx,iw)*sc+at(ix,jx)*ss wcat(ix,jx,iw)=wcat(ix,jx,iw)*sc+at(ix,jx)*cc wsgt(ix,jx,iw)=wsgt(ix,jx,iw)*sc+gt(ix,jx)*ss wcgt(ix,jx,iw)=wcgt(ix,jx,iw)*sc+gt(ix,jx)*cc do 2 kx=1,3 wsaat(ix,jx,kx,iw)=wsaat(ix,jx,kx,iw)*sc+aat(ix,jx,kx)*ss 2 wcaat(ix,jx,kx,iw)=wcaat(ix,jx,kx,iw)*sc+aat(ix,jx,kx)*cc endif 1 continue return end subroutine ifc(nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux, 1wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux, 3tsuy,tcuy,tsuz,tcuz, 4 wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0) implicit none integer iw,nw,ix,jx,kx,nw0 real*8 wsmx(*),wcmx(*),wsmy(*),wcmy(*),wsmz(*),wcmz(*), 1wsux(*),wcux(*),wsuy(*),wcuy(*),wsuz(*),wcuz(*),tsmx(*), 2tcmx(*),tsmy(*),tcmy(*),tsmz(*),tcmz(*),tsux(*),tcux(*), 3tsuy(*),tcuy(*),tsuz(*),tcuz(*), 2wsat(3,3,nw0),wcat(3,3,nw0),wsgt(3,3,nw0),wcgt(3,3,nw0), 2wsaat(3,3,3,nw0),wcaat(3,3,3,nw0) do 1 iw=1,nw wsmx(iw)=0.0d0 wcmx(iw)=0.0d0 wsmy(iw)=0.0d0 wcmy(iw)=0.0d0 wsmz(iw)=0.0d0 wcmz(iw)=0.0d0 wsux(iw)=0.0d0 wcux(iw)=0.0d0 wsuy(iw)=0.0d0 wcuy(iw)=0.0d0 wsuz(iw)=0.0d0 wcuz(iw)=0.0d0 tsmx(iw)=0.0d0 tcmx(iw)=0.0d0 tsmy(iw)=0.0d0 tcmy(iw)=0.0d0 tsmz(iw)=0.0d0 tcmz(iw)=0.0d0 tsux(iw)=0.0d0 tcux(iw)=0.0d0 tsuy(iw)=0.0d0 tcuy(iw)=0.0d0 tsuz(iw)=0.0d0 tcuz(iw)=0.0d0 do 1 ix=1,3 do 1 jx=1,3 wsat(ix,jx,iw)=0.0d0 wcat(ix,jx,iw)=0.0d0 wsgt(ix,jx,iw)=0.0d0 wcgt(ix,jx,iw)=0.0d0 do 1 kx=1,3 wcaat(ix,jx,kx,iw)=0.0d0 1 wsaat(ix,jx,kx,iw)=0.0d0 return end subroutine rwfc(ic,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux, 1wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux, 3tsuy,tcuy,tsuz,tcuz, 4 wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0) implicit none integer iw,nw,ic,nw0,ix,jx,kx real*8 wsmx(*),wcmx(*),wsmy(*),wcmy(*),wsmz(*),wcmz(*), 1wsux(*),wcux(*),wsuy(*),wcuy(*),wsuz(*),wcuz(*),tsmx(*), 2tcmx(*),tsmy(*),tcmy(*),tsmz(*),tcmz(*),tsux(*),tcux(*), 3tsuy(*),tcuy(*),tsuz(*),tcuz(*), 2wsat(3,3,nw0),wcat(3,3,nw0),wsgt(3,3,nw0),wcgt(3,3,nw0), 2wsaat(3,3,3,nw0),wcaat(3,3,3,nw0) open(9,file='FIT.SCR',form='unformatted') if(ic.eq.0)then write(9) 1(wsmx(iw),iw=1,nw),(wsmy(iw),iw=1,nw),(wsmz(iw),iw=1,nw), 1(wcmx(iw),iw=1,nw),(wcmy(iw),iw=1,nw),(wcmz(iw),iw=1,nw), 2(wsux(iw),iw=1,nw),(wsuy(iw),iw=1,nw),(wsuz(iw),iw=1,nw), 2(wcux(iw),iw=1,nw),(wcuy(iw),iw=1,nw),(wcuz(iw),iw=1,nw), 3(tsmx(iw),iw=1,nw),(tsmy(iw),iw=1,nw),(tsmz(iw),iw=1,nw), 3(tcmx(iw),iw=1,nw),(tcmy(iw),iw=1,nw),(tcmz(iw),iw=1,nw), 4(tsux(iw),iw=1,nw),(tsuy(iw),iw=1,nw),(tsuz(iw),iw=1,nw), 4(tcux(iw),iw=1,nw),(tcuy(iw),iw=1,nw),(tcuz(iw),iw=1,nw), 4(((wsat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4(((wcat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4(((wsgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4(((wcgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4((((wsaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw), 4((((wcaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw) else read(9) 1(wsmx(iw),iw=1,nw),(wsmy(iw),iw=1,nw),(wsmz(iw),iw=1,nw), 1(wcmx(iw),iw=1,nw),(wcmy(iw),iw=1,nw),(wcmz(iw),iw=1,nw), 2(wsux(iw),iw=1,nw),(wsuy(iw),iw=1,nw),(wsuz(iw),iw=1,nw), 2(wcux(iw),iw=1,nw),(wcuy(iw),iw=1,nw),(wcuz(iw),iw=1,nw), 3(tsmx(iw),iw=1,nw),(tsmy(iw),iw=1,nw),(tsmz(iw),iw=1,nw), 3(tcmx(iw),iw=1,nw),(tcmy(iw),iw=1,nw),(tcmz(iw),iw=1,nw), 4(tsux(iw),iw=1,nw),(tsuy(iw),iw=1,nw),(tsuz(iw),iw=1,nw), 4(tcux(iw),iw=1,nw),(tcuy(iw),iw=1,nw),(tcuz(iw),iw=1,nw), 4(((wsat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4(((wcat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4(((wsgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4(((wcgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw), 4((((wsaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw), 4((((wcaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw) endif close(9) return end SUBROUTINE READTTT(MX,ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,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,*)LDUM,LDUM,(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,*)LDUM,LDUM,(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,*)LDUM,LDUM,(A(IIND,I,J,K),K=1,3) RETURN END SUBROUTINE READTENI(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (nat0=300) DIMENSION P(N0,3,3),A(N0,3,3),AI(nat0,3,3),AJ(nat0,3,3) READ (15,*) NOAT if(NOAT.ne.NAT)call report('NAT<>NOAT') DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) 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,*) 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 RETURN END subroutine readatms(io,NAT,iii,MX) integer*4 MX,NAT,iii(4,MX),ia,io,idumm do 1 ia=1,NAT 1 read(io,*)idumm,(iii(ii,ia),ii=1,4) return end subroutine dtm(ii,NAT,u,MX,x,y,z) implicit none integer*4 MX,NAT,ii(4,MX),ia,i1,i2,j1,j2,ix,jx real*8 u(3,3,MX),x(MX),v1(3),v2(3),v3(3),y(MX),z(MX) c 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) 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) 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 endif 1 continue return end SUBROUTINE TTTT(MX3,ALPHA,A,G,NAT,ICO,xx,y,z,X,MX) 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(MX3,3,3),G(MX3,3,3),A(MX3,3,3,3),xx(*),X(3,MX), 1y(*),z(*) 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) v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n return end subroutine updt(na1,na2,wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz, 1wcmz,wsux,wcux,wsuy,wcuy,wsuz,wcuz,t,sc,xm,ym,zm,xmm,ymm,zmm, 2bohr,dtau,x,y,z,A,P,nat0,ux,luse) c dipole moments -> ABS,VCD implicit none logical luse(*) integer*4 na1,na2,ia,nw,nat0 real*8 wmin,dw,mx,my,mz,ux,uy,uz,t,sc,bohr,dtau,x(*),y(*),z(*), 2xm(*),ym(*),zm(*),xmm(*),ymm(*),zmm(*),cf1,cf2,vx,vy,vz,ax,ay,az, 3P(nat0,3,3),A(nat0,3,3), 2wsux(*),wsuy(*),wsuz(*),wsmx(*),wsmy(*),wsmz(*), 3wcux(*),wcuy(*),wcuz(*),wcmx(*),wcmy(*),wcmz(*) c mx=0.0d0 my=0.0d0 mz=0.0d0 ux=0.0d0 uy=0.0d0 uz=0.0d0 c constants changing vx,ax into atomic units cf1=1.0d0/bohr/dtau cf2=1.0d0/bohr/dtau/dtau do 201 ia=na1,na2 if(luse(ia))then vx=(x(ia)-xm(ia))*cf1 vy=(y(ia)-ym(ia))*cf1 vz=(z(ia)-zm(ia))*cf1 ax=(x(ia)+xmm(ia)-xm(ia)*2.0d0)*cf2 ay=(y(ia)+ymm(ia)-ym(ia)*2.0d0)*cf2 az=(z(ia)+zmm(ia)-zm(ia)*2.0d0)*cf2 mx=mx+ax*A(ia,1,1)+ay*A(ia,2,1)+az*A(ia,3,1) my=my+ax*A(ia,1,2)+ay*A(ia,2,2)+az*A(ia,3,2) mz=mz+ax*A(ia,1,3)+ay*A(ia,2,3)+az*A(ia,3,3) ux=ux+vx*P(ia,1,1)+vy*P(ia,2,1)+vz*P(ia,3,1) uy=uy+vx*P(ia,1,2)+vy*P(ia,2,2)+vz*P(ia,3,2) uz=uz+vx*P(ia,1,3)+vy*P(ia,2,3)+vz*P(ia,3,3) endif 201 continue call update2(wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz, 1wsux,wcux,wsuy,wcuy,wsuz,wcuz,mx,my,mz,ux,uy,uz,t,sc) return end subroutine rg(n,x,y,z,zim,i12,ic,nat0,maxval) implicit none integer*4 n,ia,idum,maxval,nat0,ic,ii real*8 x(*),y(*),z(*) integer*4 zim(*) integer*4 i12(maxval,nat0) 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) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind,atomic(*) PARAMETER (nat0=300,mxb=nat0,IOX=130) c nat0 maximum number of atoms c mxb maximum number of overlaped atoms for each constant c IOX max # of overlaps common/cctnrest/FF(3*nat0,3*nat0),R(3,nat0), 1IND(nat0),IB(nat0,mxb),ND(nat0), 1NDS(nat0),RS(3,nat0),FFS(3*nat0,3*nat0),IS(nat0,mxb),AM(nat0,3,3), 2INDBIG(IOX,nat0),IBONDSB(8,nat0),IBONDS(8,nat0), 3NBO(nat0),NBOB(nat0),PS(nat0,3,3),AS(nat0,3,3), 4ALPHAS(3*nat0,3,3),AAS(3*nat0,3,3,3),GS(3*nat0,3,3), 5AAT(3*nat0,3,3,3),GT(3*nat0,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,nat0),x(*),y(*),z(*) common/cname/asname(IOX,nat0,3,3), 1psname(IOX,nat0,3,3),aaname(IOX,3*nat0,3,3,3), 1gname(IOX,3*nat0,3,3), 2alphaname(IOX,3*nat0,3,3),rname(IOX,3,nat0),NATNAME(IOX),NAME, 3ffname(IOX,3*nat0,3*nat0) 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(nat0) common/usec/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,nat0,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(nat0,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') CALL READRXS(nat0,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(nat0,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(nat0,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(nat0,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(6,*)' A and G obtained from alpha as DO' CALL TTTT1(3*nat0,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*nat0,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(nat0,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(nat0,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(nat0,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO for small' CALL TTTT1(3*nat0,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*nat0,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(nat0,ALPHAS,AAT,GT,NATS) CLOSE(15) ENDIF write(6,*)'Raman tensors read' endif 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) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind PARAMETER (nat0=300,mxb=nat0,IOX=130) c nat0 maximum number of atoms c mxb maximum number of overlaped atoms for each constant c IOX max # of overlaps common/cctnrest/FF(3*nat0,3*nat0),R(3,nat0),IND(nat0), 1IB(nat0,mxb),ND(nat0), 1NDS(nat0),RS(3,nat0),FFS(3*nat0,3*nat0),IS(nat0,mxb),AM(nat0,3,3), 2INDBIG(IOX,nat0),IBONDSB(8,nat0),IBONDS(8,nat0), 3NBO(nat0),NBOB(nat0),PS(nat0,3,3),AS(nat0,3,3), 4ALPHAS(3*nat0,3,3),AAS(3*nat0,3,3,3),GS(3*nat0,3,3), 5AAT(3*nat0,3,3,3),GT(3*nat0,3,3) real*8 PB(nat0,3,3),AB(nat0,3,3),ALPHA(3*nat0,3,3), 1AA(3*nat0,3,3,3),G(3*nat0,3,3) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM 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, 1AM,INDBIG,NO,IBONDSB,NBOB,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,IERR,LWR) 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(nat0,PB,AB,R,NATi) IF(LROA.or.LRAM)THEN CALL TTTT1(3*nat0,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) implicit none INTEGER*4 nat0,mxb,IOX PARAMETER (nat0=300,mxb=nat0,IOX=130) INTEGER*4 okind,IND(nat0),IS(nat0,mxb),IB(nat0,mxb), 1INDBIG(IOX,nat0), 4IBONDS(8,nat0),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,nat0),nats_o(IOX) real*8 FFB(3*nat0,3*nat0),FFS(3*nat0,3*nat0),RB(3,nat0), 1RS(3,NATS), 1TS(3),TB(3),AM(nat0,3,3),B(3,3),TIJ(3),OTIJ(3),PB(nat0,3,3), 2AB(nat0,3,3),PS(nat0,3,3),AS(nat0,3,3),WF(IOX),ALPHA(3*nat0,3,3), 3ALPHAS(3*nat0,3,3),AA(3*nat0,3,3,3),AAS(3*nat0,3,3,3), 1G(3*nat0,3,3), 4GS(3*nat0,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 LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LNAMES,LHALTONERROR,LWR 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,nat0,3,3), 1pname(IOX,nat0,3,3),aaname(IOX,3*nat0,3,3,3), 1gname(IOX,3*nat0,3,3), 2alphaname(IOX,3*nat0,3,3),rname(IOX,3,nat0),natname(IOX),NAME, 1ffname(IOX,3*nat0,3*nat0) logical luse(nat0) common/usec/luse c 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 20 if(LWR)write(6,6010)IA,JA 6010 format(' Pair ',2i5) 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 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, 2 IERC,LINV,LSTRICT,rname,NATS,IERR) 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 20 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, 2 IERC,LINV,LSTRICT,rname,NATS,IERR) 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 20 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) 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 202 CONTINUE c NF endif 201 CONTINUE c JA endif 20 CONTINUE c IA 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 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) CHARACTER*78 TITDUM integer*4 BT(7),IBONDS(8,N0) READ(2,2000) TITDUM 2000 FORMAT(A78) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > nat0 !') 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) 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 (nat0=300) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(nat0,3),XB(nat0,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=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,*)' 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) 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 NAT,ITER,I,nat0 REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,A,XS,XB,RTOL,RMS,ANG(3),FU PARAMETER (nat0=300) COMMON/MATVAR/A(3,3),XS(nat0,3),XB(nat0,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 (nat0=300) DIMENSION P(N0,3,3),A(N0,3,3),AI(nat0,3,3),AJ(nat0,3,3), 1PV(nat0,3,3),R(3,nat0),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 (nat0=300) DIMENSION P(N0,3,3),A(N0,3,3),R(3,nat0),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 (nat0=300) DIMENSION VCD(nat0,3,3),VEL(nat0,3,3),DIP(nat0,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 (nat0=300) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,nat0) 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,nat0T,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,nat0T,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(nat0T,mxb),IS(nat0T,mxb),IND(nat0T),INDBIG(IOX,nat0T) 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 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,IPO, 1IA,JA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS,IERR) 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 IOX,mxb,nat0 PARAMETER (nat0=300,mxb=nat0) integer*4 IA,IPASS,II,NAT,INDBIG(IOX,nat0),IBA,IBONDS(8,nat0), 1NBOB(*),JB,IIA,NATS,natname(*),IX,IND(*),iu,IERC,I,IC,IERR, 2IPO,J,JA real*8 RMS,RS(3,nat0),RB(3,nat0),TB(3),TS(3),rname(IOX,3,nat0) 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 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) 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,G,AA,PB,AB,N) c project translations/rotations from the tensors IMPLICIT none integer*4 nat0,I,J,K,IA,IX,IIA,IIX,L,N,II PARAMETER (nat0=300) real*8 x(*),y(*),z(*),PB(nat0,3,3),AB(nat0,3,3),ALPHA(3*nat0,3,3), 1AA(3*nat0,3,3,3),G(3*nat0,3,3),A(3*nat0,3*nat0),TEM(3*nat0) c c N is number of atoms here CALL DOMA(x,y,z,A,n) c DO 1 J=1,3 DO 11 I=1,3*N TEM(I)=0.0d0 DO 11 IIA=1,N DO 11 IIX=1,3 11 TEM(I)=TEM(I)+PB(IIA,IIX,J)*A(I,IIX+3*(IIA-1)) do 1 IA=1,N do 1 IX=1,3 1 PB(IA,IX,J)=TEM(IX+3*(IA-1)) c DO 2 J=1,3 DO 22 I=1,3*N TEM(I)=0.0d0 DO 22 IIA=1,N DO 22 IIX=1,3 22 TEM(I)=TEM(I)+AB(IIA,IIX,J)*A(I,IIX+3*(IIA-1)) do 2 IA=1,N do 2 IX=1,3 2 AB(IA,IX,J)=TEM(IX+3*(IA-1)) c DO 3 J=1,3 DO 3 K=1,3 DO 33 I=1,3*N TEM(I)=0.0d0 DO 33 II=1,N 33 TEM(I)=TEM(I)+ALPHA(II,J,K)*A(I,II) do 3 I=1,3*N 3 ALPHA(I,J,K)=TEM(I) RETURN END C ============================================================ SUBROUTINE DOMA(x,y,z,A,nat) implicit none integer*4 nat0,n,nat,N6,ICOL,I,J,II,MX3 PARAMETER (nat0=300,MX3=3*nat0) real*8 x(*),y(*),z(*),A(MX3,MX3),TEM(MX3,MX3),S 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)=-y(I) A(3*I-1,4)= x(I) A(3*I-1,5)=-z(I) A(3*I ,5)= y(I) A(3*I-2,6)= z(I) 1 A(3*I ,6)=-x(I) N6=6 CALL ORT(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.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