program hoe implicit none integer*4 nt,nw,pr,ix,nit,nr,i,n,nat, 1ipa,ida,nbb,nproc,ipol,nc(3),ncp(2),lwt,nijj, 1niij,nijk,niijj,njjij,niiij,niijk,njjik,CHIR real*8 e0,gamma,wmin,wmax,dw,dt,bohr,ftresh,kk, 1cx(12),cpx(9),rlim,po(3),tt,temp real*8,allocatable::ff(:,:),r(:),P(:,:),A(:,:),fb(:), 1m(:),sp(:,:),f3(:,:,:),f4(:,:,:),v6(:,:), 1alpha(:,:,:),xalpha(:,:),ar(:,:,:),ai(:,:,:),bmr(:,:,:), 1bmi(:,:,:),cijjb(:),ciii(:),ciijb(:),cijkb(:), 1diijjb(:),diiii(:),diiijb(:),diijkb(:),djjijb(:),djjikb(:), 1qq(:,:) integer*4,allocatable::jb(:),ns(:),ne(:),q(:),iu(:), 1j_ijj(:),n_ijj(:),j_iij(:),n_iij(:),j_ijk(:),n_ijk(:), 1k_ijk(:),j_iijj(:),n_iijj(:),j_iiij(:),n_iiij(:),j_iijk(:), 1n_iijk(:),k_iijk(:),n_jjij(:),j_jjij(:),j_jjik(:),k_jjik(:), 1n_jjik(:) logical lgeo,lw0,lcage,lplane,lram,l3,l4,lijk,lpoint,quad, 1fcoord,fft,sft,lvcd,usea,useg,cram,lloc bohr=0.529177d0 c write(*,1010) 1010 FORMAT(' Harmonic Oscillators (Molecules) in Electric Field',/, 1 ' RP propagation',/,/, 1 ' Input: FILE.X geometry',/, 1 ' FILE.TEN dipole derivatives',/, 1 ' FILE.FC force field',/, 1 ' FTRY.INP atomic masses',/, 1 ' HOE.OPT options',/,/, 1 ' Output: ABS.PRN IR absorption',/, 1 ' VCD.PRN VCD',/) call readopt(e0,gamma,nt,wmin,wmax,nw,pr,nr,nit,dt,ipa,ida, 1ftresh,nproc,kk,lgeo,lw0,ipol,lcage,cx,nc,lplane,cpx,ncp, 1lvcd,lram,lwt,l3,l4,lijk,rlim,lpoint,po,CHIR,quad,fcoord,fft, 1sft,tt,usea,useg,cram,temp,lloc,lrt,rtc) if(nproc.ne.0)call omp_set_num_threads(nproc) open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat n=3*nat allocate(r(n),m(nat),ff(n,n),P(n,3),q(nat),qq(n,6), 1A(n,3),sp(nw,8),v6(6,n)) do 37 i=1,nat read(9,*)q(i),(r(ix+3*(i-1)),ix=1,3) do 37 ix=1,3 37 r(ix+3*(i-1))=r(ix+3*(i-1))/bohr close(9) call setv6(n,v6,r) if(ipol.ne.0)then allocate(ar(ipol,3,3),ai(ipol,3,3),bmr(ipol,3,3),bmi(ipol,3,3), 1 iu(nat),alpha(ipol,3,3),xalpha(ipol,3)) c load all polarizabilities: call ldal(alpha,xalpha,ipol) c find closest unit to each atom: call setiu(ipol,nat,iu,r,xalpha) else allocate(alpha(1,3,3),xalpha(1,3),ar(1,3,3),ai(1,3,3),bmr(1,3,3), 1 bmi(1,3,3),iu(1)) endif call READM(nat,m) call READFF(n,ff) call COUNTFF(n,ff,ftresh,nbb) allocate(fb(nbb),jb(nbb),ns(n),ne(n)) call BUFFF(n,nbb,ns,ne,jb,ff,fb,ftresh) deallocate(ff) if(l3)then allocate(f3(n,n,n)) call READF34(n,f3,'SCR.33','cubic') call COUNTF3(n,f3,ftresh,lijk,nijj,niij,nijk) allocate(ciii(n),cijjb(nijj),ciijb(niij),j_ijj(nijj),n_ijj(n), 1 j_iij(niij),n_iij(n),n_ijk(n)) if(lijk)then allocate(cijkb(nijk),j_ijk(nijk),k_ijk(nijk)) else allocate(cijkb(1),j_ijk(1),k_ijk(1)) endif call BUFF3(lijk,n,n_ijj,n_iij,n_ijk,j_iij,j_ijj, 1 j_ijk,k_ijk,f3,ftresh,ciii,cijjb,ciijb,cijkb) deallocate(f3) else allocate(ciii(1),cijjb(1),ciijb(1),j_ijj(1),n_ijj(1), 1 j_iij(1),n_iij(1),n_ijk(1),cijkb(1),j_ijk(1),k_ijk(1)) endif if(l4)then allocate(f4(n,n,n)) call READF34(n,f4,'SCR.44','quartic') call COUNTF4(n,f4,ftresh,lijk,niijj,njjij,niiij,niijk,njjik) allocate(diiii(n),diijjb(niijj),diiijb(niiij),djjijb(njjij), 1 j_iijj(niijj),n_iijj(n),j_iiij(niiij),n_iiij(n),j_jjij(njjij), 1 n_jjij(n),n_iijk(n)) if(lijk)then allocate(diijkb(niijk),j_iijk(niijk), 1 k_iijk(niijk),djjikb(njjik),j_jjik(njjik),n_jjik(n), 1 k_jjik(njjik)) else allocate(diijkb(1),j_iijk(1), 1 k_iijk(1),djjikb(1),j_jjik(1),n_jjik(1),k_jjik(1)) endif call BUFF4(lijk,n,n_iijj,n_iiij,n_jjij,n_iijk,n_jjik,j_iijj, 1 j_iiij,j_jjij,j_iijk,k_iijk,j_jjik,k_jjik, 1 f4,ftresh,diiii,diijjb,djjijb,diiijb,diijkb,djjikb) deallocate(f4) else allocate(diiii(1),diijjb(1),diiijb(1),djjijb(1), 1 j_iijj(1),n_iijj(1),j_iiij(1),n_iiij(1),j_jjij(1), 1 n_jjij(1),n_iijk(1),diijkb(1),j_iijk(1), 1 k_iijk(1),djjikb(1),j_jjik(1),n_jjik(1),k_jjik(1)) endif call READTEN(nat,P,A,.false.,r) call READQD(n,qq,r,P,quad) if(nw.eq.1)then dw=0.0d0 else dw=(wmax-wmin)/dble(nw-1) endif if(lcage)call wrf(nat,ipol,pr,nc,cx,q,r,xalpha,wmin,alpha,dw,rlim) if(lplane)call wrp(nat,ipol,pr,ncp,cpx,q,r,xalpha,wmin,alpha,dw, 1rlim) if(lpoint)call wpo(nat,ipol,xalpha,alpha,qq,P,A,po,rlim) if(fft)call runft(nt,ipol,CHIR,nat,n,po,alpha,xalpha,ar,ai,bmr, 1bmi,rlim,tt,fb,ns,ne,m,lw0,lwt,kk,P,A,gamma,dt,lgeo,q,ipa,nr, 1ida,fcoord,qq,e0,iu,jb,r,wmin,wmax) if(sft)call ft(nt,nw,ipol,CHIR,nat,n,po,alpha,xalpha,ar,ai,bmr, 1bmi,rlim,tt,fb,ns,ne,m,lw0,lwt,kk,P,A,gamma,dt,lgeo,q,ipa,nr, 1ida,fcoord,qq,e0,iu,jb,r,wmin,dw) if(lvcd)call runvcd(nt,ipol,CHIR,nat,n,nit,nw,pr, 1po,alpha,xalpha,ar,ai,bmr,lloc,crt,lrt,v6, 1bmi,rlim,tt,fb,ns,ne,m,lw0,lwt,kk,P,A,gamma,dt,lgeo,q,ipa,nr, 1ida,fcoord,qq,e0,iu,jb,r,wmin,sp,dw,l3,l4,lijk, 1ciii, cijjb,ciijb,cijkb,j_ijj,j_iij,j_ijk,k_ijk, n_ijj,n_iij, 1n_ijk,n_iijj,n_jjij,n_iiij,n_iijk,j_iijj,j_iiij,j_jjij, 1j_iijk,k_iijk,n_jjik,j_jjik,k_jjik,diiii,diijjb,djjijb,diiijb, 1diijkb,djjikb) if(lram)then if(cram.or.ipol.ne.0)then call ram(nt,ipol,CHIR,nat,n,nit,nw,pr,po,alpha,xalpha,ar,ai, 1 bmr,bmi,rlim,fb,ns,ne,m,lw0,lwt,kk,gamma,dt, 1 lgeo,q,ipa,nr,ida,fcoord,e0,iu,jb,r,wmin,sp,dw,usea,useg,temp) else call ram2(nt,CHIR,nat,n,nit,nw,pr,fb,ns,ne,m,lw0,lwt,kk, 1 gamma,dt,lgeo,q,ipa,nr,ida,e0,jb,r,wmin,sp,dw,usea,useg,temp, 1 lloc,crt,lrt,v6) endif endif end subroutine ram2(nt0,CHIR,nat,n,nit,nw,pr,fb,ns,ne,m,lw0,lwt,kk, 1gamma,dt,lgeo,q,ipa,nr,ida,e0,jb,r,wmin,sp,dw,usea,useg,temp, 1lloc,crt,lrt,v6) implicit none integer*4 nt,ip,id,i2,i3,CHIR,n,it,io,jj,ns(*),ne(*),lwt, 1i,nat,q(*),ipa,nr,ida,nfiles,iw,jb(*),ix,iwc,nt0, 1nit,nw,pr,ipx real*8 t,fb(*),m(*),ma,kk,gamma,gdt,r(*),dw,clight,dt,CM,dp,w,fc, 1e0,sip,kv,temp,f,tt,dt2,sp(nw,8),pil,t3,swt,cwt,df,cf,z2,wmin, 1di(2,2,3),w532,crt,v6(6,n) real*8,allocatable::d0(:),d1(:),fme(:),al(:,:,:),go(:,:,:), 1aa(:,:,:,:),ql(:) logical lw0,lgeo,usea,useg,lloc,lrt parameter (nfiles=12) character*7 fn(nfiles) data fn/'PNR.PRN','ANR.PRN','ATL.PRN', 1 'ADI.PRN','RAM.PRN','ROA.PRN', 1 'RAX.PRN','RAY.PRN','RAZ.PRN', 1 'ROX.PRN','ROY.PRN','ROZ.PRN'/ write(6,640) 640 format(/,' ram2: Raman spectra, with real fields',/) nt=nt0 clight=137.03604d0 z2=2.0d0 gdt=gamma*dt dt2=dt**2 CM =219474.63d0 cf=1.0d0+gamma*dt/2.0d0 df=1.0d0-gamma*dt/2.0d0 pil=4.0d0*datan(1.0d0) w532=1.0d7/532.0d0/CM kv=w532/clight allocate(al(n,3,3),go(n,3,3),aa(n,3,3,3),d0(n),d1(n),fme(n)) CALL READTTT(n,al,aa,go) if(lloc)then allocate(ql(nat)) ql=0.0d0 endif do 8 io=1,nfiles 8 open(10+io,file=fn(io)) iwc=0 c loop over frequencies: C$OMP Parallel Do Schedule(Dynamic) Default(Shared) C$OMP+ Private(w,iw,ip,id,d0,d1,t,it,tt,ipx, C$OMP+ sip,cwt,swt,i2,i3, C$OMP+ i,t3,ma,ix,io,f,jj,dp,di,fme) do 1 iw=1,nw w=wmin+dw*dble(iw-1) write(6,608)iw 608 format(i4,$) iwc=iwc+1 if(mod(iwc,20).eq.0)write(6,*) di=0.0d0 c loop over polarizations, RCPL, LCPL, scattered light: do 2 ip=1,2 sip=dble(-3+2*ip) c for CHIR=0 the same light if(CHIR.eq.0)sip=1.0d0 c loop over x y z propagation directions do 25 id=1,3 call setij(id,i2,i3) c loop over Y=i2 polarized or Z=i3 polarized, excitation light: do 26 ipx=1,2 c initialize deviations: d0=0.0d0 d1=0.0d0 c loop over time steps: t=-dt do 3 it=1,nt t=t+dt swt=e0*sin(w*t) cwt=e0*cos(w*t) c mechanical force= - FF .d: fme=0.0d0 do 36 io=1,n do 36 jj=ns(io),ne(io) 36 fme(io)=fme(io)-fb(jj)*d0(jb(jj)) c if required avoid PES slips: if(lw0)call aps(n,d0,fme,lwt) do 39 i=1,nat ma=m(i) do 39 ix=1,3 io=ix+3*(i-1) f=fme(io) c arbitrary force to keep atoms in place f=f-kk*d0(io) c force for this coordinate from electric field of the light F = P.E: if(ipx.eq.1)then c y-polarization f=f+al(io,i2,i2)*cwt-sip*al(io,i2,i3)*swt else c z-polarization f=f+al(io,i2,i3)*cwt-sip*al(io,i3,i3)*swt endif c for light polarization, add the magnetic and quadrupole forces c CHIR=0: we calculate I(with magnetic moment)-I(without) c instead of I(pol+) - I(pol-), becasue of numerical stability c CHIR=1 calculate directly I(pol+) - I(pol-)fort.15 c CHIR=2 calculate I(+,R-enantiomer)-I(+,S-enantiomer) t3=0.0d0 if(useg)then c G (E'Im(B) - Im(E)B'), G - last index magnetic if(ipx.eq.1)then c y-polarization, B along z c -p(Gyy-Gzz)cos(wt) + 2 Gyz sin(wt) t3=t3-sip*kv*(go(io,i2,i2)-go(io,i3,i3))*cwt 1 +kv*(go(io,i2,i3)+go(io,i2,i3))*swt else c z-polarization, B along y c (Gzz-Gyy)sin(wt) - 2 p Gzy cos(wt) t3=t3 +kv*(go(io,i3,i3)-go(io,i2,i2))*swt 1 -sip*kv*(go(io,i3,i2)+go(io,i3,i2))*cwt endif endif if(usea)then if(ipx.eq.1)then c y-polarization c -p(1/3)(Ay,xz + Az,xy) cos(wt) - 2 (1/3) Ay,xy sin(wt) t3=t3-sip*kv*(aa(io,i3,id,i2)+aa(io,i2,id,i3))*cwt/3.0d0 1 -kv*(aa(io,i2,i2,id)+aa(io,i2,i2,id))*swt/3.0d0 else c -(1/3)(Az,xy + Ay,xz) sin(wt) -p 2 (1/3) Az,xz cos(wt) t3=t3 -kv*(aa(io,i2,id,i3)+aa(io,i3,id,i2))*swt/3.0d0 1 -sip*kv*(aa(io,i3,id,i3)+aa(io,i3,i3,id))*cwt/3.0d0 endif endif tt=0.0d0 if((CHIR.eq.0.and.ip.eq.1).or.CHIR.eq.1)tt=t3 c XX YY ZZ XY XZ YZ c 1 2 3 4 5 6 if(CHIR.eq.2)then tt=t3 if(ip.eq.2)tt=-tt endif f=f+tt c get coordinate io in time t+dt: dp=(z2*d0(io)-d1(io)*df+f*dt2/ma)/cf c check that the coordinate does not become too big: if(dabs(dp).gt.0.5d0)call wrc(n,i,ix,io,iw,ip,id,it, 1d0,dp,f,0.0d0,fme(io),kk,lgeo,nat,q,r) c write requested path: if(iw.eq.pr.and.ip.eq.ipa.and.io.eq.nr.and.id.eq.ida 1.and.ipx.eq.1)write(11,100)it,d0(io) 100 format(i7,e14.5) c for t>nit*dt gather dissipated energy for each polarization and direction: if(it.gt.nit)di(ip,ipx,id)=di(ip,ipx,id) 1+gdt*ma*((dp-d1(io))/dt/2.0d0)**2 c energy dissipated on atom, just relative: if(lloc)ql(i)=ql(i)+ma*((dp-d1(io)))**2 d1(io)=d0(io) 39 d0(io)=dp c check that molecule is not translated or rotated if(lrt)call chkrt(n,d0,v6,crt) 3 continue c t 26 continue c ipx, planar polarizations 25 continue c id, directions 2 continue c ip,with and without chirality fc=3.1136d-11 / (dble(nt-nit)*dt*pil*e0**4)/(w*CM) 1/(1.0d0-exp(-1.440d0*w*CM/temp)) do 10 ix=1,3 c Raman and ROA for propagation along ix: sp(iw,2+ix)=0.5d0*fc*(di(2,1,ix)+di(2,2,ix)+di(1,1,ix)+di(1,2,ix)) 10 sp(iw,5+ix)=4.0d0*fc*(di(2,1,ix)+di(2,2,ix)-di(1,1,ix)-di(1,2,ix)) 1*1.0d4 if(CHIR.eq.2)then do 101 ix=1,3 101 sp(iw,5+ix)=sp(iw,5+ix)/2.0d0 endif c average strengths: sp(iw,1)=(sp(iw,3)+sp(iw,4)+sp(iw,5))/3.0d0 1 sp(iw,2)=(sp(iw,6)+sp(iw,7)+sp(iw,8))/3.0d0 c light frequency w=wmin-dw do 91 iw=1,nw w=w+dw do 91 ix=1,8 91 write(14+ix,600)w*CM,sp(iw,ix) 600 format(f12.2,e15.6) do 9 io=1,nfiles 9 close(10+io) if(lloc)call wrq(nat,ql,q,r) stop end subroutine ram(nt0,ipol,CHIR,nat,n,nit,nw,pr,po,alpha,xalpha, 1ar,ai,bmr,bmi,rlim,fb,ns,ne,m,lw0,lwt,kk,gamma,dt,lgeo,q,ipa, 1nr,ida,fcoord,e0,iu,jb,r,wmin,sp,dw,usea,useg,temp) implicit none integer*4 nt,ip,id,ipol,CHIR,n,it,io,jj,ns(*),ne(*),lwt, 1i,nat,q(*),ipa,nr,ida,nfiles,iu(*),iw,jb(*),ix,iwc,nt0, 1nit,nw,pr,ipx,j,k,l,i2,i3 real*8 rlim,po(3),t,fb(*),m(*),ma,kk,gamma,gdt,r(*),dw,sip, 1alpha(ipol,3,3),xalpha(ipol,3),dt,CM,dp,w,fc,e0,temp, 1cwt,swt,f,tt,dt2,sp(nw,8),pil,t3,ar(ipol,3,3),ai(ipol,3,3), 1akr(ipol,3,3),aki(ipol,3,3),bkr(ipol,3,3),bki(ipol,3,3), 1ckr(ipol,3,3,3),cki(ipol,3,3,3),cr(ipol,3,3,3),ci(ipol,3,3,3), 1bmr(ipol,3,3),bmi(ipol,3,3),df,cf,z2,wmin,di(2,2,3), 1Er(3) , Ei(3) , Br(3) , Bi(3), Gr(3,3), Gi(3,3), 1Epr(3), Epi(3), Bpr(3), Bpi(3), Gpr(3,3),Gpi(3,3), 1E(3) , Eii(3) , B(3) , Bii(3), G(3,3) , Gii(3,3), 1Ep(3) , Epii(3), Bp(3) , Bpii(3),Gp(3,3) ,Gpii(3,3), 1Etr(3), Eti(3), Btr(3), Bti(3), 1Et(3) , Etii(3), Bt(3), Btii(3) real*8,allocatable::d0(:),d1(:),fme(:),al(:,:,:),go(:,:,:), 1aa(:,:,:,:) logical lw0,lgeo,fcoord,usea,useg parameter (nfiles=12) character*7 fn(nfiles) data fn/'PNR.PRN','ANR.PRN','ATL.PRN', 1 'ADI.PRN','RAM.PRN','ROA.PRN', 1 'RAX.PRN','RAY.PRN','RAZ.PRN', 1 'ROX.PRN','ROY.PRN','ROZ.PRN'/ write(6,640) 640 format(/,' ram: Raman spectra, with complex fields',/) nt=nt0 z2=2.0d0 gdt=gamma*dt dt2=dt**2 CM =219474.63d0 cf=1.0d0+gamma*dt/2.0d0 df=1.0d0-gamma*dt/2.0d0 pil=4.0d0*datan(1.0d0) allocate(al(n,3,3),go(n,3,3),aa(n,3,3,3),d0(n),d1(n),fme(n)) CALL READTTT(n,al,aa,go) do 8 io=1,nfiles 8 open(10+io,file=fn(io)) iwc=0 c Excitation field Er Ei Br Bi Gr Gi c Scattered field Epr Epi Bpr Bpi Gpr Gpi c Redressed excitation E Eii B Bii G Gii c Redressed scattered Ep Epi Bp Bpii Gp Gpii c Time derivatives: c Excitation field Etr Eti Btr Bti c Scattered field c Redressed excitation Et Etii Bt Btii c Redressed scattered C$OMP Parallel Do Schedule(Dynamic) Default(Shared) C$OMP+ Private(w,iw,ip,id,d0,d1,t,it,sip,tt,ipx,cwt,swt,i2,i3, C$OMP+ i,j,k,l,t3,ma,ix,io,f,jj,dp,di,fme,ar,ai,bmr,bmi, C$OMP+ akr,aki,bkr,bki,cr,ci,ckr,cki, C$OMP+ Er ,Ei , Br ,Bi, Gr, Gi, Epr , Epi , Bpr ,Bpi , Gpr,Gpi, C$OMP+ E ,Eii ,B ,Bii, G , Gii, Ep , Epii, Bp ,Bpii, Gp ,Gpii, C$OMP+ Etr,Eti ,Btr,Bti, C$OMP+ Et ,Etii,Bt, Btii ) do 1 iw=1,nw w=wmin+dw*dble(iw-1) write(6,608)iw 608 format(i4,$) iwc=iwc+1 if(mod(iwc,20).eq.0)write(6,*) di=0.0d0 c loop over polarizations, RCPL, LCPL, scattered light: do 2 ip=1,2 sip=dble(-3+2*ip) c for CHIR=0 the same light if(CHIR.eq.0)sip=1.0d0 c loop over polarizations, X,Y, excitation light: do 2 ipx=1,2 c loop over x y z propagation directions do 25 id=1,3 call setij(id,i2,i3) c set relation between group and external field: c E(l,t) = a Eext(0,t), B(l,t) = bm Eext(0,t) if(ipol.ne.0)then if(fcoord)then c use field at point po: c scattered light ~ exp[i(k'x-wt)] call calcx(ipol,id,w,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,po,1.0d0) c excitation light ~ exp[i(-kx-wt)] call calcx(ipol,id,w,alpha,xalpha,akr,aki,bkr,bki,ckr,cki, 1 rlim,ip,CHIR,po,-1.0d0) else c use field at closets unit to each atom: call calca(ipol,id,w,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,1.0d0) call calca(ipol,id,w,alpha,xalpha,akr,aki,bkr,bki,ckr,cki, 1 rlim,ip,CHIR,-1.0d0) endif endif c initialize deviations: d0=0.0d0 d1=0.0d0 c at t=0: c planar excitation light: call setfieldram(id,Er, Ei,Etr,Eti, Br, Bi,Btr,Bti, Gr, Gi, ipx) c chiral scattered light: call fieldramp(id,Epr,Epi,Bpr,Bpi,Gpr,Gpi,sip) c loop over time steps: t=-dt do 3 it=1,nt t=t+dt swt=e0*sin(w*t) cwt=e0*cos(w*t) c mechanical force= - FF .d: fme=0.0d0 do 36 io=1,n do 36 jj=ns(io),ne(io) 36 fme(io)=fme(io)-fb(jj)*d0(jb(jj)) c if required avoid PES slips: if(lw0)call aps(n,d0,fme,lwt) do 39 i=1,nat ma=m(i) c E and B fields for this atom and this time: if(ipol.ne.0)then c E(red)=a . E, B(red)=b . E call redress(iu(i),ipol,ar, ai, bmr,bmi,Epr,Epi,Ep,Epii,Bp,Bpii) call redress(iu(i),ipol,akr,aki,bkr,bki,Er ,Ei, E, Eii, B, Bii) call redress(iu(i),ipol,akr,aki,bkr,bki,Etr,Eti,Et,Etii,Bt,Btii) Gp=Gpr Gpii=Gpi else E=Er Eii=Ei B=Br Bii=Bi Et=Etr Etii=Eti Bt=Btr Btii=Bti G=Gr Gii=Gi Ep=Epr Epii=Epi Bp=Bpr Bpii=Bpi Gp=Gpr Gpii=Gpi endif do 39 ix=1,3 io=ix+3*(i-1) if(iw.eq.pr.and.it.eq.1.and.id.eq.ida.and.io.eq.nr)call wre(w*CM, 1ip,id,Er,Ei,Br,Bi,E,Eii,B,Bii,e0) f=fme(io) c arbitrary force to keep atoms in place f=f-kk*d0(io) c force for Raman scattering: c a_io,j,k Re(Ej)Re(E'k) c -iwt -iw't c product in time t for complex numbers U=u e and V=v e : c Re(U)Re(V) ~ [Re(u)Re(v)+Im(u)Im(v)]cos(Dt)/2 c + [Im(u)Re(v)-Re(u)Im(v)]sin(Dt)/2, D = w - w' do 391 j=1,3 do 391 k=1,3 391 f=f+al(io,j,k)*( (E(j)*Ep(k)+Eii(j)*Epii(k))*cwt 1 +(Eii(j)*Ep(k)-E(j)*Epii(k))*swt ) c for light polarization, add the magnetic and quadrupole forces c CHIR=0: we calculate I(with magnetic moment)-I(without) c instead of I(pol+) - I(pol-), becasue of numerical stability c CHIR=1 calculate directly I(pol+) - I(pol-)fort.15 c CHIR=2 calculate I(+,R-enantiomer)-I(+,S-enantiomer) t3=0.0d0 if(useg)then c G ( Bt E' - Et B') , G - last index magnetic, c suppose that we have G/w already do 392 j=1,3 do 392 k=1,3 392 t3=t3+go(io,j,k)*( (Bt(k)*Ep(j)+Btii(k)*Epii(j))*cwt 1 +(Btii(k)*Ep(j)-Bt(k)*Epii(j))*swt 1 -(Et(j)*Bp(k)+Etii(j)*Bpii(k))*cwt 1 -(Etii(j)*Bp(k)-Et(j)*Bpii(k))*swt ) endif if(usea)then c (1/3) A ( E Grad' + Grad E' ), A .. first index electric do 393 j=1,3 do 393 k=1,3 do 393 l=1,3 393 t3=t3+aa(io,j,k,l)*( (E(j)*Gp(k,l)+Eii(j)*Gpii(k,l))*cwt 1 +(Eii(j)*Gp(k,l)-E(j)*Gpii(k,l))*swt 1 +(G(k,l)*Ep(j)+Gii(k,l)*Epii(j))*cwt 1 +(Gii(k,l)*Ep(j)-G(k,l)*Epii(j))*swt)/3.0d0 endif tt=0.0d0 if((CHIR.eq.0.and.ip.eq.1).or.CHIR.eq.1)tt=t3 c XX YY ZZ XY XZ YZ c 1 2 3 4 5 6 if(CHIR.eq.2)then tt=t3 if(ip.eq.2)tt=-tt endif f=f+tt c get coordinate io in time t+dt: dp=(z2*d0(io)-d1(io)*df+f*dt2/ma)/cf c check that the coordinate does not become too big: if(dabs(dp).gt.0.5d0)call wrc(n,i,ix,io,iw,ip,id,it, 1d0,dp,f,0.0d0,fme(io),kk,lgeo,nat,q,r) c write requested path: if(iw.eq.pr.and.ip.eq.ipa.and.io.eq.nr.and.id.eq.ida 1.and.ipx.eq.1)write(11,100)it,d0(io) 100 format(i7,e14.5) c for t>nit*dt gather dissipated energy for each polarization and direction: if(it.gt.nit)di(ip,ipx,id)=di(ip,ipx,id) 1+gdt*ma*((dp-d1(io))/dt/2.0d0)**2 d1(io)=d0(io) 39 d0(io)=dp 3 continue c t 25 continue c id, directions 2 continue c ip,with and without A,ipx fc=3.1136d-11 / (dble(nt-nit)*dt*pil*e0**4)/(w*CM) 1/(1.0d0-exp(-1.440d0*w*CM/temp)) do 10 ix=1,3 c Raman and ROA for propagation along ix: sp(iw,2+ix)=0.5d0*fc*(di(2,1,ix)+di(2,2,ix)+di(1,1,ix)+di(1,2,ix)) 10 sp(iw,5+ix)=4.0d0*fc*(di(2,1,ix)+di(2,2,ix)-di(1,1,ix)-di(1,2,ix)) 1*1.0d4 if(CHIR.eq.2)then do 101 ix=1,3 101 sp(iw,5+ix)=sp(iw,5+ix)/2.0d0 endif c average strengths: sp(iw,1)=(sp(iw,3)+sp(iw,4)+sp(iw,5))/3.0d0 1 sp(iw,2)=(sp(iw,6)+sp(iw,7)+sp(iw,8))/3.0d0 c light frequency w=wmin-dw do 91 iw=1,nw w=w+dw do 91 ix=1,8 91 write(14+ix,600)w*CM,sp(iw,ix) 600 format(f12.2,e15.6) do 9 io=1,nfiles 9 close(10+io) stop end subroutine runvcd(nt0,ipol,CHIR,nat,n,nit,nw,pr, 1po,alpha,xalpha,ar,ai,bmr,lloc,crt,lrt,v6, 1bmi,rlim,tt,fb,ns,ne,m,lw0,lwt,kk,P,A,gamma,dt,lgeo,q,ipa,nr, 1ida,fcoord,qq,e0,iu,jb,r,wmin,sp,dw,l3,l4,lijk, 1ciii, cijjb,ciijb,cijkb,j_ijj,j_iij,j_ijk,k_ijk, n_ijj,n_iij, 1n_ijk,n_iijj,n_jjij,n_iiij,n_iijk,j_iijj,j_iiij,j_jjij, 1j_iijk,k_iijk,n_jjik,j_jjik,k_jjik,diiii,diijjb,djjijb,diiijb, 1diijkb,djjikb) implicit none integer*4 nt,ip,id,ipol,CHIR,n,it,ipu,io,jj,ns(*),ne(*),lwt, 1i,nat,q(*),ipa,nr,ida,nfiles,iu(*),iw,jb(*),ix,iwc,nt0, 1nit,nw,pr, 1j_ijj(*),j_iij(*),j_ijk(*),k_ijk(*), n_ijj(*),n_iij(*), 1n_ijk(*),n_iijj(*),n_jjij(*),n_iiij(*),n_iijk(*), 1j_iijj(*),j_iiij(*),j_jjij(*),j_iijk(*),k_iijk(*),n_jjik(*), 2j_jjik(*),k_jjik(*) real*8 rlim,po(3),t,fb(*),m(*),ma,kk,gamma,gdt,r(*),dw,crt, 1alpha(ipol,3,3),xalpha(ipol,3),E(3),B(3),Ei(3),Eii(3), 1dt,CM,qq(n,6),dp,w,fc,e0,Er(3),Br(3),Bii(3),,v6(6,n), 1f,P(n,3),A(n,3),tt,dt2,sp(nw,8),pil,Bi(3),cr(ipol,3,3,3), 1ci(ipol,3,3,3), 1ar(ipol,3,3),ai(ipol,3,3),bmr(ipol,3,3),bmi(ipol,3,3),df,cf,z2, 1wmin,di(2,3),Gr(3,3),Gi(3,3), 1ciii(*), cijjb(*),ciijb(*),cijkb(*),diiii(*),diijjb(*),djjijb(*), 1diiijb(*),diijkb(*),djjikb(*) real*8,allocatable::d0(:),d1(:),fme(:),ql(:) logical lw0,lgeo,fcoord,l3,l4,lijk,lloc,lrt parameter (nfiles=12) character*7 fn(nfiles) data fn/'PNR.PRN','ANR.PRN','ATL.PRN', 1 'ADI.PRN','ABS.PRN','VCD.PRN', 1 'ABX.PRN','ABY.PRN','ABZ.PRN', 1 'VCX.PRN','VCY.PRN','VCZ.PRN'/ do 8 io=1,nfiles 8 open(10+io,file=fn(io)) nt=nt0 z2=2.0d0 gdt=gamma*dt dt2=dt**2 CM =219474.63d0 cf=1.0d0+gamma*dt/2.0d0 df=1.0d0-gamma*dt/2.0d0 pil=4.0d0*datan(1.0d0) allocate(d0(n),d1(n),fme(n)) if(lloc)then allocate(ql(nat)) ql=0.0d0 endif iwc=0 c loop over light frequencies C$OMP Parallel Do Schedule(Dynamic) Default(Shared) C$OMP+ Private(w,iw,ip,id,d0,d1,t,it,ipu,tt, C$OMP+ i,ma,ix,io,f,jj,dp,di,fme,ar,ai,bmr,bmi, C$OMP+ Er,Ei,Br,Bi,E,B,Eii,Bii,Gr,Gi) do 1 iw=1,nw w=wmin+dw*dble(iw-1) write(6,608)iw 608 format(i4,$) iwc=iwc+1 if(mod(iwc,20).eq.0)write(6,*) di=0.0d0 c loop over polarizations: do 2 ip=1,2 c loop over x y z propagation directions do 25 id=1,3 c set relation between group and external field: c E(l,t) = a Eext(0,t), B(l,t) = bm Eext(0,t) if(ipol.ne.0)then if(fcoord)then c use field at point po: call calcx(ipol,id,w,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,po,1.0d0) else c use field at closets unit to each atom: call calca(ipol,id,w,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,1.0d0) endif endif c initialize deviations: d0=0.0d0 d1=0.0d0 c loop over time steps: t=-dt do 3 it=1,nt t=t+dt ipu=1 if(CHIR.eq.1)ipu=ip c light at t and r=0: call setfield(id,w,t,e0,Er,Ei,Br,Bi,Gr,Gi,ipu) c mechanical force= - FF .d: fme=0.0d0 do 36 io=1,n do 36 jj=ns(io),ne(io) 36 fme(io)=fme(io)-fb(jj)*d0(jb(jj)) if(l3)call add3(n,d0,fme,ciii, cijjb,ciijb,cijkb, 1j_ijj,j_iij,j_ijk,k_ijk, n_ijj,n_iij, n_ijk,lijk) if(l4)call add4(lijk,n,d0,fme,n_iijj,n_jjij,n_iiij,n_iijk, 1j_iijj,j_iiij,j_jjij,j_iijk,k_iijk,n_jjik,j_jjik,k_jjik, 1diiii,diijjb,djjijb,diiijb,diijkb,djjikb) c if required avoid PES slips: if(lw0)call aps(n,d0,fme,lwt) do 39 i=1,nat ma=m(i) c E and B fields for this atom and this time: if(ipol.ne.0)then call redress(iu(i),ipol,ar,ai,bmr,bmi,Er,Ei,E,Eii,B,Bii) else E=Er Eii=Ei B=Br Bii=Bi endif do 39 ix=1,3 io=ix+3*(i-1) if(iw.eq.pr.and.it.eq.1.and.id.eq.ida.and.io.eq.nr)call wre(w*CM, 1ip,id,Er,Ei,Br,Bi,E,Eii,B,Bii,e0) f=fme(io) c arbitrary force to keep atoms in place f=f-kk*d0(io) c force for this coordinate from electric field of the light F = P.E: f=f+P(io,1)*E(1)+P(io,2)*E(2)+P(io,3)*E(3) c for light polarization, add the magnetic and quadrupole forces c CHIR=0: we calculate I(with magnetic moment)-I(without) c instead of I(pol+) - I(pol-), becasue of numerical stability c CHIR=1 calculate directly I(pol+) - I(pol-) c CHIR=2 calculate I(+,R-enantiomer)-I(+,S-enantiomer) tt=0.0d0 if((CHIR.eq.0.and.ip.eq.1).or.CHIR.eq.1) 1tt=-w*(A(io,1)*Bii(1)+A(io,2)*Bii(2)+A(io,3)*Bii(3)) 1 + (Gr(1,1)*qq(io,1)+Gr(1,2)*qq(io,4)+Gr(1,3)*qq(io,5) 1 + Gr(2,1)*qq(io,4)+Gr(2,2)*qq(io,2)+Gr(2,3)*qq(io,6) 1 + Gr(3,1)*qq(io,5)+Gr(3,2)*qq(io,6)+Gr(3,3)*qq(io,3))/6.0d0 c XX YY ZZ XY XZ YZ c 1 2 3 4 5 6 if(CHIR.eq.2)then tt=-w*(A(io,1)*Bii(1)+A(io,2)*Bii(2)+A(io,3)*Bii(3)) 1 + (Gr(1,1)*qq(io,1)+Gr(1,2)*qq(io,4)+Gr(1,3)*qq(io,5) 1 + Gr(2,1)*qq(io,4)+Gr(2,2)*qq(io,2)+Gr(2,3)*qq(io,6) 1 + Gr(3,1)*qq(io,5)+Gr(3,2)*qq(io,6)+Gr(3,3)*qq(io,3))/6.0d0 if(ip.eq.2)tt=-tt endif f=f+tt c get coordinate io in time t+dt: dp=(z2*d0(io)-d1(io)*df+f*dt2/ma)/cf c check that the coordinate does not become too big: if(dabs(dp).gt.0.5d0)call wrc(n,i,ix,io,iw,ip,id,it, 1d0,dp,f,P(io,1)*E(1)+P(io,2)*E(2)+P(io,3)*E(3), 1fme(io),kk,lgeo,nat,q,r) c write requested path: if(iw.eq.pr.and.ip.eq.ipa.and.io.eq.nr.and.id.eq.ida) 1write(11,100)it,d0(io) 100 format(i7,e14.5) c for t>nit*dt gather dissipated energy for each polarization and direction: if(it.gt.nit)di(ip,id)=di(ip,id)+gdt*ma*((dp-d1(io))/dt/2.0d0)**2 c energy absorbed on atom, just relative: if(lloc)ql(i)=ql(i)+ma*((dp-d1(io)))**2 d1(io)=d0(io) 39 d0(io)=dp if(lrt)call chkrt(n,d0,v6,crt) 3 continue c t 25 continue c id, directions 2 continue c ip,with and without A fc=4.0d0 / (1.4217d-3*dble(nt-nit)*dt*pil*e0**2) do 10 ix=1,3 c dipole and rotational strength for propagation along ix: sp(iw,2+ix)=0.5d0*fc*(di(2,ix)+di(1,ix)) 10 sp(iw,5+ix)=4.0d0*fc*(di(2,ix)-di(1,ix)) if(CHIR.eq.2)then do 101 ix=1,3 101 sp(iw,5+ix)=sp(iw,5+ix)/2.0d0 endif c average strengths: sp(iw,1)=(sp(iw,3)+sp(iw,4)+sp(iw,5))/3.0d0 1 sp(iw,2)=(sp(iw,6)+sp(iw,7)+sp(iw,8))/3.0d0 c light frequency w=wmin-dw do 91 iw=1,nw w=w+dw do 91 ix=1,8 91 write(14+ix,600)w*CM,sp(iw,ix) 600 format(f12.2,e15.6) do 9 io=1,nfiles 9 close(10+io) if(lloc)call wrq(nat,ql,q,r) stop end subroutine wrq(nat,ql,q,r) implicit none integer*4 nat,ia,ix,q(*) real*8 ql(*),r(*),wmin,wmax,bohr,qa bohr=0.529177d0 wmin=ql(1) wmax=ql(1) do 1 ia=2,nat if(ql(ia).gt.wmax)wmax=ql(ia) 1 if(ql(ia).lt.wmin)wmin=ql(ia) open(99,file='XL.X') write(99,991)nat 991 format(' relative absorption/scattering in atomic charges',/,i6) qa=0.0d0 do 2 ia=1,nat ql(ia)=(ql(ia)-wmin)/(wmax-wmin)-0.5d0 2 qa=qa+ql(ia) qa=qa/dble(nat) do 3 ia=1,nat ql(ia)=ql(ia)-qa 3 write(99,992)q(ia),(bohr*r(ix+3*(ia-1)),ix=1,3),ql(ia) 992 format(i3,3f12.6,7(' 0'),f10.5) close(99) end subroutine setfield(idd,w,t,e0,E,Ei,B,Bi,Gr,Gi,pol) c E=E0 exp(-iwt), rot E=-dB/dt implicit none integer*4 id,i2,i3,pol,idd real*8 w,t,E(3),B(3),clight,e0,phase,Ei(3),Bi(3),s, 1e0s,e0c,Gr(3,3),Gi(3,3) id=idd if(pol.eq.1)then s=1.0d0 else if(pol.eq.2)then s=-1.0d0 else if(pol.eq.0)then s=0.0d0 else write(6,*)pol call report('Unknown polarization requested') endif endif endif phase=-w*t e0c=e0*cos(phase) e0s=e0*sin(phase) clight=137.036d0 call setij(id,i2,i3) c Electric field E(id)=0.0d0 E(i2)= e0c E(i3)= s*e0s Ei(id)=0.0d0 Ei(i2)= e0s Ei(i3)=-s*e0c c magnetic field B(id)=0.0d0 B(i2)=-E(i3)/clight B(i3)= E(i2)/clight Bi(id)=0.0d0 Bi(i2)=-Ei(i3)/clight Bi(i3)= Ei(i2)/clight c electric field gradient Gr=0.0d0 Gi=0.0d0 Gr(id,i2)=-Ei(i2)*w/clight Gr(id,i3)=-Ei(i3)*w/clight Gi(id,i2)= E( i2)*w/clight Gi(id,i3)= E( i3)*w/clight return end subroutine fieldramp(idd,E,Ei,B,Bi,Gr,Gi,s) c scattered radiation, time-phase from freq. difference w, c CPL c E=E0 exp(-iwt), rot E=-dB/dt implicit none integer*4 id,i2,i3,idd real*8 E(3),B(3),clight,Ei(3),Bi(3),s,k, 1Gr(3,3),Gi(3,3),CM,w532 CM=219474.63d0 w532=1.0d7/532.0d0/CM clight=137.036d0 k=w532/clight id=idd call setij(id,i2,i3) E=0.0d0 Ei=0.0d0 B=0.0d0 Bi=0.0d0 Gr=0.0d0 Gi=0.0d0 c i(kx-wt) c Electric field p e E(i2) =1.0d0 Ei(i3)=s c magnetic field Bi(i2)=-Ei(i3)/clight B(i3 )= E( i2)/clight c electric field gradient Gi(id,i2)= k*E(i2) Gr(id,i3)=-k*Ei(i3) return end subroutine setfieldram(idd,E,Ei,Et,Eti,B,Bi,Bt,Bti,Gr,Gi,pol) c excitation radiation, planar c E=E0 exp(-iwt), rot E=-dB/dt implicit none integer*4 id,i2,i3,pol,idd real*8 E(3),B(3),w532,clight,Ei(3),Bi(3),k, 1Gr(3,3),Gi(3,3),CM,Et(3),Eti(3),Bt(3),Bti(3) CM=219474.63d0 w532=1.0d7/532.0d0/CM clight=137.036d0 k=w532/clight id=idd call setij(id,i2,i3) E =0.0d0 Ei=0.0d0 B =0.0d0 Bi=0.0d0 Et =0.0d0 Eti=0.0d0 Bt =0.0d0 Bti=0.0d0 Gr=0.0d0 Gi=0.0d0 if(pol.eq.1)then c i(-kx-wt) c i2 planar polarization E = p e E(i2)=1.0d0 Eti(i2)=-w532 Gi(id,i2)=-k B(i3)= -1.0d0/clight Bti(i3)=k else if(pol.eq.2)then c i3 planar polarization E(i3)=1.0d0 Eti(i3)=-w532 Gi(id,i3)=-k B(i2)= 1.0d0/clight Bti(i2)=-k else write(6,*)pol call report('Unknown polarization requested') endif endif return end subroutine setdf(idd,tt,t,e0,E,B,Bt,G,pol) c pulse c E=E0 delta(t-t0-r/v) implicit none integer*4 id,i2,i3,pol,idd real*8 t,E(3),B(3),Bt(3),clight,e0,s, 1G(3,3),tt,delta,dd,pi,deltat,t0 t0=3.0d0*tt id=idd if(pol.eq.1)then s=1.0d0 else if(pol.eq.2)then s=-1.0d0 else if(pol.eq.0)then s=0.0d0 else write(6,*)pol call report('Unknown polarization requested') endif endif endif pi=4.0d0*datan(1.0d0) dd=((t-t0)/tt)**2 if(dd.lt.20.0d0)then delta=e0/tt/dsqrt(pi)*exp(-dd) deltat=-2.0d0*(t-t0)/tt**2*delta else delta=0.0d0 deltat=0.0d0 endif clight=137.036d0 call setij(id,i2,i3) c Electric field E=0.0d0 E( i2)= delta E( i3)=-s*delta c magnetic field B=0.0d0 Bt=0.0d0 B( i2)=-s*delta /clight B( i3)= delta /clight Bt(i2)=-s*deltat/clight Bt(i3)= deltat/clight c electric field gradient G=0.0d0 G(id,i2)= -deltat/clight G(id,i3)= -s*deltat/clight return end subroutine setmi(mi,id) implicit none integer*4 id,i2,i3 real*8 mi(3,3),clight clight=137.036d0 c Bext = M . Eext mi=0.0d0 call setij(id,i2,i3) mi(i2,i3)=-1.0d0/clight mi(i3,i2)= 1.0d0/clight return end subroutine calca(ipol,d,w,alpha,xalpha,ar,ai,br,bi,cr,ci, 1rlim,ips,CHIR,ks) c field at each polarized unit implicit none integer*4 ipol,id,ia,ix,jx,d,ip,CHIR,ips,kx real*8 ai(ipol,3,3),ar(ipol,3,3),br(ipol,3,3),bi(ipol,3,3), 1w,xalpha(ipol,3),rlim,alpha(ipol,3,3),cr(ipol,3,3,3), 1ci(ipol,3,3,3),ks, 1a1i(3,3),a1r(3,3),b1i(3,3),b1r(3,3),c1i(3,3,3),c1r(3,3,3) id=d ip=ips do 1 ia=1,ipol call calcpo(xalpha(ia,1),xalpha(ia,2),xalpha(ia,3),ipol, 1id,w,alpha,xalpha,a1r,a1i,b1r,b1i,c1r,c1i,rlim,1,ip,CHIR,ks) do 1 ix=1,3 do 1 jx=1,3 ar(ia,ix,jx)=a1r(ix,jx) ai(ia,ix,jx)=a1i(ix,jx) br(ia,ix,jx)=b1r(ix,jx) bi(ia,ix,jx)=b1i(ix,jx) do 1 kx=1,3 cr(ia,ix,jx,kx)=c1r(ix,jx,kx) 1 ci(ia,ix,jx,kx)=c1i(ix,jx,kx) return end subroutine calcx(ipol,d,w,alpha,xalpha,ar,ai,br,bi,cr,ci,rlim,ips, 1CHIR,po,ks) c fields at each polarized unit are the same, as in po implicit none integer*4 ipol,id,ia,ix,jx,d,ip,CHIR,ips,kx real*8 ai(ipol,3,3),ar(ipol,3,3),br(ipol,3,3),bi(ipol,3,3), 1ci(ipol,3,3,3),cr(ipol,3,3,3),c1r(3,3,3),c1i(3,3,3), 1w,xalpha(ipol,3),rlim,alpha(ipol,3,3),po(3),ks, 1a1i(3,3),a1r(3,3),b1i(3,3),b1r(3,3) id=d ip=ips call calcpo(po(1),po(2),po(3),ipol, 1id,w,alpha,xalpha,a1r,a1i,b1r,b1i,c1r,c1i,rlim,1,ip,CHIR,ks) do 1 ia=1,ipol do 1 ix=1,3 do 1 jx=1,3 ar(ia,ix,jx)=a1r(ix,jx) ai(ia,ix,jx)=a1i(ix,jx) br(ia,ix,jx)=b1r(ix,jx) bi(ia,ix,jx)=b1i(ix,jx) do 1 kx=1,3 cr(ia,ix,jx,kx)=c1r(ix,jx,kx) 1 ci(ia,ix,jx,kx)=c1i(ix,jx,kx) return end subroutine calcpo(x,y,z,ipol,d,w,alpha,xalpha,ar,ai,br,bi, 1cr,ci,rlim,ic,ip,CHIR,ks) c fields at (x y z) c ic=1 adding original field c ip,CHIR - polarization, chirality c ks ... sign of the wave vector implicit none integer*4 ipol,id,I,ix,jx,kx,ic,d,CHIR,ip,ki real*8 tr(3,3),ti(3,3),k,clight,atr(3,3),ati(3,3),s,ks, 1tgr(3,3,3),tgi(3,3,3),cr(3,3,3),ci(3,3,3), 1ai(3,3),ar(3,3),r,r2,r3,r4,r5,r6,r7,br(3,3),bi(3,3),rlim, 1rar(3,3),rai(3,3),w,xalpha(ipol,3),rij(3),xyz(3), 1alpha(ipol,3,3),wr(3,3),wi(3,3),sp,cp,x,y,z,mi(3,3) c c if CHIR=2 and opposite polarization, make for enantiomeric system c R -> -R if(CHIR.eq.2.and.ip.eq.2)then s=-1.0d0 else s= 1.0d0 endif id=d clight=137.036d0 k=w/clight ar=0.0d0 ai=0.0d0 br=0.0d0 bi=0.0d0 cr=0.0d0 ci=0.0d0 xyz(1)=x*s xyz(2)=y*s xyz(3)=z*s if(ic.eq.1)then c adding ext. field tensor with a phase factor exp(ik.ks.r) call setmi(mi,id) cp=cos(k *xyz(id)) sp=sin(k*ks*xyz(id)) do 100 ix=1,3 ar(ix,ix)=cp ai(ix,ix)=sp do 100 jx=1,3 cr(ix,jx,jx)=-k*sp ci(ix,jx,jx)= k*cp br(ix,jx)=mi(ix,jx)*cp 100 bi(ix,jx)=mi(ix,jx)*sp endif do 1 I=1,ipol c positions of polarizable groups with respect to xyz: do 111 ix=1,3 111 rij(ix)=(s*xalpha(I,ix)-xyz(ix)) r2=rij(1)**2+rij(2)**2+rij(3)**2 r=dsqrt(r2) if(r.gt.rlim)then r3=r *r2 r4=r2*r2 r5=r2*r3 r6=r3*r3 r7=r2*r5 c r x alpha: do 11 ix=1,3 rar(1,ix)=rij(2)*alpha(I,3,ix)-rij(3)*alpha(I,2,ix) rar(2,ix)=rij(3)*alpha(I,1,ix)-rij(1)*alpha(I,3,ix) 11 rar(3,ix)=rij(1)*alpha(I,2,ix)-rij(2)*alpha(I,1,ix) c allow real alpha so far: rai=0.0d0 c distance tensor t, its derivative tg, and w tr =0.0d0 ti =0.0d0 tgr=0.0d0 tgi=0.0d0 do 2 ix=1,3 c tij = (k^2r^2-1+ikr)/r^3 delta(i,j) -ri rj (r^2 k^2 - 3 + 3ikr)/r^5 tr(ix,ix)=tr(ix,ix)+(r2*k**2-1.0d0)/r3 ti(ix,ix)=ti(ix,ix)+k/r2 do 21 kx=1,3 tgr(ix,ix,kx)=9.0d0/r5-3.0d0*k**2/r3 21 tgi(ix,ix,kx)=-8.0d0*k/r4 do 2 jx=1,3 tr(ix,jx)=tr(ix,jx)-rij(ix)*rij(jx)*(r2*k**2-3.0d0)/r5 ti(ix,jx)=ti(ix,jx)-rij(ix)*rij(jx)*k*3.0d0/r4 do 22 kx=1,3 tgr(ix,ix,kx)=tgr(ix,jx,kx) 1 -(15.0d0/r7-3.0d0*k**2/r5)*rij(ix)*rij(ix)*rij(jx)*rij(kx) 22 tgi(ix,ix,kx)=tgi(ix,ix,kx) 1 +12.0d0*k/r6*rij(ix)*rij(ix)*rij(jx)*rij(kx) c w=(k/c)(k/r^2+i/r^3)* r x alpha wr(ix,jx)=(k/clight)*(k/r2*rar(ix,jx) - rai(ix,jx)/r3) 2 wi(ix,jx)=(k/clight)*(k/r2*rai(ix,jx) + rar(ix,jx)/r3) c product with polarizability, t.alpha (first index new field): atr=0.0d0 ati=0.0d0 do 3 ix=1,3 do 3 jx=1,3 do 3 kx=1,3 atr(ix,jx)=atr(ix,jx)+tr(kx,ix)*alpha(I,kx,jx) 3 ati(ix,jx)=ati(ix,jx)+ti(kx,ix)*alpha(I,kx,jx) c multiply by phase factor exp[i(k.r_I + k.r_Il)], add cp=cos(k*(ks*s*xalpha(I,id)+r)) sp=sin(k+(ks*s*xalpha(I,id)+r)) do 4 ix=1,3 do 4 jx=1,3 br(ix,jx)=br(ix,jx)+wr(ix,jx)*cp-wi(ix,jx)*sp bi(ix,jx)=bi(ix,jx)+wr(ix,jx)*cp+wr(ix,jx)*sp ar(ix,jx)=ar(ix,jx)+atr(ix,jx)*cp-ati(ix,jx)*sp ai(ix,jx)=ai(ix,jx)+atr(ix,jx)*sp+ati(ix,jx)*cp c c .. first index x-derivative c c .. second index new field c c .. third index old field do 4 kx=1,3 ci(ix,jx,kx)=ci(ix,jx,kx)+k*rij(ix)/r* 1 (atr(jx,kx)*cp-ati(jx,kx)*sp) cr(ix,jx,kx)=cr(ix,jx,kx)-k*rij(ix)/r* 1 (atr(jx,kx)*sp+ati(jx,kx)*cp) do 4 ki=1,3 cr(ix,jx,kx)=cr(ix,jx,kx)+tgr(jx,ki,ix)*alpha(I,ki,kx)*cp 4 ci(ix,jx,kx)=ci(ix,jx,kx)+tgr(jx,ki,ix)*alpha(I,ki,kx)*sp endif 1 continue return end subroutine setiu(ipol,nat,iu,r,xalpha) implicit none integer*4 ipol,nat,iu(*),ia,I,umin real*8 r(*),xalpha(ipol,3),xia,yia,zia,rmin,d do 1 ia=1,nat xia=r(1+3*(ia-1)) yia=r(2+3*(ia-1)) zia=r(3+3*(ia-1)) umin=1 rmin=1.0d90 c what is the closest unit to atom ia? do 2 I=1,ipol d=(xia-xalpha(I,1))**2+(yia-xalpha(I,2))**2+(zia-xalpha(I,3))**2 if(d.lt.rmin)then rmin=d umin=I endif 2 continue 1 iu(ia)=umin return end subroutine redres2(u,ipol,am,bm,E,B,Ea) implicit none integer*4 ipol,u,ix real*8 am(ipol,3,3),bm(ipol,3,3), 1E(3),B(3),Ea(3) c set field for this time and atom/unit do 1 ix=1,3 Ea (ix)=am(u,ix,1)*E(1)+am(u,ix,2)*E(2)+am(u,ix,3)*E(3) 1 B (ix)=bm(u,ix,1)*E(1)+bm(u,ix,2)*E(2)+bm(u,ix,3)*E(3) return end subroutine redress(u,ipol,ar,ai,br,bi,Er,Ei,E,Eii,B,Bii) implicit none integer*4 ipol,u,ix real*8 ar(ipol,3,3),ai(ipol,3,3),br(ipol,3,3),bi(ipol,3,3), 1E(3),Ei(3),B(3),Er(3),Bii(3),Eii(3) c set field for this time and atom/unit do 1 ix=1,3 c E=a.Ei, real part: E (ix)=ar(u,ix,1)*Er(1)-ai(u,ix,1)*Ei(1) 1 +ar(u,ix,2)*Er(2)-ai(u,ix,2)*Ei(2) 1 +ar(u,ix,3)*Er(3)-ai(u,ix,3)*Ei(3) Eii(ix)=ai(u,ix,1)*Er(1)+ar(u,ix,1)*Ei(1) 1 +ai(u,ix,2)*Er(2)+ar(u,ix,2)*Ei(2) 1 +ai(u,ix,3)*Er(3)+ar(u,ix,3)*Ei(3) c B=b.Ei, real and imaginary parts: B (ix)=br(u,ix,1)*Er(1)-bi(u,ix,1)*Ei(1) 1 +br(u,ix,2)*Er(2)-bi(u,ix,2)*Ei(2) 1 +br(u,ix,3)*Er(3)-bi(u,ix,3)*Ei(3) 1 Bii(ix)=bi(u,ix,1)*Er(1)+br(u,ix,1)*Ei(1) 1 +bi(u,ix,2)*Er(2)+br(u,ix,2)*Ei(2) 1 +bi(u,ix,3)*Er(3)+br(u,ix,3)*Ei(3) return end subroutine wrf(nat,ipol,pr,nc,cx,q,r,xalpha,wmin,alpha,dw,rlim) implicit none integer*4 nat,q(*),nc(3),i,ix,iy,iz,id,pr,jx,ipol real*8 r(*),cx(12),xalpha(ipol,3),wmin,w,alpha(ipol,3,3), 1ar(3,3),ai(3,3),bmr(3,3),bmi(3,3),Er(3),Ei(3),Br(3),Bi(3), 1tr(3),ti(3),dw,x,y,z,rlim,Gr(3,3),Gi(3,3),cr(3,3,3),ci(3,3,3) real*8,allocatable::ii(:) character*8 cn(3) data cn/'intx.cub','inty.cub','intz.cub'/ c write(6,*)'Making intensity cube' if(ipol.eq.0)then write(6,*)'ipol = 0, nothing to do' return endif w=wmin+dw*dble(pr-1) allocate(ii(nc(3))) do 3 id=1,3 call setfield(id,w,0.0d0,1.0d0,Er,Ei,Br,Bi,Gr,Gi,1) open(3,file=cn(id)) write(3,3000)nat+ipol,(cx(ix),ix=1,3),nc(1),(cx(i),i=4,6), 1nc(2),(cx(i),i=7,9),nc(3),(cx(i),i=10,12) 3000 format(' Nuclear density',/,' SCF Total Density',/, 13(i5,3f12.6,/),i5,3f12.6) do 2 i=1,nat 2 write(3,3001)q(i),dble(q(i)),(r(ix+3*(i-1)),ix=1,3) 3001 format(i5,4f12.6) do 4 i=1,ipol 4 write(3,3001)107,dble(107),(xalpha(i,ix),ix=1,3) do 1 ix=1,nc(1) do 1 iy=1,nc(2) do 11 iz=1,nc(3) x=cx(1)+dble(ix-1)*cx(4)+dble(iy-1)*cx(7)+dble(iz-1)*cx(10) y=cx(2)+dble(ix-1)*cx(5)+dble(iy-1)*cx(8)+dble(iz-1)*cx(11) z=cx(3)+dble(ix-1)*cx(6)+dble(iy-1)*cx(9)+dble(iz-1)*cx(12) call calcpo(x,y,z,ipol,id,w,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1rlim,0,1,0,1.0d0) do 12 jx=1,3 tr(jx)=ar(jx,1)*Er(1)-ai(jx,1)*Ei(1) 1 +ar(jx,2)*Er(2)-ai(jx,2)*Ei(2) 1 +ar(jx,3)*Er(3)-ai(jx,3)*Ei(3) 12 ti(jx)=ar(jx,1)*Ei(1)+ai(jx,1)*Er(1) 1 +ar(jx,2)*Ei(2)+ai(jx,2)*Er(2) 1 +ar(jx,3)*Ei(3)+ai(jx,3)*Er(3) 11 ii(iz)=(tr(1)*tr(1)+tr(2)*tr(2)+tr(3)*tr(3) 1 +ti(1)*ti(1)+ti(2)*ti(2)+ti(3)*ti(3))/1.0d0 1 write(3,4000)(ii(iz),iz=1,nc(3)) 4000 format(6E13.5) close(3) 3 write(6,*)cn(id) return end subroutine wpo(nat,ipol,xalpha,alpha,qq,P,A,pv,rlim) c read F.INP and calculate field at point=po, and also the spectrum implicit none integer*4 nat,i,ii,ix,iy,iz,id,jx,ipol,nfiles,iw,nm,ip,n,L real*8 xalpha(ipol,3),w,alpha(ipol,3,3),B(3),Bii(3),rlim, 1ar(3,3),ai(3,3),bmr(3,3),bmi(3,3),Er(3),Ei(3),Br(3),Bi(3), 1tr(3),ti(3),z0,m(3),u(3),WR,FD,debye,CM,AMU,clight,po,pi, 1qq(3*nat,6),QLD(6),Q(3,3),BOHR,auk,aulambda,P(3*nat,3),sp, 1A(3*nat,3),FC,rs,ds,dsi(3),rsi(3),Vr(2),Vi(2),wau,pv(3), 1Gr(3,3),Gi(3,3),Egr(3,3),Egi(3,3),Vdr(2),Vdi(2),Vmr(2),Vmi(2), 1cr(3,3,3),ci(3,3,3) real*8,allocatable::s(:,:),ew(:) parameter (nfiles=8) character*9 cn(nfiles) data cn/'DOG00.TAB','DOG0X.TAB','DOG0Y.TAB','DOG0Z.TAB', 1 'DOGH0.TAB','DOGHX.TAB','DOGHY.TAB','DOGHZ.TAB'/ z0=0.0d0 debye=2.541732d0 CM=219474.63d0 AMU=1822.887D0 clight=137.03604d0 BOHR=0.52917715d0 FD =debye*dsqrt(CM/AMU/2.0d0) FC=debye*dsqrt(2.0d0/CM/AMU)/clight pi=4.0d0*atan(1.0d0) po=dsqrt(2.0d0*pi) c do 4 i=1,nfiles open(9+i,file=cn(i)) 4 WRITE (9+i,501) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2/,'-----------------------------------------------------------') n=3*nat allocate(s(n,n),ew(n)) call reads(s,n,nm,ew) c loop over transition frequencies: do 1 iw=1,nm c wavenumber /cm-1: w=ew(iw) wau=w/CM wr=dsqrt(dabs(w)) c wavelength in atomic units: aulambda=1.0d0/w*1.0d8/BOHR c wave vector in atomic units auk=2.0d0*pi/aulambda DO 57 ix=1,6 QLD(ix)=0.0d0 DO 57 L=1,n 57 QLD(ix) = QLD(ix) + qq(L,ix)*s(L,iw) Q(1,1)=QLD(1) Q(2,2)=QLD(2) Q(3,3)=QLD(3) Q(1,2)=QLD(4) Q(2,1)=QLD(4) Q(1,3)=QLD(5) Q(3,1)=QLD(5) Q(2,3)=QLD(6) Q(3,2)=QLD(6) DO 59 ix=1,3 DO 59 jx=1,3 59 Q(ix,jx)=Q(ix,jx)*FD/WR*0.5d0 c transition electric and magnetic dipoles in debyes: do 3 ix=1,3 u(ix)=z0 m(ix)=z0 do 31 ii=1,n u(ix)=u(ix)+s(ii,iw)*P(ii,ix) 31 m(ix)=m(ix)+s(ii,iw)*A(ii,ix) if(dabs(w).gt.1.0d0)u(ix)=u(ix)/WR*FD 3 m(ix)=m(ix)*WR*FC ds=u(1)*u(1)+u(2)*u(2)+u(3)*u(3) rs=u(1)*m(1)+u(2)*m(2)+u(3)*m(3) WRITE (10,523)iw,w,ds,rs,rs,po*w*ds 523 FORMAT(I5,f16.7,4G15.6) do 58 ix=1,3 iy=ix+1 if(iy.eq.4)iy=1 iz=iy+1 if(iz.eq.4)iz=1 ds= 1.5d0*(u(iy)*u(iy) +u(iz)*u(iz) ) rs=-auk*1.5d0*(u(iy)*Q(iz,ix)-u(iz)*Q(ix,iy)) 1+ 1.5d0*(u(iy)*m(iy) +u(iz)*m(iz) ) 58 WRITE (10+ix,523)iw,w,ds,rs,rs,po*w*ds c our definition of magnetic dipole: do 33 ix=1,3 33 m(ix)=m(ix)*clight do 12 id=1,3 c calculate the dressing matrices at point pv: call calcpo(pv(1),pv(2),pv(3),ipol,id,wau,alpha,xalpha,ar,ai, 1bmr,bmi,cr,ci,rlim,1,1,0,1.0d0) do 11 ip=1,2 c calculate external field in t=0 for direction id polarization ip call setfield(id,wau,z0,1.0d0,Er,Ei,Br,Bi,Gr,Gi,ip) do 2 ix=1,3 c electric field EF = A . E0: tr(ix)= ar( ix,1)*Er(1)-ai( ix,1)*Ei(1) 1 +ar( ix,2)*Er(2)-ai( ix,2)*Ei(2) 1 +ar( ix,3)*Er(3)-ai( ix,3)*Ei(3) ti(ix)= ar( ix,1)*Ei(1)+ai( ix,1)*Er(1) 1 +ar( ix,2)*Ei(2)+ai( ix,2)*Er(2) 1 +ar( ix,3)*Ei(3)+ai( ix,3)*Er(3) c magnetic field BF = B . E0: B(ix) =bmr(ix,1)*Er(1)-bmi(ix,1)*Ei(1) 1 +bmr(ix,2)*Er(2)-bmi(ix,2)*Ei(2) 1 +bmr(ix,3)*Er(3)-bmi(ix,3)*Ei(3) Bii(ix)=bmi(ix,1)*Er(1)+bmr(ix,1)*Ei(1) 1 +bmi(ix,2)*Er(2)+bmr(ix,2)*Ei(2) 1 +bmi(ix,3)*Er(3)+bmr(ix,3)*Ei(3) c electric field gradient do 2 jx=1,3 Egr(ix,jx)=Gr(ix,jx) 2 Egi(ix,jx)=Gi(ix,jx) c if(id.eq.1.and.iw.eq.1)call wre(w,ip,id,Er,Ei,Br,Bi,tr,ti,B,Bii, 11.0d0) c <0|V|1> element, V=-uE-mB-QGE: Vdr(ip)=-sp(u,tr) Vdi(ip)=-sp(u,ti) Vmr(ip)= sp(m,Bii) Vmi(ip)=-sp(m,B) if(id.eq.1.and.iw.eq.1)then write(6,609)Vdr(ip),Vmr(ip),Vdi(ip),Vmi(ip) 609 format(' Vdr:',E17.8,/,' Vmr:',E17.8,/, 1 ' Vdi:',E17.8,/,' Vmi:',E17.8,/) endif Vr(ip)=Vdr(ip)+Vmr(ip) Vi(ip)=Vdi(ip)+Vmi(ip) do 13 ix=1,3 do 13 jx=1,3 Vr(ip)=Vr(ip)-Q(ix,jx)*Egr(ix,jx) 13 Vi(ip)=Vi(ip)-Q(ix,jx)*Egi(ix,jx) 11 continue c ip if(id.eq.1.and.iw.eq.1)then write(6,*)'Differences:' write(6,609)Vdr(1)**2-Vdr(2)**2,Vmr(1)**2-Vmr(2)**2, 1 Vdi(1)**2-Vdi(2)**2,Vmi(1)**2-Vmi(2)**2 endif dsi(id)=(Vr(2)**2+Vi(2)**2+Vr(1)**2+Vi(1)**2)/4.0d0*3.0d0 rsi(id)=(Vr(2)**2+Vi(2)**2-Vr(1)**2-Vi(1)**2)/8.0d0*3.0d0 12 WRITE (14+id,523)iw,w,dsi(id),rsi(id),rsi(id),po*w*dsi(id) ds=(dsi(1)+dsi(2)+dsi(3))/3.0d0 rs=(rsi(1)+rsi(2)+rsi(3))/3.0d0 1 WRITE (14,523)iw,w,ds,rs,rs,po*w*ds do 5 i=1,nfiles WRITE(9+i,1819) 1819 FORMAT(60('-')) 5 close(9+i) stop end subroutine reads(s,n,nm,e) implicit none integer*4 n,nm,nat,i,ia,k real*8 s(n,n),e(n) open(9,file='F.INP',status='old') read(9,*)nm,nat,nat do 8 i=1,nat 8 read(9,*) read(9,*) do 3 ia=1,nat do 3 i=1,nm 3 read(9,*)(s(3*(ia-1)+k,i),k=1,2),(s(3*(ia-1)+k,i),k=1,3) read(9,*) read(9,*)(e(i),i=1,nm) write(6,*)'F.INP read' close(9) return end subroutine wrp(nat,ipol,pr,nc,cpx,q,r,xalpha,wmin,alpha,dw,rlim) implicit none integer*4 nat,q(*),nc(2),i,ix,iy,id,pr,jx,ipol real*8 r(*),cpx(9),xalpha(ipol,3),wmin,w,alpha(ipol,3,3),rlim, 1ar(3,3),ai(3,3),bmr(3,3),bmi(3,3),Er(3),Ei(3),Br(3),Bi(3), 1tr(3),ti(3),dw,x,y,z,c3(3),rx,ry,c10(3),c2p(3),bohr,ii,an, 1Gr(3,3),Gi(3,3),cr(3,3,3),ci(3,3,3) character*8 cn(3) data cn/'intx.txt','inty.txt','intz.txt'/ c write(6,*)'Making intensity plane' if(ipol.eq.0)then write(6,*)'ipol = 0, nothing to do' return endif w=wmin+dw*dble(pr-1) bohr=0.529177d0 c origin: cpx(1),cpx(2),cpx(3) c vectors defining plane: c1=cpx(4),cpx(5),cpx(6) c c2=cpx(7),cpx(8),cpx(9) c c3 = c1 x c2 c c2p = c3 x c1, normalized c plane x .. along c1 c plane y .. along c2p c3(1)=cpx(5)*cpx(9)-cpx(6)*cpx(8) c3(2)=cpx(6)*cpx(7)-cpx(4)*cpx(9) c3(3)=cpx(4)*cpx(8)-cpx(5)*cpx(7) an=1.0d0/dsqrt(c3(1)**2+c3(2)**2+c3(3)**2) c3(1)=c3(1)*an c3(2)=c3(2)*an c3(3)=c3(3)*an an=1.0d0/dsqrt(cpx(4)**2+cpx(5)**2+cpx(6)**2) c10(1)=cpx(4)*an c10(2)=cpx(5)*an c10(3)=cpx(6)*an c2p(1)=c3(2)*c10(3)-c3(3)*c10(2) c2p(2)=c3(3)*c10(1)-c3(1)*c10(3) c2p(3)=c3(1)*c10(2)-c3(2)*c10(1) open(9,file='lsta.txt') do 2 i=1,nat c projection into the plane x=r(1+3*(i-1))-cpx(1) y=r(2+3*(i-1))-cpx(2) z=r(3+3*(i-1))-cpx(3) rx=x*c10(1)+y*c10(2)+z*c10(3) ry=x*c2p(1)+y*c2p(2)+z*c2p(3) 2 write(9,3001)q(i),rx*bohr,ry*bohr do 4 i=1,ipol x=xalpha(i,1)-cpx(1) y=xalpha(i,2)-cpx(2) z=xalpha(i,3)-cpx(3) rx=x*c10(1)+y*c10(2)+z*c10(3) ry=x*c2p(1)+y*c2p(2)+z*c2p(3) 4 write(9,3001)107,rx*bohr,ry*bohr 3001 format(i5,2f12.6) close(9) do 3 id=1,3 call setfield(id,w,0.0d0,1.0d0,Er,Ei,Br,Bi,Gr,Gi,1) open(3,file=cn(id)) do 1 ix=1,nc(1) do 1 iy=1,nc(2) x=dble(ix-1)*cpx(4)+dble(iy-1)*cpx(7) y=dble(ix-1)*cpx(5)+dble(iy-1)*cpx(8) z=dble(ix-1)*cpx(6)+dble(iy-1)*cpx(9) rx=x*c10(1)+y*c10(2)+z*c10(3) ry=x*c2p(1)+y*c2p(2)+z*c2p(3) x=cpx(1)+x y=cpx(2)+y z=cpx(3)+z call calcpo(x,y,z,ipol,id,w,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1rlim,0,1,0,1.0d0) do 12 jx=1,3 tr(jx)=ar(jx,1)*Er(1)-ai(jx,1)*Ei(1) 1 +ar(jx,2)*Er(2)-ai(jx,2)*Ei(2) 1 +ar(jx,3)*Er(3)-ai(jx,3)*Ei(3) 12 ti(jx)=ar(jx,1)*Ei(1)+ai(jx,1)*Er(1) 1 +ar(jx,2)*Ei(2)+ai(jx,2)*Er(2) 1 +ar(jx,3)*Ei(3)+ai(jx,3)*Er(3) ii=tr(1)*tr(1)+tr(2)*tr(2)+tr(3)*tr(3) 1 +ti(1)*ti(1)+ti(2)*ti(2)+ti(3)*ti(3) 1 write(3,4000)rx*bohr,ry*bohr,ii/1.0d0 4000 format(6E13.5) close(3) 3 write(6,*)cn(id) return end subroutine aps(n,d0,fme,lwt) implicit none integer*4 n,lwt,io real*8 d0(n),fme(*),ss,nd,ff ss=0.0d0 nd=0.0d0 do 38 io=1,n nd=nd+d0(io)**2 38 ss=ss+d0(io)*fme(io) if(nd.gt.1.0d-10.and.ss.gt.0.0d0)then ff=0.0d0 if(lwt.eq.1)ff=1.0d0 if(lwt.eq.2)ff=2.0d0 if(lwt.eq.3)ff=2.0d0*(1.0d0-exp(-nd/0.0025d0)) if(lwt.lt.4)then ff=ff*ss/nd do 40 io=1,n 40 fme(io)=fme(io)-d0(io)*ff endif if(lwt.eq.4)then do 41 io=1,n 41 fme(io)=-fme(io) endif if(lwt.eq.5)then ff=2.0d0*exp(-nd/0.0025d0) - 1.0d0 do 42 io=1,n 42 fme(io)=fme(io)*ff endif endif return end subroutine reqd(fn,q,nat) implicit none character*(*) fn integer*4 i,ia,ix,nat real*8 q(3*nat,6) open(9,file=fn) do 1 i=1,9 1 read(9,*) do 2 ia=1,nat do 2 ix=1,3 2 read(9,93)(q(ix+3*(ia-1),i),i=1,6) 93 format(9x,3e13.4) close(9) write(6,*)fn//' read' return end subroutine READQD(n,qq,r,P,quad) implicit none integer*4 n,L,J,i real*8 qq(n,6),r(*),x,y,z,P(n,3) logical lex,quad if(quad)then inquire(file='FILE.QEN',exist=lex) if(lex)then call reqd('FILE.QEN',qq,n/3) else DO 2204 L=1,n/3 DO 2204 J=1,3 i=J+3*(L-1) x=r(1+3*(L-1)) y=r(2+3*(L-1)) z=r(3+3*(L-1)) qq(i,1)=P(i,1)*x+P(i,1)*x qq(i,2)=P(i,2)*y+P(i,2)*y qq(i,3)=P(i,3)*z+P(i,3)*z qq(i,4)=P(i,1)*y+P(i,2)*x qq(i,5)=P(i,1)*z+P(i,3)*x 2204 qq(i,6)=P(i,2)*z+P(i,3)*y write(6,*)' Quadrupole derivs from APT' endif else write(6,*)' Quadrupole derivs not used' qq=0.0d0 endif return end subroutine wrc(n,i,ix,io,iw,ip,id,it,d0,dp, 1f,fe,fme,kk,lgeo,nat,q,r) implicit none integer*4 n,i,ix,io,iw,ip,id,it,nat,ia,jx,q(*) real*8 d0(n),dp,f,fe,kk,r(*), 1bohr,fme logical lgeo bohr=0.529177d0 write(6,776)i,ix,io,iw,ip,id,it 776 format(' Coordinate too big',/, 1' Atom:',i4,' xyz: ',i1,' (',i4,'), iw ip id it:',i6,2i4,i9) write(6,777)d0(io),d0(io),dp, 1f,fe,fme, 1-kk*d0(io) 777 format(' dM d0 dP:',3e14.5,/, 1' F :',e14.5,/, 1' electrical :',e14.5,/, 1' mechanical :',e14.5,/, 1' arbitrary :',e14.5,/) if(lgeo)then open(40,file='CC.X') write(40,400)nat 400 format(' Distorted geometry',/,i6) do 1 ia=1,nat 1 write(40,401)q(ia),(d0(jx+3*(ia-1))+r(jx+3*(ia-1))*bohr,jx=1,3) 401 format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0') close(40) write(6,*)'CC.X' endif stop end subroutine readopt(e0,gamma,nt,wmin,wmax,nw,pr,nr,nit,dt,ipa, 1ida,ftresh,nproc,kk,lgeo,lw0,ipol,lcage,cx,nc,lplane,cpx,ncp, 1lvcd,lram,lwt,l3,l4,lijk,rlim,lpoint,po,CHIR,quad,fcoord,fft, 1sft,tt,usea,useg,cram,temp,lloc,lrt,rtc) implicit none integer*4 nt,nw,pr,nit,nr,ipa,ida,nproc,iar,xr,ipol,ix,nc(3), 1ncp(2),lwt,omp_get_num_procs,CHIR real*8 CM,e0,gamma,wmin,wmax,dt,fsau,ftresh,kk,bohr,cx(12), 1cpx(9),rlim,po(3),tt,temp,rtc character*5 key logical lgeo,lw0,lcage,lplane,lram,l3,l4,lijk,lpoint,quad, 1fcoord,fft,sft,lvcd,usea,useg,cram,lloc,lrt character*1 xyz(3) data xyz/'x','y','z'/ bohr=0.529177d0 CM=219474.63d0 fsau=41.34137d0 temp=293.15d0 rtc=0.1d0 lrt=.true. cram=.false. lloc=.false. usea=.true. useg=.true. tt=10.0d0 cx=0.0d0 CHIR=0 dt=0.5d0 e0=1.0d-3 fcoord=.false. fft=.false. sft=.false. lvcd=.true. ftresh=1.0d-4 gamma=1.0d-4 ipa=1 ida=1 ipol=0 kk=1.0d-3 lgeo=.false. lw0=.false. lwt=1 lram=.false. lpoint=.false. l3=.false. l4=.false. lijk=.false. lcage=.false. lplane=.false. nproc=1 nr=1 nc=0 nt=10000 nit=5000 nw=100 po=0.0d0 pr=1 quad=.true. rlim=5.0d0 wmax=1850 wmin=1450 open(9,file='HOE.OPT') 1 read(9,90,end=99,err=99)key 90 format(a5) c polarization for which path is written if(key(1:4).eq.'APOL')read(9,*)ipa c plot the intensity in a defined cage if(key(1:4).eq.'CAGE')then read(9,*)lcage read(9,*)(cx(ix),ix=1,3) read(9,*)nc(1),(cx(ix),ix=4,6) read(9,*)nc(2),(cx(ix),ix=7,9) read(9,*)nc(3),(cx(ix),ix=10,12) do 2 ix=1,12 2 cx(ix)=cx(ix)/bohr endif c force complex Raman path: if(key(1:4).eq.'CRAM')read(9,*)cram c type of accounting for chirlaity if(key(1:4).eq.'CHIR')read(9,*)CHIR c rotation and translation coupling for the correction: if(key(1:3).eq.'CRT')read(9,*)crt c delta function width if(key(1:4).eq.'DELT')read(9,*)tt c electric field,au if(key(1:2).eq.'E0')read(9,*)e0 c force coordinate to POINT if(key.eq.'FCOOR')read(9,*)fcoord c Fast Fourier transform: if(key(1:3).eq.'FFT')read(9,*)fft c damping constant, au if(key.eq.'GAMMA')read(9,*)gamma c write geometries if(key(1:3).eq.'GEO')read(9,*)lgeo c number of polarizabilities (if <>0, POL.LST expected if(key(1:4).eq.'IPOL')read(9,*)ipol c arbitrary constant to keep atoms in place if(key(1:2).eq.'KK')read(9,*)kk c use all cubic and all iijk quartic constants if(key(1:4).eq.'LIJK')read(9,*)lijk c cubic constants: if(key(1:2).eq.'L3')read(9,*)l3 c quartic constants: if(key(1:2).eq.'L4')read(9,*)l4 c write atomic energies for mode localization: if(key(1:3).eq.'LOC')read(9,*)lloc c correct rotation and translation: if(key(1:3).eq.'LRT')read(9,*)lrt c avoid negative slips if(key(1:3).eq.'LW0')read(9,*)lw0 c avoid negative slip facon: if(key(1:3).eq.'LWT')read(9,*)lwt c step when the integration starts if(key(1:4).eq.'NINT')read(9,*)nit c number of processors for parallel run: if(key.eq.'NPROC')read(9,*)nproc c coordinate to record if(key(1:2).eq.'NR')read(9,*)nr c number of point if(key(1:2).eq.'NT')read(9,*)nt c number of light frequencies if(key(1:2).eq.'NW')read(9,*)nw c path to record if(key(1:4).eq.'PATH')read(9,*)pr c write the field in a plane if(key(1:5).eq.'PLANE')then read(9,*)lplane read(9,*)(cpx(ix),ix=1,3) read(9,*)ncp(1),(cpx(ix),ix=4,6) read(9,*)ncp(2),(cpx(ix),ix=7,9) do 3 ix=1,9 3 cpx(ix)=cpx(ix)/bohr endif c use field at a point in ana analytical formula: if(key.eq.'POINT')then read(9,*)lpoint read(9,*)po(1),po(2),po(3) endif c use the quadrupole derivatives if(key(1:4).eq.'QUAD')read(9,*)quad c do Raman spectrum if(key(1:5).eq.'RAMAN')read(9,*)lram c direction for which path is written if(key(1:4).eq.'RDIR')read(9,*)ida c distance limit for polarization interactions: if(key(1:4).eq.'RLIM')read(9,*)rlim c time step, fs: if(key(1:4).eq.'STEP')read(9,*)dt c Slow Fourier transform: if(key(1:3).eq.'SFT')read(9,*)sft c Temperature for Raman, in K: if(key(1:4).eq.'TEMP')read(9,*)temp c treshold for the force field if(key.eq.'THRES')read(9,*)ftresh c use A for ROA: if(key(1:4).eq.'USEA')read(9,*)usea c use G for ROA: if(key(1:4).eq.'USEG')read(9,*)useg c VCD,IR: if(key(1:3).eq.'VCD')read(9,*)lvcd c minimal and maximal light frequencies if(key.eq.'WLMIN'.or.key(1:4).eq.'WMIN')read(9,*)wmin if(key.eq.'WLMAX'.or.key(1:4).eq.'WMAX')read(9,*)wmax goto 1 99 close(9) iar=(nr+2)/3 xr=nr-3*(iar-1) write(6,62)e0,gamma,nt, wmin,wmax,nw, dt,nit, rlim,CHIR,quad, 1ftresh,kk,crt,lrt,nproc, 1kk,lgeo,lw0,ipol, lvcd,lram,cram, usea,useg,lloc, 1fcoord,fft,sft, tt, l3,l4,lijk,lcage,lplane,lwt, 1pr,iar,xyz(xr),nr,ipa,xyz(ida),lpoint,po 62 format(' E0 :',e12.4,', Gamma:',e12.4,', NT:',i9,/, 1 ' Wmin:',f12.2,' cm-1, Wmax:',f12.2,' cm-1, Nw:',i9,/, 1 ' dt:',f12.2,' fs ',12x ,' Nint:',i9,/, 1 ' rlim:',f12.2,' A CHIR:',i12 ,' Quadrup.:',l9,/, 1 'Ftres:',e12.3,' kk:',e12.4,' Nproc:',i9,/, 1 ' KK:',e12.3,' CRT:',e12.4,' LRT:',i9,/, 1 'Geoms:',l12 ,' LW0:',l12 ,' ipol:',i9,/, 1 ' VCD:',l12 ,' RAMAN:',l12 ,' CRAM:',l9,/, 1 ' useA:',l12 ,' useG:',l12 ,' LOC:',l9,/, 1 'Fcoor:',l12 ,' FFT:',l12 ,' SFT:',l9,/, 1 'tt/FT:',f12.2,' fs',/, 1 ' L3:',l12 ,' L4:',l12 ,' LIJK:',l9,/, 1 ' Cage:',l12 ,' Plane:',l12 ,' LWT:',i9,/,/, 1 ' Record: PATH(iw)',i4,', Atom',i4,a1,' (NR',i4, 1 '), APOL',i2,', RDIR',1x,a1,/,/, 1 ' Point:',l2,3f12.4,' A',/) c all in atomic units: wmin=wmin/CM wmax=wmax/CM tt=tt*fsau dt=dt*fsau rlim=rlim/bohr do 4 ix=1,3 4 po(ix)=po(ix)/bohr if(nproc.gt.omp_get_num_procs())then write(6,*)nproc,omp_get_num_procs() call report('not enough processors') endif return end SUBROUTINE READFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) OPEN(20,FILE='FILE.FC',STATUS='OLD') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' CLOSE(20) RETURN END SUBROUTINE READF34(N,F3,fn,s) IMPLICIT none integer*4 N,n3,i,j,k,id real*8 F3(N,N,N),av,avd character*(*) fn,s F3=0.0d0 n3=0 av=0.0d0 avd=0.0d0 id=0 OPEN(20,FILE=fn,STATUS='OLD') 1 read(20,*,end=33,err=33)i,j,k,F3(i,j,k) if(i.eq.j.and.j.eq.k)then id=id+1 avd=avd+dabs(F3(i,j,k)) endif av=av+dabs(F3(i,j,k)) n3=n3+1 goto 1 33 close(20) write(6,*)n3,' ',s//' constants read' if(n3.gt.0)write(6,600)av/dble(n3) 600 format(' Average absolute value: ',E12.4) if(id.gt.0)write(6,601)avd/dble(id) 601 format(' Average diagonal absolute value:',E12.4) return end SUBROUTINE READTEN(NAT,Pout,Aout,LAPT,X) IMPLICIT none INTEGER*4 NAT,L,I,J,ID,IG REAL*8 Pout(NAT*3,3),Aout(NAT*3,3),X(*),EPS,SU real*8,allocatable::AI(:,:,:),AJ(:,:,:),PV(:,:,:),R(:,:), 1A(:,:,:),P(:,:,:) LOGICAL LAPT OPEN(15,FILE='FILE.TEN') allocate(AI(NAT,3,3),AJ(NAT,3,3),PV(NAT,3,3),R(3,NAT), 1A(NAT,3,3),P(NAT,3,3)) DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I+3*(L-1)) READ (15,*) DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SU=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SU=SU+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)=SU WRITE(6,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (AI(L,J,I),I=1,3) DO 230 L=1,NAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) WRITE(6,*)'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) DO 5 L=1,NAT DO 5 I=1,3 DO 5 J=1,3 Pout(I+3*(L-1),J)=P(L,I,J) 5 Aout(I+3*(L-1),J)=A(L,I,J) C CLOSE(15) RETURN END FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA DIMENSION VCD(NAT,3,3),VEL(NAT,3,3),DIP(NAT,3,3),C(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 subroutine READM(nat,m) implicit none integer*4 i,nat real*8 m(nat),AMU,mw AMU=1822.0d0 open(9,file='FTRY.INP', status='old') read(9,*) read(9,*) read(9,90)(m(i),i=1,nat) 90 format(6F12.6) read(9,*) read(9,90)(m(i),i=1,nat) close(9) mw=0.0d0 do 1 i=1,nat mw=mw+m(i) 1 m(i)=m(i)*AMU WRITE(*,6)mw 6 format(' Mw = ',f12.2,' g/mol') return end subroutine setij(id,i2,i3) implicit none integer*4 id,i2,i3 i2=id+1 if(i2.gt.3)i2=1 i3=i2+1 if(i3.gt.3)i3=1 return end subroutine COUNTFF(n,ff,ftresh,nbb) implicit none integer*4 n,nbb,i,j real*8 ff(n,n),ftresh,dd,df nbb=0 dd=0.0d0 df=0.0d0 do 1 i=1,n dd=dd+dabs(ff(i,i)) nbb=nbb+1 do 2 j=1,i-1 2 df=df+dabs(ff(i,j)) do 1 j=1,n 1 if(i.ne.j.and.dabs(ff(i,j)).ge.ftresh)nbb=nbb+1 dd=dd/dble(n) df=df/dble(n*(n-1)) write(6,700)nbb,n**2,ftresh,dd,df 700 format(i6,' FF elements of',i12,' bigger than',e10.3,/, 1' Average | diagonal element|:',e10.3,/, 1' Average |off-diagonal element|:',e10.3) return end subroutine BUFFF(n,nbb,ns,ne,jb,ff,fb,ftresh) implicit none integer*4 n,ns(n),ne(n),nbb,ibb,i,j,jb(nbb) real*8 ff(n,n),ftresh,fb(nbb) ibb=0 do 1 i=1,n ibb=ibb+1 ns(i)=ibb fb(ibb)=ff(i,i) jb(ibb)=i do 11 j=1,n if(j.ne.i.and.dabs(ff(i,j)).ge.ftresh)then ibb=ibb+1 fb(ibb)=ff(i,j) jb(ibb)=j endif 11 continue ne(i)=ibb 1 continue return end subroutine COUNTF3(n,c,tr,lijk,nijj,niij,nijk) implicit none integer*4 n,i,j,k,nijj,niij,nijk,niii real*8 c(n,n,n),tr logical lijk niii=0 nijj=0 niij=0 nijk=0 do 1 i=1,n niii=niii+1 do 1 j=1,n if(j.ne.i)then if(dabs(c(i,j,j)).ge.tr)nijj=nijj+1 if(dabs(c(i,i,j)).ge.tr)niij=niij+1 if(lijk)then do 11 k=1,n 11 if(j.ne.k.and.i.ne.k.and. 1 dabs(c(i,j,k)).ge.tr)nijk=nijk+1 endif endif 1 continue write(6,701)'iii',niii,n write(6,701)'ijj',nijj,n*(n-1) write(6,701)'iij',niij,n*(n-1) write(6,701)'ijk',nijk,n*(n-1)*(n-2) 701 format(1x,a3,':',i6,' elements of',i12,' kept') write(6,700)tr 700 format(' Threshold:',e10.3,/) return end subroutine BUFF3(lijk,n,n_ijj,n_iij,n_ijk,j_iij,j_ijj, 1j_ijk,k_ijk,c,tr,ciii,cijjb,ciijb,cijkb) implicit none integer*4 n,n_ijj(*),n_iij(*),n_ijk(*),j_iij(*), 1j_ijj(*),j_ijk(*),i,j,k,k_ijk(*) real*8 c(n,n,n),tr,ciii(*),ciijb(*),cijjb(*),cijkb(*) logical lijk do 1 i=1,n ciii(i)=c(i,i,i) n_ijj(i)=0 n_iij(i)=0 n_ijk(i)=0 do 1 j=1,n if(j.ne.i)then if(dabs(c(i,j,j)).ge.tr)then n_ijj(i)=n_ijj(i)+1 cijjb(n_ijj(i))=c(i,j,j) j_ijj(n_ijj(i))=j endif if(dabs(c(i,i,j)).ge.tr)then n_iij(i)=n_iij(i)+1 ciijb(n_iij(i))=c(i,i,j) j_iij(n_iij(i))=j endif if(lijk)then do 12 k=1,n if(j.ne.k.and.i.ne.k)then if(dabs(c(i,j,k)).ge.tr)then n_ijk(i)=n_ijk(i)+1 cijkb(n_ijk(i))=c(i,j,k) j_ijk(n_ijk(i))=j k_ijk(n_ijk(i))=k endif endif 12 continue endif endif 1 continue return end subroutine COUNTF4(n,c,tr,lijk,niijj,njjij,niiij,niijk,njjik) implicit none integer*4 n,i,j,k,niijj,njjij,niijk,niiii,niiij,njjik real*8 c(n,n,n),tr logical lijk niiii=0 niijj=0 niiij=0 njjij=0 niijk=0 njjik=0 do 1 i=1,n niiii=niiii+1 do 1 j=1,n if(j.ne.i)then if(dabs(c(i,j,j)).ge.tr)niijj=niijj+1 if(dabs(c(j,i,j)).ge.tr)njjij=njjij+1 if(dabs(c(i,i,j)).ge.tr)niiij=niiij+1 if(lijk)then do 11 k=1,n if(j.ne.k.and.i.ne.k)then if(dabs(c(i,j,k)).ge.tr)niijk=niijk+1 if(dabs(c(j,i,k)).ge.tr)njjik=njjik+1 endif 11 continue endif endif 1 continue write(6,701)'iiii',niiii,n write(6,701)'iijj',niijj,n*(n-1) write(6,701)'iiij',niiij,n*(n-1) write(6,701)'jjij',njjij,n*(n-1) write(6,701)'iijk',niijk,n*(n-1)*(n-2) write(6,701)'jjik',njjik,n*(n-1)*(n-2) 701 format(1x,a4,':',i6,' elements of',i12,' kept') write(6,700)tr 700 format(' Threshold:',e10.3,/) return end subroutine add3(n,d0,fme, 1ciii, cijjb,ciijb,cijkb, 1 j_ijj,j_iij,j_ijk,k_ijk, 1 n_ijj,n_iij, n_ijk,lijk) implicit none integer*4 n,n_ijj(*),n_iij(*),j_iij(*),n_ijk(*),j_ijk(*),k_ijk(*), 1io,jj,j_ijj(*) real*8 d0(n),fme(n),ciii(n),cijjb(*),ciijb(*),cijkb(*),xo,ft logical lijk do 363 io=1,n ft=0.0d0 xo=d0(io) ft=ft-0.5d0*ciii(io)*xo**2 do 3632 jj=1,n_iij(io) 3632 ft=ft-ciijb(jj)*d0(j_iij(jj))*xo do 3631 jj=1,n_ijj(io) 3631 ft=ft-0.5d0*cijjb(jj)*d0(j_ijj(jj))**2 if(lijk)then do 3633 jj=1,n_ijk(io) 3633 ft=ft-0.5d0*cijkb(jj)*d0(j_ijk(jj))*d0(k_ijk(jj)) endif 363 fme(io)=fme(io)+ft return end subroutine BUFF4(lijk,n,n_iijj,n_iiij,n_jjij,n_iijk,n_jjik,j_iijj, 1 j_iiij,j_jjij,j_iijk,k_iijk,j_jjik,k_jjik, 1 d,tr,diiii,diijjb,djjijb,diiijb,diijkb,djjikb) implicit none integer*4 n,n_iijj(*),n_iiij(*),n_jjij(*),n_iijk(*),j_iijj(*), 1j_iiij(*),j_jjij(*),j_iijk(*),k_iijk(*),i,j,k,n_jjik(*), 1j_jjik(*),k_jjik(*) real*8 d(n,n,n),tr,diiii(*),diijjb(*),diijkb(*), 1djjikb(*),djjijb(*),diiijb(*) logical lijk do 1 i=1,n diiii(i)=d(i,i,i) n_iijj(i)=0 n_iiij(i)=0 n_jjij(i)=0 n_iijk(i)=0 n_jjik(i)=0 do 1 j=1,n if(j.ne.i)then if(dabs(d(i,j,j)).ge.tr)then n_iijj(i)=n_iijj(i)+1 diijjb(n_iijj(i))=d(i,j,j) j_iijj(n_iijj(i))=j endif if(dabs(d(i,i,j)).ge.tr)then n_iiij(i)=n_iiij(i)+1 diiijb(n_iiij(i))=d(i,i,j) j_iiij(n_iiij(i))=j endif if(dabs(d(j,i,j)).ge.tr)then n_jjij(i)=n_jjij(i)+1 djjijb(n_jjij(i))=d(j,i,j) j_jjij(n_jjij(i))=j endif if(lijk)then do 12 k=1,n if(j.ne.k.and.i.ne.k)then if(dabs(d(i,j,k)).ge.tr)then n_iijk(i)=n_iijk(i)+1 diijkb(n_iijk(i))=d(i,j,k) j_iijk(n_iijk(i))=j k_iijk(n_iijk(i))=k endif if(dabs(d(j,i,k)).ge.tr)then n_jjik(i)=n_jjik(i)+1 djjikb(n_jjik(i))=d(j,i,k) j_jjik(n_jjik(i))=j k_jjik(n_jjik(i))=k endif endif 12 continue endif endif 1 continue return end subroutine add4(lijk,n,d0,fme,n_iijj,n_jjij,n_iiij,n_iijk, 1j_iijj,j_iiij,j_jjij,j_iijk,k_iijk,n_jjik,j_jjik,k_jjik, 1diiii,diijjb,djjijb,diiijb,diijkb,djjikb) implicit none integer*4 n,io,n_iijj(*),n_jjij(*),n_iiij(*),n_iijk(*), 1j_iijj(*),j_iiij(*),j_jjij(*),j_iijk(*),k_iijk(*),jo,ko, 1n_jjik(*),j_jjik(*),k_jjik(*),jj real*8 xo,d0(n),fme(n),ft,diiii(*),diijjb(*),djjijb(*), 1diiijb(*),diijkb(*),djjikb(*) logical lijk do 1 io=1,n ft=0.0d0 xo=d0(io) ft=ft-diiii(io)*xo**3/6.0d0 do 4631 jj=1,n_iijj(io) 4631 ft=ft-diijjb(jj)*d0(j_iijj(jj))**2*xo do 4632 jj=1,n_jjij(io) 4632 ft=ft-djjijb(jj)*d0(j_jjij(jj))**3/6.0d0 do 4634 jj=1,n_iiij(io) 4634 ft=ft-0.5d0*diiijb(jj)*d0(j_iiij(jj))*xo**2 if(lijk)then do 3 jj=1,n_iijk(io) jo=j_iijk(jj) ko=k_iijk(jj) 3 ft=ft-0.5d0*diijkb(jj)*xo*d0(jo)*d0(ko) do 4 jj=1,n_jjik(io) jo=j_jjik(jj) ko=k_jjik(jj) 4 ft=ft-0.5d0*djjikb(jj)*d0(jo)**2*d0(ko) endif 1 fme(io)=fme(io)+ft return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine ldal(alpha,xalpha,ipol) implicit none integer*4 ipol,a,b,i real*8 alpha(ipol,3,3),xalpha(ipol,3) open(9,file='POL.LST') do 1 i=1,ipol read(9,*)((alpha(i,a,b),a=1,b),b=1,3),(xalpha(i,a),a=1,3) alpha(i,2,1)=alpha(i,1,2) alpha(i,3,1)=alpha(i,1,3) 1 alpha(i,3,2)=alpha(i,2,3) return end subroutine wre(w,ip,id,Er,Ei,Br,Bi,E,Eii,B,Bii,an) implicit none integer*4 ip,id,i real*8 w,Er(3),Ei(3),Br(3),Bi(3),E(3),Eii(3),B(3),Bii(3),EN,BN, 1an write(6,601)w,ip,id 601 format(' w:',f12.2' cm-1 pol:',i2,', dir:',i2,':') write(6,'(a)')' External field:' EN=dsqrt(Er(1)**2+Ei(1)**2+Er(2)**2+Ei(2)**2+Er(3)**2+Ei(3)**2) BN=dsqrt(Br(1)**2+Bi(1)**2+Br(2)**2+Bi(2)**2+Br(3)**2+Bi(3)**2) write(6,607)(Er(i)/an,i=1,3),(Ei(i)/an,i=1,3), 1 (Br(i)/an,i=1,3),(Bi(i)/an,i=1,3),EN/an,BN/an 607 format(' E:',3E11.3,' + i',3E11.3,/,' B:',3E11.3,' + i',3E11.3,/, 1' |E|:',f15.6,' |B|:',f15.6) write(6,'(a)')' Re-dressed:' EN=dsqrt(E(1)**2+Eii(1)**2+E(2)**2+Eii(2)**2+E(3)**2+Eii(3)**2) BN=dsqrt(B(1)**2+Bii(1)**2+B(2)**2+Bii(2)**2+B(3)**2+Bii(3)**2) write(6,607)(E(i)/an,i=1,3),(Eii(i)/an,i=1,3), 1 (B(i)/an,i=1,3),(Bii(i)/an,i=1,3),EN/an,BN/an return end function sp(a,b) implicit none real*8 a(3),b(3),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine ft(nt0,nw,ipol,CHIR,nat,n,po,alpha,xalpha,ar,ai,bmr, 1bmi,rlim,tt,fb,ns,ne,m,lw0,lwt,kk,P,A,gamma,dt,lgeo,q,ipa,nr, 1ida,fcoord,qq,e0,iu,jb,r,wmin,dw) implicit none integer*4 nt,nt0,ip,id,ipol,CHIR,n,it,ipu,io,jj,ns(*),ne(*),lwt, 1i,nat,q(*),ipa,nr,ida,nfiles,iu(*),iw,jb(*),ix,nw real*8 rlim,po(3),t,fb(*),m(*),ma,kk,gamma,gdt,r(*),t3,dw, 1alpha(ipol,3,3),xalpha(ipol,3),E(3),B(3),Bt(3), 1dt,G(3,3),pi,CM,qq(n,6),dp,re,im,w,fc,e0, 1f,P(n,3),A(n,3),tt,dt2,ftn,sp(8),c,s,cp,ss,s0,c0, 1ar(ipol,3,3),ai(ipol,3,3),bmr(ipol,3,3),bmi(ipol,3,3),df,cf,z2, 1wmin,Ea(3),cr(ipol,3,3,3),ci(ipol,3,3,3) real*8,allocatable::d0(:),d1(:),dttr(:),di(:,:,:),fme(:),dtti(:) logical lw0,lgeo,fcoord parameter (nfiles=12) character*7 fn(nfiles) data fn/'PNR.PRN','ANR.PRN','ATL.PRN', 1 'ADI.PRN','ABS.PRN','VCD.PRN', 1 'ABX.PRN','ABY.PRN','ABZ.PRN', 1 'VCX.PRN','VCY.PRN','VCZ.PRN'/ do 8 io=1,nfiles 8 open(10+io,file=fn(io)) write(6,*)'sft started' z2=2.0d0 pi=4.0d0*datan(1.0d0) gdt=gamma*dt cf=1.0d0+gamma*dt/2.0d0 df=1.0d0-gamma*dt/2.0d0 dt2=dt**2 CM =219474.63d0 nt=nt0 allocate(d0(n),d1(n),dttr(nw),dtti(nw),fme(n),di(nw,2,3)) di=0.0d0 C$OMP Parallel Do Schedule(Dynamic) Default(Shared) C$OMP+ Private(ip,id,ar,ai,bmr,bmi,cr,ci,d0,d1,dttr,dtti,t,it,ipu,ss, C$OMP+ E,B,Bt,G,fme,ftn,i,ma,Ea,ix,io,f,t3,dp,w,iw,re,im,c,s,c0,s0,cp) do 2 ip=1,2 c loop over x y z propagation directions do 25 id=1,3 write(6,602)ip,id 602 format(2i2,$) c set relation between group and external field: c E(l,t) = a Eext(0,t), B(l,t) = bm Eext(0,t) if(ipol.ne.0)then if(fcoord)then c use field at point po: call calcx(ipol,id,0.0d0,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,po,1.0d0) else c use field at closets unit to each atom: call calca(ipol,id,0.0d0,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,1.0d0) endif endif c initialize deviations: d0=0.0d0 d1=0.0d0 dttr=0.0d0 dtti=0.0d0 c loop over time steps: t=-dt do 3 it=1,nt t=t+dt ipu=1 if(CHIR.eq.1)ipu=ip c light at t and r=0: call setdf(id,tt,t,e0,E,B,Bt,G,ipu) c mechanical force= - FF .d: fme=0.0d0 do 36 io=1,n ftn=0.0d0 do 361 jj=ns(io),ne(io) 361 ftn=ftn-fb(jj)*d0(jb(jj)) 36 fme(io)=fme(io)+ftn if(lw0)call aps(n,d0,fme,lwt) do 39 i=1,nat ma=m(i) if(ipol.ne.0)then call redres2(iu(i),ipol,ar,bmr,E,B,Ea) else Ea=E endif do 39 ix=1,3 io=ix+3*(i-1) f=fme(io)-kk*d0(io)+P(io,1)*Ea(1)+P(io,2)*Ea(2)+P(io,3)*Ea(3) t3=0.0d0 if((CHIR.eq.0.and.ip.eq.1).or.CHIR.eq.1) 1t3=-A(io,1)*Bt(1)+A(io,2)*Bt(2)+A(io,3)*Bt(3) 1 + (G(1,1)*qq(io,1)+G(1,2)*qq(io,4)+G(1,3)*qq(io,5) 1 + G(2,1)*qq(io,4)+G(2,2)*qq(io,2)+G(2,3)*qq(io,6) 1 + G(3,1)*qq(io,5)+G(3,2)*qq(io,6)+G(3,3)*qq(io,3))/6.0d0 c XX YY ZZ XY XZ YZ c 1 2 3 4 5 6 if(CHIR.eq.2)then t3=-A(io,1)*Bt(1)+A(io,2)*Bt(2)+A(io,3)*Bt(3) 1 + (G(1,1)*qq(io,1)+G(1,2)*qq(io,4)+G(1,3)*qq(io,5) 1 + G(2,1)*qq(io,4)+G(2,2)*qq(io,2)+G(2,3)*qq(io,6) 1 + G(3,1)*qq(io,5)+G(3,2)*qq(io,6)+G(3,3)*qq(io,3))/6.0d0 if(ip.eq.2)t3=-t3 endif f=f+t3 dp=(z2*d0(io)-d1(io)*df+f*dt2/ma)/cf if(dabs(dp).gt.0.5d0)call wrc(n,i,ix,io,0,ip,id,it, 1d0,dp,f,P(io,1)*Ea(1)+P(io,2)*Ea(2)+P(io,3)*Ea(3), 1fme(io),kk,lgeo,nat,q,r) if(ip.eq.ipa.and.io.eq.nr.and.id.eq.ida)write(11,100)it,d0(io) 100 format(i7,e14.5) c slow ft: s=dp*sin((wmin-dw)*t) c=dp*cos((wmin-dw)*t) s0=sin(dw*t) c0=cos(dw*t) do 777 iw=1,nw ss=c*s0+s*c0 cp=c*c0-s*s0 dttr(iw)=dttr(iw)+cp dtti(iw)=dtti(iw)+ss c=cp 777 s=ss d1(io)=d0(io) 39 d0(io)=dp 3 continue c t w=wmin-dw do 22 iw=1,nw w=w+dw re=dttr(iw) im=dtti(iw) 22 di(iw,ip,id)=di(iw,ip,id)+w**2*gdt*ma*(re**2+im**2) 25 continue c id, directions write(6,*) 2 continue c ip,with and without A w=wmin-dw do 911 iw=1,nw w=w+dw fc=4.0d0 / (1.4217d-3*dble(nt)*dt*pi*e0**2) do 10 ix=1,3 c dipole and rotational strength for propagation along ix: sp(2+ix)=0.5d0*fc*(di(iw,2,ix)+di(iw,1,ix)) 10 sp(5+ix)=4.0d0*fc*(di(iw,2,ix)-di(iw,1,ix)) if(CHIR.eq.2)then do 101 ix=1,3 101 sp(5+ix)=sp(5+ix)/2.0d0 endif sp(1)=(sp(3)+sp(4)+sp(5))/3.0d0 sp(2)=(sp(6)+sp(7)+sp(8))/3.0d0 do 91 ix=1,8 91 write(14+ix,600)w*CM,sp(ix) 600 format(f12.2,e15.6) 911 continue do 9 io=1,nfiles 9 close(10+io) stop end subroutine runft(nt0,ipol,CHIR,nat,n,po,alpha,xalpha,ar,ai,bmr, 1bmi,rlim,tt,fb,ns,ne,m,lw0,lwt,kk,P,A,gamma,dt,lgeo,q,ipa,nr, 1ida,fcoord,qq,e0,iu,jb,r,wmin,wmax) implicit none integer*4 nt,nt0,ip,id,ipol,CHIR,n,it,ipu,io,jj,ns(*),ne(*),lwt, 1i,nat,q(*),ipa,nr,ida,nfiles,iu(*),iw,jb(*),ix real*8 rlim,po(3),t,fb(*),m(*),ma,kk,gamma,gdt,r(*),t3, 1alpha(ipol,3,3),xalpha(ipol,3),E(3),B(3),Bt(3), 1dt,G(3,3),pi,CM,qq(n,6),dp,re,im,w,fc,e0, 1f,P(n,3),A(n,3),tt,dt2,ftn,sp(8), 1ar(ipol,3,3),ai(ipol,3,3),bmr(ipol,3,3),bmi(ipol,3,3),df,cf,z2, 1wmin,wmax,Ea(3),cr(ipol,3,3,3),ci(ipol,3,3,3) real*8,allocatable::d0(:),d1(:),dtt(:),dtr(:,:),di(:,:,:),fme(:) logical lw0,lgeo,fcoord parameter (nfiles=12) character*7 fn(nfiles) data fn/'PNR.PRN','ANR.PRN','ATL.PRN', 1 'ADI.PRN','ABS.PRN','VCD.PRN', 1 'ABX.PRN','ABY.PRN','ABZ.PRN', 1 'VCX.PRN','VCY.PRN','VCZ.PRN'/ do 8 io=1,nfiles 8 open(10+io,file=fn(io)) z2=2.0d0 pi=4.0d0*datan(1.0d0) gdt=gamma*dt cf=1.0d0+gamma*dt/2.0d0 df=1.0d0-gamma*dt/2.0d0 dt2=dt**2 CM =219474.63d0 nt=2**int(log(dble(nt0))/log(2.0d0)) write(6,605)nt 605 format(' Number of time steps changed to',i12) allocate(di(nt,2,3)) di=0.0d0 allocate(d0(n),d1(n),dtr(n,2*nt),dtt(2*nt),fme(n)) write(6,*)'Dynamics' C$OMP Parallel Do Schedule(Dynamic) Default(Shared) C$OMP+ Private(ip,id,ar,ai,bmr,bmi,cr,ci,d0,d1,t,it,ipu,E,B,Bt,G, C$OMP+ fme,ftn,i,ma,Ea,ix,io,f,t3,dp,dtr,dtt,iw,w,re,im) do 2 ip=1,2 c loop over x y z propagation directions do 25 id=1,3 write(6,602)ip,id 602 format(2i2,$) c set relation between group and external field: c E(l,t) = a Eext(0,t), B(l,t) = bm Eext(0,t) if(ipol.ne.0)then if(fcoord)then c use field at point po: call calcx(ipol,id,0.0d0,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,po,1.0d0) else c use field at closets unit to each atom: call calca(ipol,id,0.0d0,alpha,xalpha,ar,ai,bmr,bmi,cr,ci, 1 rlim,ip,CHIR,1.0d0) endif endif c initialize deviations: d0=0.0d0 d1=0.0d0 c loop over time steps: t=-dt do 3 it=1,nt t=t+dt ipu=1 if(CHIR.eq.1)ipu=ip c light at t and r=0: call setdf(id,tt,t,e0,E,B,Bt,G,ipu) c mechanical force= - FF .d: fme=0.0d0 do 36 io=1,n ftn=0.0d0 do 361 jj=ns(io),ne(io) 361 ftn=ftn-fb(jj)*d0(jb(jj)) 36 fme(io)=fme(io)+ftn if(lw0)call aps(n,d0,fme,lwt) do 39 i=1,nat ma=m(i) if(ipol.ne.0)then call redres2(iu(i),ipol,ar,bmr,E,B,Ea) else Ea=E endif do 39 ix=1,3 io=ix+3*(i-1) f=fme(io)-kk*d0(io)+P(io,1)*Ea(1)+P(io,2)*Ea(2)+P(io,3)*Ea(3) t3=0.0d0 if((CHIR.eq.0.and.ip.eq.1).or.CHIR.eq.1) 1t3=-A(io,1)*Bt(1)+A(io,2)*Bt(2)+A(io,3)*Bt(3) 1 + (G(1,1)*qq(io,1)+G(1,2)*qq(io,4)+G(1,3)*qq(io,5) 1 + G(2,1)*qq(io,4)+G(2,2)*qq(io,2)+G(2,3)*qq(io,6) 1 + G(3,1)*qq(io,5)+G(3,2)*qq(io,6)+G(3,3)*qq(io,3))/6.0d0 c XX YY ZZ XY XZ YZ c 1 2 3 4 5 6 if(CHIR.eq.2)then t3=-A(io,1)*Bt(1)+A(io,2)*Bt(2)+A(io,3)*Bt(3) 1 + (G(1,1)*qq(io,1)+G(1,2)*qq(io,4)+G(1,3)*qq(io,5) 1 + G(2,1)*qq(io,4)+G(2,2)*qq(io,2)+G(2,3)*qq(io,6) 1 + G(3,1)*qq(io,5)+G(3,2)*qq(io,6)+G(3,3)*qq(io,3))/6.0d0 if(ip.eq.2)t3=-t3 endif f=f+t3 dp=(z2*d0(io)-d1(io)*df+f*dt2/ma)/cf if(dabs(dp).gt.0.5d0)call wrc(n,i,ix,io,0,ip,id,it, 1d0,dp,f,P(io,1)*Ea(1)+P(io,2)*Ea(2)+P(io,3)*Ea(3), 1fme(io),kk,lgeo,nat,q,r) if(ip.eq.ipa.and.io.eq.nr.and.id.eq.ida)write(11,100)it,d0(io) 100 format(i7,e14.5) c record path for FT: dtr(io,2*it-1)=dp dtr(io,2*it )=0.0d0 d1(io)=d0(io) 39 d0(io)=dp 3 continue c t do 1 io=1,n do 11 it=1,2*nt 11 dtt(it)=dtr(io,it) call four1(dtt,nt,1) do 21 it=1,2*nt 21 dtr(io,it)=dtt(it) 1 continue do 22 iw=1,nt w=z2*pi*dble(iw)/dble(nt)/dt do 22 io=1,n re=dtr(io,2*iw-1) im=dtr(io,2*iw) 22 di(iw,ip,id)=di(iw,ip,id)+w**2*gdt*ma*(re**2+im**2) 25 continue c id, directions write(6,*) 2 continue c ip,with and without A do 911 iw=1,nt w=z2*pi*dble(iw)/dble(nt)/dt if(w.gt.wmin.and.w.lt.wmax)then fc=4.0d0 / (1.4217d-3*dble(nt)*dt*pi*e0**2) do 10 ix=1,3 c dipole and rotational strength for propagation along ix: sp(2+ix)=0.5d0*fc*(di(iw,2,ix)+di(iw,1,ix)) 10 sp(5+ix)=4.0d0*fc*(di(iw,2,ix)-di(iw,1,ix)) if(CHIR.eq.2)then do 101 ix=1,3 101 sp(5+ix)=sp(5+ix)/2.0d0 endif sp(1)=(sp(3)+sp(4)+sp(5))/3.0d0 sp(2)=(sp(6)+sp(7)+sp(8))/3.0d0 do 91 ix=1,8 91 write(14+ix,600)w*CM,sp(ix) 600 format(f12.2,e15.6) endif 911 continue do 9 io=1,nfiles 9 close(10+io) stop end SUBROUTINE four1(data,nn,isign) c Fast Fourier Tranform (FFT) from Numerical recepies, with some c changes for real*8 precision implicit none c nn .. number of complex numbers (2*nn real ones) c data(1) Real f0 c data(2) Im f0 c data(3) Real f1 c data(4) Im f1 etc c return the same ordering. f0 zero-frequency c Replaces data(1:2*nn) by its discrete Fourier transform, if isign is c input as 1; or replaces data(1:2*nn) by nn times its inverse discrete c Fourier transform, if isign is input as -1. c data is a complex array of length nn or, equivalently, a real array c of length 2*nn. nn MUST be an integer power of 2 c (this is not checked for!). INTEGER*4 isign,nn REAL*8 data(2*nn) INTEGER*4 i,istep,j,m,mmax,n REAL*8 tempi,tempr,theta,wi,wpi,wpr,wr,wtemp c n=2*nn j=1 c c bit reversal section: do 11 i=1,n,2 if(j.gt.i)then c exchange two complex numbers: tempr=data(j) tempi=data(j+1) data(j)=data(i) data(j+1)=data(i+1) data(i)=tempr data(i+1)=tempi endif m=nn 1 if ((m.ge.2).and.(j.gt.m)) then j=j-m m=m/2 goto 1 endif 11 j=j+m c c Danielson-Lanczos section: mmax=2 c Outer loop executed log2 nn times. 2 if (n.gt.mmax) then istep=2*mmax c Initialize for the trigonometric recurrence: theta=6.28318530717959d0/(isign*mmax) wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0 wi=0.d0 c Two nested inner loops: do 13 m=1,mmax,2 do 12 i=m,n,istep j=i+mmax c Danielson-Lanczos formula: tempr=wr*data(j)-wi*data(j+1) tempi=wr*data(j+1)+wi*data(j) data(j)=data(i)-tempr data(j+1)=data(i+1)-tempi data(i)=data(i)+tempr 12 data(i+1)=data(i+1)+tempi c Trigonometric recurrence: wtemp=wr wr=wr*wpr-wi*wpi+wr 13 wi=wi*wpr+wtemp*wpi+wi mmax=istep c Not yet done: goto 2 endif c All done: return END SUBROUTINE READTTT(n,ALPHA,A,G) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(n,3,3),A(n,3,3,3),G(n,3,3) open(15,file='FILE.TTT') READ(15,*) READ(15,*)NAT READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) CLOSE(15) RETURN END subroutine setv6(n,v6,r) c set unit trnaslational and rotational vectors implicit none integer*4 n,nat,i1,i2,i3,ia real*8 v6(6,n),r(n),vn nat=n/3 v6=0.0d0 do 1 i1=1,3 call setij(i1,i2,i3) do 1 ia=1,nat v6(i1 ,i1+3*(ia-1))= 1.0d0 v6(i1+3,i2+3*(ia-1))= r(i3+3*(ia-1)) 1 v6(i1+3,i3+3*(ia-1))=-r(i2+3*(ia-1)) do 2 i1=1,6 vn=0.0d0 do 21 i2=1,n 21 vn=vn+v6(i1,i2)**2 vn=1.0d0/dsqrt(vn) do 2 i2=1,n 2 v6(i1,i2)=v6(i1,i2)*vn return end subroutine chkrt(n,d0,v6,crt) c check that molecule is not translated or rotated implicit none integer*4 n,i,i6 real*8 d0(n),v6(6,n),sp(6),crt do 1 i6=1,6 sp(i6)=0.0d0 do 1 i=1,n 1 sp(i6)=sp(i6)+d0(i)*v6(i6,i) do 2 i=1,n do 2 i6=1,6 2 d0(i)=d0(i)-sp(i6)*v6(i6,i)*crt return end