PROGRAM NEW4FT IMPLICIT none integer*4 MX3,MX,NV0,np0,ns0,nsim,it0,mclock PARAMETER (MX3=3*21000,MX=MX3/3,NV0=5000,np0=4000,ns0=100) integer*4 j33(MX3),NC(NV0,MX),NA,NM(MX),NOAT,NOSC,I,J,IA, 1K,nn,it,ii,NP,nmd,nrep,nrep0,nhalf,ix,i1,i2,i3,is,nff,i33 real*8 sti(MX3*ns0),str(MX3*ns0),str0(MX3*ns0),sti0(MX3*ns0), 1AXR(MX3*ns0),AXI(MX3*ns0),st2r(MX3*ns0),st2i(MX3*ns0), 1st3r(MX3*ns0),st3i(MX3*ns0), 1smmr(MX3*ns0),smmi(MX3*ns0),smr(MX3*ns0),smi(MX3*ns0), 1sm3r(MX3*ns0),sm3i(MX3*ns0), 2a33(MX3),ZM(MX3),AMULT(MX3),P(MX3,3),A(MX3,3),ZIM(MX),CHAR(MX), 3sbi(3*np0*ns0),sbr(3*np0*ns0),s0r(np0*ns0),s0i(np0*ns0), 3sai(3*np0*ns0),sar(3*np0*ns0),sff(np0), 3SCFAC(NV0),tr(3),ti(3),ur(3),ui(3),mr(3),mi(3),qr(3),qi(3), 1s1,s2,s3,s4,cti,ctr,ei,tir,wi,s55, 4tii,t,ANORM,sf,CONST,c1,CHATOT,BM,alpha,al0,el0,random,ce,se, 5fi,WMIN,WMAX,DELW,dt,dd,d,con,wj,FC,FD,c2,pi logical lex,LG,LDW,LEB,LVEL,LWS,LSEC,LT3 character*5 s5 CONST=4.359828d0/0.5291772d0/0.5291772d0 c1=1302.828d0 pi=4.0d0*atan(1.0d0) it0=mclock() C WRITE(*,60000) 60000 FORMAT(' SPECTRA GENERATION BY TIME-PROPAGATION OF',/, 1 ' CONTINUEOUS FOURIER TRANSFORMATION',/, 1 ' PC version - in cartesian coordinates',/, 2 ' Scaling of masses!',/,/, 2 ' For giant force fields',/,/, 3 ' Petr Bour, Prague 2008',/,/) C OPEN(17,FILE='FILE.X') READ(17,*) READ(17,*)NOAT CLOSE(17) NA=3*NOAT if(NA.gt.MX3)call report('too many atoms') C C Read-in scalefactors: OPEN(15,FILE='FTRY.INP') READ(15,*) READ(15,*)NOSC if(NOSC.gt.NV0)call report('too many scale factors') DO 222 I=1,NOSC 222 READ(15,*) SCFAC(I),NM(I),(NC(I,J),J=1,NM(I)) 1040 FORMAT(6G12.6) WRITE(*,*)' Scale factors read in...' C C Read-in atomic charges for FPC: READ(15,*) (CHAR(I),I=1,NOAT) CHATOT=0.0d0 DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHAR(I) WRITE(*,777)CHATOT,NOAT 777 FORMAT(' TOTAL CHARGE',F20.6,/I6,' atoms') C C Read-in atomic masses: READ(15,*) READ(15,1040)(ZIM(I),I=1,NOAT) BM=0.0d0 DO 102 IA=1,NOAT 102 BM=BM+ZIM(IA) DO 220 I=1,NOAT DO 220 J=1,3 220 ZM(3*(I-1)+J)=1/SQRT(ZIM(I)) WRITE(*,778)BM 778 FORMAT(' TOTAL MASS ',F20.6) CLOSE(15) WRITE(*,*)' FTRY.INP READ IN' write(6,6009)mclock()-it0 6009 format(' Time: ',i10,' msec') C c scale if required: DO 201 I=1,NA 201 AMULT(I)=1.0d0 IF(NOSC.GT.0)THEN DO 202 I=1,NOSC t=1.0d0/SQRT(SCFAC(I)) DO 202 J=1,NM(I) DO 202 K=1,3 202 AMULT(3*(NC(I,J)-1)+K)=t ENDIF c WRITE(*,*)' Making scaled mass-weighed FF from AFILE.FC' OPEN(20,FILE='AFILE.FC',STATUS='OLD',form='unformatted') OPEN(21,FILE='A3.SCR',form='unformatted') do 20 i=1,NA fi=ZM(i)*AMULT(i)*CONST read(20)nn,(j33(j),j=1,nn),(a33(j),j=1,nn) write(21)nn,(j33(j),j=1,nn),(a33(j)*fi*ZM(j33(j)) 1*AMULT(j33(j)) ,j=1,nn) 20 continue close(20) close(21) close(22) WRITE(*,*)' Cartesian FF rewritten into A3.SCR' write(6,6009)mclock()-it0 c c set the options: WMIN=100.0d0 WMAX=4000.0d0 DELW=10.0d0 NP=4000 nmd=10000 dt=0.01d0 LG=.true. LDW=.true. LEB=.false. LVEL=.false. LSEC=.true. LT3=.false. LWS=.false. alpha=0.0d0 al0=1500.0d0 nrep=0 nrep0=1000 nhalf=1 nsim=10 inquire(file='NEW4FT.OPT',exist=lex) if(lex)then open(77,file='NEW4FT.OPT') 111 read(77,770,end=112,err=112)s5 770 format(a5) if(s5(1:4).eq.'WMIN')read(77,*)WMIN if(s5(1:4).eq.'WMAX')read(77,*)WMAX if(s5(1:4).eq.'DELW')read(77,*)DELW if(s5(1:3).eq.'NMD')read(77,*)NMD if(s5(1:2).eq.'DT')read(77,*)dt if(s5(1:2).eq.'NP')read(77,*)NP if(s5(1:2).eq.'LG')read(77,*)LG if(s5(1:3).eq.'LDW')read(77,*)LDW if(s5(1:3).eq.'LEB')read(77,*)LEB if(s5(1:3).eq.'LVE')read(77,*)LVEL if(s5(1:3).eq.'LSE')read(77,*)LSEC if(s5(1:3).eq.'LT3')read(77,*)LT3 if(s5(1:3).eq.'LWS')read(77,*)LWS if(s5(1:3).eq.'AL0')read(77,*)al0 if(s5(1:4).eq.'NREP')read(77,*)nrep0 if(s5(1:4).eq.'ALPH')read(77,*)alpha if(s5(1:4).eq.'NHAL')read(77,*)nhalf if(s5(1:4).eq.'NSIM')read(77,*)nsim goto 111 112 close(77) endif if(nhalf.eq.0)nhalf=1 write(6,6001)WMIN,WMAX,DELW,NP,nmd,nrep0,dt,alpha,al0,nsim,nhalf 6001 format(' Spectral interval ',f10.1,' - ',f10.1,' cm-1',/, 1 ' bandwidth ',f10.2,' cm-1',/, 2 ' number of spectral points ',i15,/, 2 ' number of MD points ',i15,/, 2 ' write spectra each step ',i15,/, 2 ' integration time step ',f15.3,/, 2 ' damping coeffs ',f15.3,/, 2 ' damping limit ',f15.3,/, 2 ' number of simultanous spec',i15,/, 2 ' pure propagation MD steps ',i15) write(6,*)' Velocity integration : ',LVEL write(6,*)' Second derivs on fly : ',LSEC write(6,*)' Third derivs on fly : ',LT3 write(6,*)' Write all spectra : ',LWS write(6,*)' Double-ener integration: ',LDW write(6,*)' Debug DEBUG.TXT : ',LEB if(LG)then write(6,*)'Gaussian peaks' d=DELW/2.0d0/dsqrt(log(2.0d0)) else write(6,*)'Lorentzian peaks' d=DELW/2.0d0 endif if(NP.gt.np0)call report('Too many points') if(nsim.gt.ns0)call report('Too many simultaneous spectra') el0=al0**2/c1**2 c c pre-calculate convolution points dd=(WMAX-WMIN)/dble(NP-1) if(LDW)then nff=0 wj=-dd do 11 j=1,NP wj=wj+dd con=sf(d,LG,0.0d0,wj) if(con.gt.1.0d-3)then nff=nff+1 sff(nff)=con else goto 12 endif 11 continue 12 continue write(6,*)nff,' convolution points' endif c c read dipole derivatives call READTEN(MX3,P,A,NOAT) c c Make mass-weighed derivatives for the propagation FD=dsqrt(388.9219d0) FC=131.14d-8 do 10 i=1,NA do 10 ix=1,3 P(i,ix)=P(i,ix)*ZM(i)*FD 10 A(i,ix)=A(i,ix)*ZM(i)*FC write(6,6009)mclock()-it0 c c initialize power spectra s0r real part s0i imaginary part call vz(s0r,0.0d0,NP*nsim) call vz(s0i,0.0d0,NP*nsim) c initialize dipole spectra call vz(sbi,0.0d0,NP*3*nsim) call vz(sbr,0.0d0,NP*3*nsim) call vz(sai,0.0d0,NP*3*nsim) call vz(sar,0.0d0,NP*3*nsim) c c Initialize trial vector-real and imaginary part,str, sti call vz(sti,0.0d0,NA*nsim) do 80 i=1,NA*nsim 80 str(i)=random(0)-0.5d0 do 81 i=1,nsim 81 call normvc(str,sti,NA*(i-1)+1,NA*i,ANORM) c c Initialize second,third derivatives and old vectors call vz(st3i,0.0d0,NA*nsim) call vz(st3r,0.0d0,NA*nsim) call vz(st2i,0.0d0,NA*nsim) call vz(st2r,0.0d0,NA*nsim) do 83 i=1,NA*nsim sm3r(i)=str(i) sm3i(i)=sti(i) smmr(i)=str(i) smmi(i)=sti(i) smr(i)=str(i) 83 smi(i)=sti(i) c c do the propagation: open(49,file='A3.SCR',form='unformatted') if(LEB)open(89,file='DEBUG.TXT') t=0.0d0 nrep=0 do 3 it=1,nmd nrep=nrep+1 t=t+dt c c mutliply st by FF, store in AX: call vz(AXR,0.0d0,NA*nsim) call vz(AXI,0.0d0,NA*nsim) rewind 49 DO 4 I=1,NA read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 4 ii=1,nn J=j33(ii) do 4 is=1,nsim i1=NA*(is-1) AXR(i1+J)=AXR(i1+J)+a33(ii)*str(i1+I) AXR(i1+I)=AXR(i1+I)+a33(ii)*str(i1+J) AXI(i1+J)=AXI(i1+J)+a33(ii)*sti(i1+I) 4 AXI(i1+I)=AXI(i1+I)+a33(ii)*sti(i1+J) c c new vectors do 5 i=1,NA*nsim str(i)=str(i)-AXI(i)*dt 5 sti(i)=sti(i)+AXR(i)*dt if(alpha.gt.1.0d-9)then do 53 i=1,NA*nsim str(i)=str(i)-alpha*(AXR(i)-el0*str(i))*dt 53 sti(i)=sti(i)-alpha*(AXI(i)-el0*sti(i))*dt endif if(LSEC)then do 52 i=1,NA*nsim str(i)=str(i)+st2r(i) 52 sti(i)=sti(i)+st2i(i) endif if(LT3)then do 54 i=1,NA*nsim str(i)=str(i)+st3r(i) 54 sti(i)=sti(i)+st3i(i) endif c c normalize do 51 i=1,nsim call normvc(str,sti,NA*(i-1)+1,NA*i,ANORM) 51 if(dabs(ANORM-1.0d0).gt.0.3d0)write(6,*)i,'norm too high',ANORM c c calculate second derivatives*dt^2*0.5 c calculate third derivatives*dt^3/6 do 84 i=1,NA*nsim st3r(i)=(3.0d0*(smmr(i)-smr(i))+str(i)-sm3r(i))/6.0d0 st3i(i)=(3.0d0*(smmi(i)-smi(i))+sti(i)-sm3i(i))/6.0d0 st2r(i)=(str(i)+smmr(i))*0.5d0-smr(i) 84 st2i(i)=(sti(i)+smmi(i))*0.5d0-smi(i) c c remember this and previous vectors do 85 i=1,NA*nsim sm3r(i)=smmr(i) sm3i(i)=smmi(i) smmr(i)=smr(i) smmi(i)=smi(i) smr(i)=str(i) 85 smi(i)=sti(i) c if(it.eq.nhalf)then write(6,*)it,' step, time-integration started' c write(6,*)'damping switched off, reference vector stored' c alpha=0.0d0 do 9 i=1,NA*nsim str0(i)=str(i) 9 sti0(i)=sti(i) endif if(it.gt.nhalf)then c do 711 is=1,nsim i1=NA*(is-1) c time-dependent dipole moments and autocorrelation function: call vz(ur,0.0d0,3) call vz(ui,0.0d0,3) call vz(mr,0.0d0,3) call vz(mi,0.0d0,3) cti=0.0d0 ctr=0.0d0 do 71 i=1,NA if(LVEL)then ctr=ctr+AXR(i1+i) cti=cti+AXI(i1+i) else ctr=ctr+str(i1+i)*str0(i1+i)-sti(i1+i)*sti0(i1+i) cti=cti+str(i1+i)*sti0(i1+i)+sti(i1+i)*str0(i1+i) endif do 71 ix=1,3 mr(ix)=mr(ix)+A(i,ix)*str(i1+i) mi(ix)=mi(ix)+A(i,ix)*sti(i1+i) ur(ix)=ur(ix)+P(i,ix)*str(i1+i) 71 ui(ix)=ui(ix)+P(i,ix)*sti(i1+i) if(LEB.and.is.eq.1)write(89,3939)t,ctr,cti,ur,ui 3939 format(9f12.6) c c accumulate FT - double w-integration to disperse it ctr=ctr*dt cti=cti*dt do 73 ix=1,3 mr(ix)=mr(ix)*dt mi(ix)=mi(ix)*dt ur(ix)=ur(ix)*dt 73 ui(ix)=ui(ix)*dt i2=np*(is-1) i3=3*np*(is-1) wi=WMIN-dd do 8 i=1,np wi=wi+dd ei=wi**2/c1**2 ce=cos(ei*t) se=sin(ei*t) tir=ce*ctr+se*cti tii=ce*cti-se*ctr do 74 ix=1,3 qr(ix)=ce*mr(ix)+se*mi(ix) qi(ix)=ce*mi(ix)-se*mr(ix) tr(ix)=ce*ur(ix)+se*ui(ix) 74 ti(ix)=ce*ui(ix)-se*ur(ix) if(LDW)then do 82 j=i,min(np,i+nff-1) con=sff(j-i+1) i33=i3+3*(j-1) 82 call sadd(s0r,s0i,con,tir,tii,qr,qi,tr,ti,sar,sai,sbr,sbi, 1 i2+j,i33) do 76 j=max(1,i-nff+1),i-1 con=sff(i-j+1) i33=i3+3*(j-1) 76 call sadd(s0r,s0i,con,tir,tii,qr,qi,tr,ti,sar,sai,sbr,sbi, 1 i2+j,i33) else s0r(i2+i)=s0r(i2+i)+tir s0i(i2+i)=s0i(i2+i)+tii do 75 ix=1,3 sar(i3+3*(i-1)+ix)=sar(i3+3*(i-1)+ix)+qr(ix) sai(i3+3*(i-1)+ix)=sai(i3+3*(i-1)+ix)+qi(ix) sbr(i3+3*(i-1)+ix)=sbr(i3+3*(i-1)+ix)+tr(ix) 75 sbi(i3+3*(i-1)+ix)=sbi(i3+3*(i-1)+ix)+ti(ix) endif 8 continue c 711 continue c c if(LDW)then c2=1.0d0/dble(NA)*2.0d0/pi*2.0d0*d/c1**2/dble(nsim) else c2=1.0d0/dble(NA)*2.0d0/pi/t/dble(nsim) endif c write the spectrum if(nrep.ge.nrep0)then nrep=0 if(LWS)then do 712 is=1,nsim i2=np*(is-1) i3=3*np*(is-1) write(s5,1000)is 1000 format(i5) open(9,file=s5//'S.PRN') open(19,file=s5//'SB.PRN') open(39,file=s5//'SBR.PRN') open(29,file=s5//'SC.PRN') open(59,file=s5//'SCR.PRN') wi=WMIN-dd do 6 i=1,NP wi=wi+dd s1=0.0d0 s2=0.0d0 s3=0.0d0 s4=0.0d0 s55=0.0d0 do 63 ix=1,3 s1=sbr(i3+3*(i-1)+ix)*sbr(i3+3*(i-1)+ix) 63 s4=sbr(i3+3*(i-1)+ix)*sar(i3+3*(i-1)+ix) s2=s0r(i2+i)**2 if(LDW)then s1=s1*108.7d0*wi**2*c2 s4=s4*435.0d0*wi**3*c2 else s1=s1*108.7d0*wi*c2 s4=s4*435.0d0*wi**2*c2 endif if(s2.gt.1.0d-10)s3=s1/dsqrt(s2) if(s2.gt.1.0d-10)s55=s4/dsqrt(s2) write(19,494)wi,s1 write(39,494)wi,s4 write(59,494)wi,s55 write(29,494)wi,s3 6 write(9,494)wi,s2 close(9) close(19) close(29) close(59) 712 close(39) endif open(9,file='S.PRN') open(19,file='SB.PRN') open(39,file='SBR.PRN') open(29,file='SC.PRN') open(59,file='SCR.PRN') wi=WMIN-dd do 64 i=1,NP wi=wi+dd s1=0.0d0 s2=0.0d0 s3=0.0d0 s4=0.0d0 s55=0.0d0 do 61 is=1,nsim i2=np*(is-1)+i s2=s2+s0r(i2)**2 i3=3*np*(is-1)+3*(i-1) do 61 ix=1,3 s1=s1+sbr(i3+ix)*sbr(i3+ix) 61 s4=s4+sbr(i3+ix)*sar(i3+ix) if(LDW)then s1=s1*108.7d0*wi**2*c2 s4=s4*435.0d0*wi**3*c2 else s1=s1*108.7d0*wi*c2 s4=s4*435.0d0*wi**2*c2 endif if(s2.gt.1.0d-10)s3=s1/dsqrt(s2) if(s2.gt.1.0d-10)s55=s4/dsqrt(s2) write(19,494)wi,s1 write( 9,494)wi,s2 write(29,494)wi,s3 write(39,494)wi,s4 64 write(59,494)wi,s55 494 format(g15.6,' ',g15.6) close(9) close(19) close(29) close(39) close(59) write(6,*)it,' Spectra written' write(6,6009)mclock()-it0 endif endif 3 continue close(49) if(LEB)close(89) c stop END c SUBROUTINE VZ(V,Z,N) REAL*8 V(*),Z INTEGER*4 I,N DO 1 I=1,N 1 V(I)=Z RETURN END subroutine report(s) character*(*) s write(6,*)s stop end subroutine normvc(A,B,N1,N2,ANORM) real*8 A(*),ANORM,B(*),f integer*4 N1,N2,I ANORM=0.0d0 DO 1 I=N1,N2 1 ANORM=ANORM+A(I)**2+B(I)**2 f=1.0d0/DSQRT(ANORM) DO 2 I=N1,N2 A(I)=A(I)*f 2 B(I)=B(I)*f return end function sf(d,g,x,x0) implicit none logical g real*8 pi,spi,sf,d,x,x0,dd pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1000.0d0)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif return end SUBROUTINE READTEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3),A(N0,3) open(15,FILE='FILE.TEN') READ (15,*) NAT WRITE(*,*)NAT,' atoms in FILE.TEN' DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(3*(L-1)+J,I),I=1,3) WRITE(6,*)' Electric dipole derivatives read in ' DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(3*(L-1)+J,I),I=1,3) WRITE(6,*)' MAgnetic dipole derivatives read in ' 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,*) close(15) RETURN END subroutine sadd(s0r,s0i,con,tir,tii,qr,qi,tr,ti,sar,sai,sbr,sbi, 1i2,i33) implicit none real*8 s0r(*),s0i(*),sar(*),sai(*),sbr(*),sbi(*),con,qr(*), 1qi(*),tr(*),ti(*),tii,tir integer*4 i2,i33 s0r(i2)=s0r(i2)+con*tir s0i(i2)=s0i(i2)+con*tii sar(i33+1)=sar(i33+1)+con*qr(1) sai(i33+1)=sai(i33+1)+con*qi(1) sbr(i33+1)=sbr(i33+1)+con*tr(1) sbi(i33+1)=sbi(i33+1)+con*ti(1) sar(i33+2)=sar(i33+2)+con*qr(2) sai(i33+2)=sai(i33+2)+con*qi(2) sbr(i33+2)=sbr(i33+2)+con*tr(2) sbi(i33+2)=sbi(i33+2)+con*ti(2) sar(i33+3)=sar(i33+3)+con*qr(3) sai(i33+3)=sai(i33+3)+con*qi(3) sbr(i33+3)=sbr(i33+3)+con*tr(3) sbi(i33+3)=sbi(i33+3)+con*ti(3) return end