PROGRAM NEW4FT2 IMPLICIT none integer*4 MX3,MX,NV0,np0,nt0 PARAMETER (MX3=3*21000,MX=MX3/3,NV0=5000,np0=10000,nt0=1000000) integer*4 j33(MX3),NC(NV0,MX),NA,NM(MX),NOAT,NOSC,I,J,IA,ip, 1K,nn,it,ii,NP,nmd,nrep,nrep0,i1,i2,irow,icol,nb,j33b(MX3) real*8 a33(MX3),ZIM(MX),CHAR(MX),AMULT(MX3),SCFAC(NV0),ZM(MX3), 1CONST,c1,CHATOT,BM,random,a33b(MX3), 2fi,WMIN,WMAX,DELW,dt,dd,d,s0r(np0),s0i(np0),sti(MX3),stem(MX3), 3str(MX3),ctr,ei,tir,wi,tii,t,AXR(MX3),ANORM,sf, 4str0(MX3),sti0(MX3),con,wj,P(MX3,3),A(MX3,3),sbi(3*np0), 5sbr(3*np0),xi,xa,tr(nt0),XMAX,XMIN,ani,anr,tau,time, 6sdi(np0),sdr(np0),ses(np0),ces(np0),bi,br,vecr(MX3) real*8 uxr,uyr,uzr,ce,se,txr,tyr,tzr,txi,tyi,tzi,sum,vecc logical lex,LG,LDW,LEB,LDB,lex2 character*5 s5 CONST=4.359828d0/0.5291772d0/0.5291772d0 c1=1302.828d0 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-3+J)=1/SQRT(ZIM(I)) WRITE(*,778)BM 778 FORMAT(' TOTAL MASS ',F20.6) CLOSE(15) WRITE(*,*)' FTRY.INP READ IN' 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' c c set the options: XMIN=100.0d0 XMAX=4000.0d0 WMIN=100.0d0 WMAX=4000.0d0 DELW=10.0d0 NP=4000 nmd=10000 dt=0.01d0 LG=.false. LDW=.false. LDB=.false. LEB=.false. nrep=0 nrep0=100 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.'XMIN')read(77,*)XMIN if(s5(1:4).eq.'XMAX')read(77,*)XMAX 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.'LDB')read(77,*)LDB if(s5(1:4).eq.'NREP')read(77,*)nrep0 goto 111 112 close(77) endif write(6,6001)WMIN,WMAX,DELW,NP,nmd,nrep0,dt 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) write(6,*)' Double-time integration: ',LDB 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(nmd.gt.nt0)call report('Too many MD steps') if(NP.gt.np0)call report('Too many points') c c read dipole derivatives call READTEN(MX3,P,A,NOAT) c c initialize power spectra c s0r real part c s0i imaginary part call vz(s0r,0.0d0,NP) call vz(s0i,0.0d0,NP) c initialize dipole spectra call vz(sbi,0.0d0,NP*3) call vz(sbr,0.0d0,NP*3) c c Trial vector-real and imaginary part,str, sti, remember initial call vz(sti0,0.0d0,NA) call vz(str0,0.0d0,NA) do 10 i=1,NA sti0(i)=random(0) 10 str0(i)=random(0) call normv(str0,1,NA,ANORM) call normv(sti0,1,NA,ANORM) do 9 i=1,NA str(i)=str0(i) 9 sti(i)=sti0(i) c c make squared matrix: inquire(file='A32.SCR',exist=lex2) if(lex2)then write(6,*)'Squared matrix found on disk' else write(6,*)'Making squared matrix' open(49,file='A3.SCR',form='unformatted') open(50,file='A32.SCR',form='unformatted') ip=0 do 41 irow=1,NA if(irow*100/NA.gt.ip)then ip=irow*100/NA write(6,*)ip,'%' endif rewind 49 DO 43 I=1,irow-1 43 read(49) read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) call vz(vecr,0.0d0,NA) do 44 ii=1,nn vecr(j33(ii))=a33(ii) 44 if(j33(ii).eq.irow)vecr(j33(ii))=vecr(j33(ii))+a33(ii) nb=0 do 42 icol=1,NA rewind 49 DO 431 I=1,icol-1 431 read(49) read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) sum=0.0d0 do 46 ii=1,nn vecc=a33(ii) if(j33(ii).eq.icol)vecc=vecc*2.0d0 46 sum=sum+vecc*vecr(j33(ii)) if(icol.eq.irow)sum=sum/2.0d0 if(sum.gt.1.0d-15)then nb=nb+1 j33b(nb)=icol a33b(nb)=sum endif 42 continue 41 write(50)nb,(j33b(ii),ii=1,nb),(a33b(ii),ii=1,nb) close(49) write(6,*)'Done' close(50) endif c c do the propagation: open(49,file='A32.SCR',form='unformatted') if(LEB)open(39,file='DEBUG.TXT') t=0.0d0 nrep=0 do 3 it=1,nmd t=t+dt c c mutliply st by FF, store in AX: call vz(AXR,0.0d0,NA) 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) AXR(J)=AXR(J)+a33(ii)*str(I) 4 AXR(I)=AXR(I)+a33(ii)*str(J) c c new vector S(t+dt)=2S(t)-S(t-dt)-F^2S(t)dt^2 do 5 i=1,NA stem(i)=str(i) 5 str(i)=2*str(i)-sti(i)-AXR(i)*dt**2 do 51 i=1,NA 51 sti(i)=stem(i) c c normalize call normv(str,1,NA,ANORM) if(dabs(ANORM-1.0d0).gt.0.3d0) 1write(6,*)'norm too high',ANORM c c autocorrelation function: ctr=0.0d0 c dipole moment uxr=0.0d0 uyr=0.0d0 uzr=0.0d0 do 7 i=1,NA uxr=uxr+P(i,1)*str(i) uyr=uyr+P(i,2)*str(i) uzr=uzr+P(i,3)*str(i) 7 ctr=ctr+str(i)*str0(i) tr(it)=ctr if(LEB)write(39,3939)t,ctr,uxr,uyr,uzr 3939 format(9f12.6) c c accumulate FT - double w-integration to disperse it ctr=ctr*dt uxr=uxr*dt uyr=uyr*dt uzr=uzr*dt dd=(WMAX-WMIN)/dble(NP-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 tii=-se*ctr txr=ce*uxr txi=-se*uxr tyr=ce*uyr tyi=-se*uyr tzr=ce*uzr tzi=-se*uzr if(LDW)then wj=WMIN-dd do 81 j=1,np wj=wj+dd con=sf(d,LG,wi,wj) s0i(j)=s0i(j)+con*tir 81 s0r(j)=s0r(j)+con*tii else s0i(i)=s0i(i)+tir s0r(i)=s0r(i)+tii sbi(3*(i-1)+1)=sbi(3*(i-1)+1)+txr sbr(3*(i-1)+1)=sbr(3*(i-1)+1)+txi sbi(3*(i-1)+2)=sbi(3*(i-1)+2)+tyr sbr(3*(i-1)+2)=sbr(3*(i-1)+2)+tyi sbi(3*(i-1)+3)=sbi(3*(i-1)+3)+tzr sbr(3*(i-1)+3)=sbr(3*(i-1)+3)+tzi endif 8 continue c c write the spectrum nrep=nrep+1 if(nrep.eq.nrep0)then nrep=0 open(9,file='S.PRN') open(19,file='SB.PRN') wi=WMIN-dd do 6 i=1,NP wi=wi+dd write(19,494)wi,sbr(3*(i-1)+1)**2+sbi(3*(i-1)+1)**2+ 1 sbr(3*(i-1)+2)**2+sbi(3*(i-1)+2)**2+ 1 sbr(3*(i-1)+3)**2+sbi(3*(i-1)+3)**2 6 write(9,494)wi,s0i(i)**2+s0r(i)**2 494 format(g15.6,' ',g15.6) close(9) close(19) write(6,*)it,' S.PRN and SB.PRN written' if(LDB)then write(6,6009)XMIN,XMAX 6009 format(' Band diagonalization integration',/, 1 ' xmin = ',f10.2,' xmax = ',f10.2) xi=XMIN**2/c1**2 xa=XMAX**2/c1**2 call vz(sdr,0.0d0,NP) call vz(sdi,0.0d0,NP) time=0.0d0 do 61 i1=1,it time=time+dt wi=WMIN-dd do 62 i=1,NP wi=wi+dd ei=(wi/c1)**2*time ces(i)=cos(ei)*dt**2 62 ses(i)=sin(ei)*dt**2 tau=dt*i1 do 61 i2=-i1+1,it-i1 tau=tau-dt if(i2.eq.0)then anr=xa-xi ani=0.0d0 else anr=(-sin(xi*tau)+sin(xa*tau))/tau ani=( cos(xi*tau)-cos(xa*tau))/tau endif br=anr*tr(i1+i2) bi=ani*tr(i1+i2) do 61 i=1,NP sdr(i)=sdr(i)+br*ces(i)+bi*ses(i) 61 sdi(i)=sdi(i)+bi*ces(i)-br*ses(i) open(9,file='SD.PRN') wi=WMIN-dd do 63 i=1,NP wi=wi+dd 63 write(9,494)wi,sdr(i)**2+sdi(i)**2 close(9) write(6,*)it,' SD.PRN written' endif endif 3 continue close(49) if(LEB)close(39) 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 normv(A,N1,N2,ANORM) real*8 A(*),ANORM,f integer*4 N1,N2,I ANORM=0.0d0 DO 1 I=N1,N2 1 ANORM=ANORM+A(I)**2 f=1.0d0/DSQRT(ANORM) DO 2 I=N1,N2 2 A(I)=A(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