PROGRAM NEW4FTR c with Raman extension IMPLICIT none integer*4 MX3,MX,NV0,np0,nsim,it0,mclock,nb0,nproc0 PARAMETER (MX=21000,MX3=MX*3,NV0=5,np0=4000, 1nb0=8000000,nproc0=99) c NV0 - number of scale factors c MX - number of atoms c np0 - number of spectral points c c big buffer - allocate dynamically later integer*4 jbb(nb0) real*8 abb(nb0) c c local variables: integer*4 j33(MX3),mm(MX3),nn(MX3),i9,i7,io,iti,itj,itkl,l, 1NC(NV0,MX),NA,NM(MX),NOAT,NOSC,I,J,IA,iel, 1K,it,ii,NP,nmd,nhalf,ix,is,nff,iok,inn,jn, 1nproc,omp_get_num_procs,omp_get_max_threads,OMP_GET_THREAD_NUM, 1nstart(nproc0),nend(nproc0),nd,nrest,ipr,it00,seed real*8 sti(MX3),str(MX3),str0(MX3),a33(MX3),AXR(MX3),AXI(MX3), 1st2r(MX3),st2i(MX3),st3r(MX3),st3i(MX3),smmr(MX3),smmi(MX3), 2smr(MX3),smi(MX3),sm3r(MX3),sm3i(MX3),ZM(MX3),AMULT(MX3),P(MX3,3), 3A(MX3,3),ZIM(MX),CHAR(MX),sbi(3*np0),sbr(3*np0),s0r(np0),s0i(np0), 4st(5*np0),sai(3*np0),sar(3*np0),sff(np0), 5SCFAC(NV0),tr(3),ti(3),ur(3),rtr,ann, 6ui(3),mr(3),mi(3),qr(3),qi(3),s1,s2,s3,cti,ctr,ei,tir,wi, 4tii,t,ANORM,sf,CONST,c1,CHATOT,BM,random,ce,se, 5fi,WMIN,WMAX,DELW,dt,dd,d,con,wj,FC,FD,c2,pi,ALPHA(MX3,9), 6G(MX3,9),AT(MX3,27),talr(9),tali(9),alr(9),ali(9),tgr(9),tgi(9), 7gi(9),gr(9),tar(27),tai(27),ai(27),ar(27), 8salr(9*np0),sali(9*np0),sgr(9*np0),sgi(9*np0),satr(27*np0), 9sati(27*np0),SAL0,SAL1,SAG0,SAG1,SA1,EPS,stemp,s4,sp logical lex,LG,LDW,LEB,LVEL,LSEC,LT3,lrtr character*5 s5 CONST=4.359828d0/0.5291772d0/0.5291772d0 c1=1302.828d0 pi=4.0d0*atan(1.0d0) it0=mclock() it00=-1 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)ii,(j33(j),j=1,ii),(a33(j),j=1,ii) write(21)ii,(j33(j),j=1,ii),(a33(j)*fi*ZM(j33(j)) 1*AMULT(j33(j)) ,j=1,ii) 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=.false. LEB=.false. LVEL=.false. LSEC=.true. LT3=.false. lrtr=.false. nhalf=1 nsim=10 nproc=1 seed=0 stemp=298.0d0 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:4).eq.'NHAL')read(77,*)nhalf if(s5(1:4).eq.'NSIM')read(77,*)nsim if(s5(1:4).eq.'#PRO')read(77,*)nproc if(s5(1:4).eq.'RAND')read(77,*)seed if(s5(1:4).eq.'LRTR')read(77,*)lrtr if(s5(1:4).eq.'TEMP')read(77,*)stemp goto 111 112 close(77) endif if(nhalf.eq.0)nhalf=1 write(6,6001)WMIN,WMAX,DELW,NP,nmd,dt,nsim,nhalf, 1omp_get_num_procs(),omp_get_max_threads(),nproc,seed 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 ' integration time step ',f15.3,/, 2 ' number of simultanous spec',i15,/, 2 ' pure propagation MD steps ',i15,/, 2 ' available processors ',i15,/, 2 ' maximum of threads ',i15,/, 2 ' requested processors ',i15,/, 2 ' random seed ',i15,/) if(nproc.gt.omp_get_num_procs()) 1call report('not enough processors') if(nproc.gt.omp_get_max_threads())then write(6,*)'Number of threads will increase' call omp_set_num_threads(nproc) endif call omp_set_dynamic(0) call omp_set_nested(0) c c precalculate parallel distribution: if(nproc.gt.nproc0)call report('too many processors requested') nd=nsim/nproc nrest=nsim-nd*nproc nstart(1)=1 nend(1)=nd if(nrest.ge.1)nend(1)=nend(1)+1 do 7 ipr=2,nproc nstart(ipr)=nend(ipr-1)+1 nend(ipr)=nstart(ipr)+nd-1 7 if(nrest.ge.ipr)nend(ipr)=nend(ipr)+1 write(6,6007) 6007 format(' Parallel distribution:') do 72 ipr=1,nproc 72 write(6,60071)ipr,nstart(ipr),nend(ipr),nend(ipr)-nstart(ipr)+1 60071 format(i2,i5,' - ',i5,i10) write(6,*)' Velocity integration : ',LVEL write(6,*)' Second derivs on fly : ',LSEC write(6,*)' Third derivs on fly : ',LT3 write(6,*)' Double-ener integration: ',LDW write(6,*)' Debug DEBUG.TXT : ',LEB write(6,*)' SQRT(Random) : ',lrtr 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') 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 read POLARIZABILITY derivatives call READTTT(MX3,ALPHA,AT,G,NOAT) c c Make mass-weighed derivatives for the propagation FD=dsqrt(388.9219d0) FC=131.14d-8 do 10 i=1,NA do 101 ix=1,3 P(i,ix)=P(i,ix)*ZM(i)*FD 101 A(i,ix)=A(i,ix)*ZM(i)*FC do 103 ix=1,9 ALPHA(i,ix)=ALPHA(i,ix)*ZM(i) 103 G(i,ix)=G(i,ix)*ZM(i) do 10 ix=1,27 10 AT(i,ix)=AT(i,ix)*ZM(i) write(6,6009)mclock()-it0 c c load the matrix into memory - non-zero elements only c where the line starts : mm c how many elements it has: nn iok=0 iel=0 open(49,file='A3.SCR',form='unformatted') do 441 i=1,NA read(49)nn(i),(j33(ii),ii=1,nn(i)),(a33(ii),ii=1,nn(i)) mm(i)=iel do 441 ii=1,nn(i) iel=iel+1 if(iel.gt.nb0)iok=1 if(iok.eq.0)jbb(iel)=j33(ii) 441 if(iok.eq.0)abb(iel)=a33(ii) close(49) if(iok.gt.0)then write(6,6055)nb0,iel 6055 format(' Current: ',i10,' Needed: ',i10) call report('memory exceeded') endif write(6,4900)(iel*8+iel*4+NA*8)/1000 4900 format(' Matrix loaded into memory',i12,' kb') c c initialize spectra do 2994 i=1,5*NP 2994 st(i)=0.0d0 c c initialize random function do 2993 i=1,13*seed sti(i)=random(0) 2993 if(sti(i).gt.100.0d0)write(6,*)'blbost' c C$OMP PARALLEL DEFAULT(SHARED) C$OMP+ PRIVATE(s0r,s0i,sbi,sbr,sai,sar,sti,str,st3i,st3r,st2i, C$OMP+ i,st2r,sm3r,sm3i,smmr,smmi,t,it,AXR,AXI,ii,J,is, C$OMP+ con,ur,ui,mr,mi,cti,ctr,wi,ei,ce,se,tir,tii,ix,ANORM,ipr, C$OMP+ smr,smi,str0,qr,qi,tr,ti,ix,sali,salr,sgi,sgr,sati,satr, C$OMP+ talr,tali,tgr,tgi,tar,tai,ali,alr,gi,gr,ai,ar,SAL0,SAL1, C$OMP+ SAG0,SAG1,K,L,SA1,itkl,i9,i7,iti,itj,io,inn,jn,ann) ipr=OMP_GET_THREAD_NUM()+1 write(6,*)' Integration started, thread ',ipr do 899 is=nstart(ipr),nend(ipr) if(ipr.eq.1.and.is.eq.1)it00=mclock() if(ipr.eq.1.and.is.eq.2)then it00=mclock()-it00 write(6,6005)it00,it00*(nend(ipr)-nstart(ipr)+1) 6005 format('One spectrum time:',i10,' total estimate:',i11) endif do 8991 i=1,NP s0r(i)=0.0d0 8991 s0i(i)=0.0d0 do 89911 i=1,3*NP sbi(i)=0.0d0 sbr(i)=0.0d0 sai(i)=0.0d0 89911 sar(i)=0.0d0 do 89912 i=1,9*NP sali(i)=0.0d0 salr(i)=0.0d0 sgi(i)=0.0d0 89912 sgr(i)=0.0d0 do 89913 i=1,27*NP sati(i)=0.0d0 89913 satr(i)=0.0d0 c c Initialize trial vector-real and imaginary part,str, sti do 2992 i=1,(ipr-1)*(131+13*is) sti(i)=random(0) 2992 if(sti(i).gt.100.0d0)write(6,*)'blbost' ANORM=0.0d0 do 8992 i=1,NA sti(i)=0.0d0 rtr=random(0)-0.5d0 if(lrtr)then str(i)=dsqrt(dabs(rtr)) if(rtr.lt.0.0d0)str(i)=-str(i) else str(i)=rtr endif 8992 ANORM=ANORM+str(i)**2 do 81 i=1,NA str(i)=str(i)/dsqrt(ANORM) 81 str0(i)=str(i) c c Initialize second,third derivatives and old vectors do 8993 i=1,NA st3r(i)=0.0d0 st3i(i)=0.0d0 st2r(i)=0.0d0 st2i(i)=0.0d0 sm3r(i)=str(i) sm3i(i)=sti(i) smmr(i)=str(i) smmi(i)=sti(i) smr(i)=str(i) 8993 smi(i)=sti(i) c c do the propagation: t=0.0d0 do 3 it=1,nmd t=t+dt c c mutliply st by FF, store in AX: do 4993 i=1,NA AXR(i)=0.0d0 4993 AXI(i)=0.0d0 DO 4 I=1,NA DO 4 ii=1,nn(I) inn=mm(I)+ii jn=jbb(inn) ann=abb(inn) AXR(jn)=AXR(jn)+ann*str(I) AXR(I )=AXR(I )+ann*str(jn) AXI(jn)=AXI(jn)+ann*sti(I) 4 AXI(I )=AXI(I )+ann*sti(jn) do 5 i=1,NA str(i)=str(i)-AXI(i)*dt 5 sti(i)=sti(i)+AXR(i)*dt if(LSEC)then do 52 i=1,NA str(i)=str(i)+st2r(i) 52 sti(i)=sti(i)+st2i(i) endif if(LT3)then do 54 i=1,NA str(i)=str(i)+st3r(i) 54 sti(i)=sti(i)+st3i(i) endif c c normalize ANORM=0.0d0 do 4995 i=1,NA 4995 ANORM=ANORM+str(i)**2+sti(i)**2 ann=1.0d0/dsqrt(ANORM) do 8995 i=1,NA sti(i)=sti(i)*ann 8995 str(i)=str(i)*ann 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 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 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 time-dependent dipole moments and autocorrelation function: do 71 ix=1,3 ur(ix)=0.0d0 ui(ix)=0.0d0 mi(ix)=0.0d0 mr(ix)=0.0d0 do 73 i=1,NA mr(ix)=mr(ix)+A(i,ix)*str(i) mi(ix)=mi(ix)+A(i,ix)*sti(i) ur(ix)=ur(ix)+P(i,ix)*str(i) 73 ui(ix)=ui(ix)+P(i,ix)*sti(i) mr(ix)=mr(ix)*dt mi(ix)=mi(ix)*dt ur(ix)=ur(ix)*dt 71 ui(ix)=ui(ix)*dt do 7111 ix=1,9 alr(ix)=0.0d0 ali(ix)=0.0d0 gr(ix)=0.0d0 gi(ix)=0.0d0 ar(ix)=0.0d0 ai(ix)=0.0d0 ar(ix+9)=0.0d0 ai(ix+9)=0.0d0 ar(ix+18)=0.0d0 ai(ix+18)=0.0d0 do 7311 i=1,NA alr(ix)=alr(ix)+ALPHA(i,ix)*str(i) ali(ix)=ali(ix)+ALPHA(i,ix)*sti(i) gr(ix)=gr(ix)+G(i,ix)*str(i) gi(ix)=gi(ix)+G(i,ix)*sti(i) ar(ix)=ar(ix)+AT(i,ix)*str(i) ai(ix)=ai(ix)+AT(i,ix)*sti(i) ar(ix+9)=ar(ix+9)+AT(i,ix+9)*str(i) ai(ix+9)=ai(ix+9)+AT(i,ix+9)*sti(i) ar(ix+18)=ar(ix+18)+AT(i,ix+18)*str(i) 7311 ai(ix+18)=ai(ix+18)+AT(i,ix+18)*sti(i) alr(ix)=alr(ix)*dt ali(ix)=ali(ix)*dt gr(ix)=gr(ix)*dt gi(ix)=gi(ix)*dt ar(ix)=ar(ix)*dt ai(ix)=ai(ix)*dt ar(ix+9)=ar(ix+9)*dt ai(ix+9)=ai(ix+9)*dt ar(ix+18)=ar(ix+18)*dt 7111 ai(ix+18)=ai(ix+18)*dt cti=0.0d0 ctr=0.0d0 if(LVEL)then do 711 i=1,NA ctr=ctr+AXR(i) 711 cti=cti+AXI(i) else do 712 i=1,NA ctr=ctr+str(i)*str0(i) 712 cti=cti +sti(i)*str0(i) endif ctr=ctr*dt cti=cti*dt c accumulate FT - double w-integration to disperse it: 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) do 7411 ix=1,9 talr(ix)=ce*alr(ix)+se*ali(ix) tali(ix)=ce*ali(ix)-se*alr(ix) tgr(ix)=ce*gr(ix)+se*gi(ix) 7411 tgi(ix)=ce*gi(ix)-se*gr(ix) do 7422 ix=1,27 tar(ix)=ce*ar(ix)+se*ai(ix) 7422 tai(ix)=ce*ai(ix)-se*ar(ix) if(LDW)then do 82 j=i,min(np,i+nff-1) 82 call sadd(s0r,s0i,sff(j-i+1),tir,tii,qr,qi,tr,ti,sar,sai,sbr, 1 sbi,j,3*(j-1),talr,tali,tgr,tgi,tar,tai,sali,salr,sgi,sgr,sati, 2 satr) do 76 j=max(1,i-nff+1),i-1 76 call sadd(s0r,s0i,sff(i-j+1),tir,tii,qr,qi,tr,ti,sar,sai,sbr, 1 sbi,j,3*(j-1),talr,tali,tgr,tgi,tar,tai,sali,salr,sgi,sgr,sati, 2 satr) else call sadd(s0r,s0i,1.0d0,tir,tii,qr,qi,tr,ti,sar,sai,sbr, 1 sbi,i,3*(i-1),talr,tali,tgr,tgi,tar,tai,sali,salr,sgi,sgr,sati, 2 satr) endif 8 continue c 3 continue c c add the spectrum to the total: wi=WMIN-dd do 8997 i=1,NP wi=wi+dd c 1 correlation c 2 abs direct c 3 vcd direct c 4 Raman-180 c 5 ROA-ICP-180 st(i)=st(i)+s0r(i)**2 do 89971 ix=1,3 st( NP+i)=st( NP+i)+sbr(3*(i-1)+ix)*sbr(3*(i-1)+ix) 89971 st(2*NP+i)=st(2*NP+i)+sbr(3*(i-1)+ix)*sar(3*(i-1)+ix) SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 i9=9*(i-1) i7=27*(i-1) DO 6 II=1,3 DO 6 J=1,3 iti=II+3*(II-1) itj=J+3*(J-1) io=II+3*(J-1) c SAL0=SAL0+ALPHA(IQ,I,I)* ALPHA(IQ,J,J) SAL0=SAL0+salr(i9+iti)*salr(i9+itj)+ 1 sali(i9+iti)*sali(i9+itj) c SAL1=SAL1+ALPHA(IQ,I,J)* ALPHA(IQ,I,J) SAL1=SAL1+salr(i9+io )*salr(i9+io )+ 1 sali(i9+io )*sali(i9+io ) c SAG0=SAG0+ALPHA(IQ,I,I)* G(IQ,J,J) SAG0=SAG0+salr(i9+iti )*sgr(i9+itj )+ 1 sali(i9+iti )*sgi(i9+itj ) c SAG1=SAG1+ALPHA(IQ,I,J)* G(IQ,I,J) SAG1=SAG1+salr(i9+io )*sgr(i9+io )+ 1 sali(i9+io )*sgi(i9+io ) DO 6 K=1,3 DO 6 L=1,3 c SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0 itkl=K+3*(L-1)+9*(J-1) 6 SA1=SA1+EPS(II,K,L)*(salr(i9+io)*satr(i7+itkl)+ 1 sali(i9+io)*sati(i7+itkl))/3.0d0 c c 180 backscattering IL+IR st(3*NP+i)=st(3*NP+i)+7.0d0*SAL1+SAL0 c c 180 backscattering IR-IL 8997 st(4*NP+i)=st(4*NP+i)+8.0d0*(3.0d0*SAG1-SAG0+SA1) if(is.eq.nend(ipr))then write(6,91)ipr,nstart(ipr),nend(ipr) 91 format('Thread',i3,' spectra ',i4,' - ',i4,' added.') endif 899 continue C$OMP END PARALLEL 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/dble(nmd)/dt/dble(nsim) endif open(91,file='S.PRN') open(92,file='D.PRN') open(93,file='R.PRN') open(94,file='RAM.PRN') open(95,file='ROA.PRN') wi=WMIN-dd do 64 i=1,NP wi=wi+dd s1=st(i) s2=st(NP+i) s3=st(2*NP+i) s4=st(3*NP+i) sp=st(4*NP+i) if(LDW)then s2=s2*108.7d0*wi**2*c2 s3=s3*435.0d0*wi**3*c2 s4=s4*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp))) sp=sp*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp))) else s2=s2*108.7d0*wi*c2 s3=s3*435.0d0*wi**2*c2 s4=s4*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp))) sp=sp*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp))) endif write(91,494)wi,s1 write(92,494)wi,s2 write(93,494)wi,s3 write(94,494)wi,s4 64 write(95,494)wi,sp 494 format(g15.6,' ',g15.6) close(91) close(92) close(93) close(94) close(95) write(6,*)' Spectra written' write(6,6009)mclock()-it0 stop END c subroutine report(s) character*(*) s write(6,*)s stop 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,talr,tali,tgr,tgi,tar,tai,sali,salr,sgi,sgr,sati, 2 satr) implicit none real*8 s0r(*),s0i(*),sar(*),sai(*),sbr(*),sbi(*),con,qr(*), 1qi(*),tr(*),ti(*),tii,tir,talr(*),tali(*),tar(*),tai(*),tgr(*), 2tgi(*),sali(*),salr(*),sgr(*),sgi(*),sati(*),satr(*) integer*4 i2,i33,ix s0r(i2)=s0r(i2)+con*tir s0i(i2)=s0i(i2)+con*tii do 1 ix=1,3 sar(i33+ix)=sar(i33+ix)+con*qr(ix) sai(i33+ix)=sai(i33+ix)+con*qi(ix) sbr(i33+ix)=sbr(i33+ix)+con*tr(ix) 1 sbi(i33+ix)=sbi(i33+ix)+con*ti(ix) do 2 ix=1,9 salr(i33*3+ix)=salr(i33*3+ix)+con*talr(ix) sali(i33*3+ix)=sali(i33*3+ix)+con*tali(ix) sgr( i33*3+ix)=sgr( i33*3+ix)+con*tgr( ix) 2 sgi( i33*3+ix)=sgi( i33*3+ix)+con*tgi( ix) do 3 ix=1,27 satr(i33*9+ix)=satr(i33*9+ix)+con*tar(ix) 3 sati(i33*9+ix)=sati(i33*9+ix)+con*tai(ix) return end C SUBROUTINE READTTT(MX3,ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(MX3,9),G(3*MX3,9),A(3*MX3,27) logical lex inquire(file='FILE.TTT',exist=lex) if(.not.lex)return open(15,file='FILE.TTT') READ(15,*) READ(15,*)NATC if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file') 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,IXDUM,(ALPHA(IIND,I+3*(J-1)),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,IXDUM,(G(IIND,I+3*(J-1)),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,IXDUM,(A(IIND,I+3*(J-1)+9*(K-1)),K=1,3) CLOSE(15) RETURN END FUNCTION EPS(I,J,K) IMPLICIT INTEGER*4 (I-N) REAL*8 EPS 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)THEN EPS= 1.0D0 ELSE IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0D0 ENDIF else IF(J.EQ.1.AND.K.EQ.2)THEN EPS= 1.0D0 ELSE IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0D0 ENDIF endif endif RETURN END