program dynamim implicit none c nw0 ... maximum number of spectral points c maxval maximal valence numbers (number of bonds from one atom) c K0 maximum number of named kinds of molecules c IOX max # of overlaps integer*4 maxval,nw0,nn0,mm,nmol,npair,nnames,K0,IOX,ia,iw,iprint, 1n,nw,ix,NAT1,istep,i,iind,j,jx,k,l,n1,n2,okind,NO,NATS,IWG, 1nsa,ii,mxb,pkind,irep,nm,natx,i1,i2,i3 parameter (nw0=2000,maxval=8,K0=5,IOX=4000) character*10 snam(K0) character*1 lg logical ltti,ltei,lexp,lper,lpair,ldress,lp27,lp28,LWR,LABS,LVCD, 1LROA,LRAM,LSTRICT,LINV,lexca,LAPT,LALF,LRDO,LNAMES,LPOLYMER, 1LHALTONERROR,lcct,lproject,lani,lconv,lp29,lcell,lincr,lnnames, 1lpol,lpolg,lpola,lcidr,lave real*8 dt,fr,spec(nw0),stem(nw0),t,w,dbb, 1ramconst,dw,wmax,tp,sol,wmin,res,sc,fac,dtau,bohr,pi,kelvin, 1vx,vy,vz,cf2,POLM(K0,9),sa1,sal0,sal1,sag0,sag1,ax,ay,az,ac(3), 1bc(3),cc(3),s180,d180,ydy,ydx,tf,rams(4,nw0),perx,pery,perz,px2, 1py2,pz2,po(3,3),gpo(3,3),apo(3,3,3), 1bigpol(3,3),biggpol(3,3),bigapol(3,3,3),POLGM(K0,9),POLAM(K0,27) integer*4,allocatable:: i12(:,:),iii(:,:),atomic(:),isr(:), 1al(:,:),alp(:,:),IIM(:,:,:),natm(:),natmp(:),imk(:), 1ND(:),NDS(:),INDBIG(:,:),IBONDSB(:,:),IND(:), 1IBONDS(:,:),NBO(:),NBOB(:),NATNAME(:),NDNAME(:,:) c wsu 1 el _x sin c wsu 2 el _x cos c wsu 3 el _y sin c wsu 4 el _y cos c wsu 5 el _z sin c wsu 6 el _z cos c wsu 7 mag_x sin c wsu 8 mag_x cos c wsu 9 mag_y sin c wsu10 mag_y cos c wsu11 mag_z sin c wsu12 mag_z cos c for all molecules and spectral points c index = i12 + (iw-1)*12 +(mol-1)*nw0*12 c current (x,y,z), previous (xm,ym,zm) and previous c previous(xmm,ymm,zmm) atomic coordinates c rm are the molecular coordinates real*8,allocatable::x(:),y(:),z(:), 1xm(:),ym(:),zm(:),xmm(:),ymm(:),zmm(:),pchgft(:),pol(:), 2P(:,:,:),A(:,:,:),ALPHAI(:,:,:),AAI(:,:,:,:),GI(:,:,:),coo(:,:), 4ALPHA1(:,:,:),AA1(:,:,:,:),G1(:,:,:),ALPHA(:,:,:),opol(:,:,:), 4ALP29(:,:),ALPHM(:,:),ALPHAM(:,:,:),GM(:,:,:),AAM(:,:,:), 5APTM(:,:,:),AM(:,:,:),AA(:,:,:),G(:,:,:),APTI(:,:,:),AI(:,:,:), 6u(:,:,:),Xwork(:,:),uu(:),gt(:,:),at(:,:),aat(:,:),gopol(:,:,:), 7wsu(:),wsat(:,:),wcat(:,:),wsgt(:,:),wcgt(:,:), 2wsaat(:,:),wcaat(:,:),alave(:),lq(:),gq(:),aq(:),eq(:), 1R(:,:),RS(:,:),PS(:,:,:),AS(:,:,:),ALPHAS(:,:,:),AAS(:,:,:,:), 1GS(:,:,:),AAT_c(:,:,:,:),GT_c(:,:,:),asname(:,:,:,:), 1psname(:,:,:,:),aaname(:,:,:,:,:),gname(:,:,:,:),aopol(:,:,:,:), 1alphaname(:,:,:,:),rname(:,:,:),polname(:,:,:),gpolname(:,:,:), 1apolname(:,:,:,:),lqo(:),gqo(:),aqo(:) common/cctnopts/NO,NATS, 1LWR,LABS,LVCD,IWG,LROA,LRAM,LSTRICT,LINV,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,n1,n2,okind logical,allocatable:: luse(:) write(6,6009) 6009 format(/,/, 1' DYNAMIM: spectra generation from molecular trajectories',/,/, 1' Petr Bour',/, 1' Academy of Sciences, Prague, 2015',/,/, 1' Input: FILE.X MD trajectories',/, 1' <> tensor files',/, 1' DYNAMIM.OPT options',/, 1' Output: .PRN spectra',//) call readopt3(nmol,lpair) nn0=nw0*nmol fr=0.085645d0 nm=0 open(10,file='FILE.X',status='old') read(10,*) read(10,*)n write(6,*)n,' atoms,',nmol,' molecules' allocate(x(n),y(n),z(n),xm(n),ym(n),coo(n,3), 1zm(n),xmm(n),ymm(n),zmm(n),pchgft(n),pol(n),P(n,3,3),A(n,3,3), 2ALPHAI(3*n,3,3),AAI(3*n,3,3,3),GI(3*n,3,3),ALPHA1(3*n,3,3), 3AA1(3*n,3,3,3),G1(3*n,3,3),opol(n,3,3), 4ALP29(9,3*n),ALPHM(9,3*n),ALPHAM(9,K0,3*n),GM(9,K0,3*n), 6AAM(27,K0, 3*n),APTM(n,9,K0),AM(n,9,K0), 7APTI(n,3,3),AI(n,3,3),u(3,3,n),Xwork(3,n), 8i12(maxval,n),iii(4,n),atomic(n),isr(n),gopol(n,3,3), 1IIM(K0,4,n),luse(n),natm(nmol),natmp(nmol),imk(nmol), 6uu(6*nmol),gt(9,nmol),at(9,nmol),aat(27,nmol),aopol(n,3,3,3)) call rg(n,x,y,z,atomic,i12,0,maxval) close(10) if(lpair)then npair=nmol*(nmol-1)/2 write(6,6001)npair 6001 format(i4,' pairs') else npair=nmol endif allocate(al(npair,n),alp(npair,n)) call readopt1(K0,n,NSA,iprint,nmol,nnames,isr,natm,imk,al, 1dt,kelvin,bohr,pi,dbb,perx,pery,perz,tp,sol,dw,wmin,wmax,res, 1ac,bc,cc,lper,lcell,lp27,lp28,lp29,lproject,lnnames,lpair,LABS, 1LVCD,LRAM,LROA,ltti,ldress,ltei,lani,lconv,lg,snam,lincr,lcct, 1lexp,lexca,px2,py2,pz2,sc,dtau,nw,nw0,lcidr,lave,irep,natx) allocate(ALPHA(9,nmol,3*natx),AA(27,nmol,3*natx),G(9,nmol,3*natx)) c do 314 ia=1,n pchgft(ia)=0.0d0 314 luse(ia)=.false. call setactive(NSA,luse,isr,n) c if(lexca)call readcharge(n,pchgft) if(lexp)call readatpolars(n,pol) if(lexp.and.lpair)call fillpairs(nmol,natm,natmp,npair,n,al,alp) c read polarizability derivatives and polarizabilities in c internal coordinates (related to an atom and bonds): if(ltti)call inttt(lnnames,snam,NAT1,n,ALPHAI,AAI,GI,iii, 1K0,IIM,ALPHAM,GM,AAM,po,gpo,apo,POLM,POLGM,POLAM, 1nnames,fr) c read dipole derivatives in internal coordinates: if(ltei)call intte(lnnames,nnames,snam,n,APTI,AI,NAT1,iii, 1K0,IIM,APTM,AM) c writeout some parameters: call ctrl3(iprint,dt,dtau,n,NSA,dw,wmin,wmax,nw,res,lper,perx, 1pery,perz,lnnames,lincr,lp27,lp28,lp29,lcell,LRAM,LROA,LABS,LVCD, 1lani,lconv,ltei,ltti,lcct,lexp,lexca,lproject,ldress,lg) c c predefine convolution function call dc(res,lg,dw,wmin,nw) dw =dw *tp*sol wmax=wmax*tp*sol wmin=wmin*tp*sol c c initiate frequency coefficients: allocate(wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0)) if(lincr)then c restart from previous calculation: call ifcr(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) else c new accumulations: call ifc(wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) istep=0 endif do 30 ia=1,n xm(ia)=x(ia) ym(ia)=y(ia) zm(ia)=z(ia) xmm(ia)=x(ia) ymm(ia)=y(ia) 30 zmm(ia)=z(ia) c input dipole derivatives and polarizabilities in c cartesian coordinates: if(lcct)then mxb=n allocate(R(3,n),ND(n),NDS(n), 1 INDBIG(IOX,n),IBONDSB(8,n),IND(n), 2 IBONDS(8,n),NBO(n),NBOB(n),PS(n,3,3),AS(n,3,3),ALPHAS(3*n,3,3), 3 AAS(3*n,3,3,3),GS(3*n,3,3),AAT_c(3*n,3,3,3),GT_c(3*n,3,3), 1 RS(3,n),polname(IOX,3,3),gpolname(IOX,3,3),apolname(IOX,3,3,3), 1 NATNAME(IOX)) allocate(aaname(IOX,3*natx,3,3,3),gname(IOX,3*natx,3,3), 1 alphaname(IOX,3*natx,3,3),asname(IOX,natx,3,3), 1 rname(IOX,3,natx),psname(IOX,natx,3,3),NDNAME(IOX,natx)) call cctn_tinker2(n,atomic,x,y,z, 1 R,ND,NDS,RS,INDBIG,IBONDSB,IBONDS,NBO,NBOB, 2 PS,AS,ALPHAS,AAS,GS,AAT_c,GT_c, 3 asname,psname,aaname,gname,alphaname,rname,NATNAME, 4 NDNAME,IOX,luse,po,apo,gpo,lpol,lpolg,lpola, 1 polname,gpolname,apolname,fr,natx) endif allocate(lq(9*3*n),gq(9*3*n),aq(27*3*n),eq(3*n)) lq=0 aq=0 gq=0 if(lave)then c average molecular polarizabilities c equilibrium and derivs c A, B, C = X, Y, Z c type index c 1 1 .. 3 a_AA x a_AA c 2 4 .. 6 a_AA x G_AA c 3 7 .. 15 a_AB c 4 16 .. 24 a_AB x a_AB c 5 25 .. 33 a_AA x a_BB c 6 34 .. 42 G_AB c 7 43 .. 51 a_AB x G_AB c 8 52 .. 60 a_AA x G_BB c 9 61 .. 87 a_AA x A_B,AC c 10 88 ..114 a_AA x A_A,BC c 11 115 ..141 a_AB x A_A,AC c Two-molecule averages c 12 142 ..150 ai_AB x aj_AB c 13 151 ..159 2ai_AB x Gj_AB + ai_AA x (Gj_AA - Gj_BB) c 14 160 ..186 2 R_ijC ai_AA x aj_AB c 15 187 ..213 (ai_AA Aj_B,AC + ai_AA Aj_A,BC - 2 aj_AB Ai_A,AC )/3 pkind=213 allocate(alave(pkind)) alave=0.0d0 call inpder(3*n,nm,lq,gq,aq,eq) c restart from previous calculation: if(lincr)call rpca(pkind,0,alave,istep) else nm=1 endif allocate(lqo(nmol*nm*9),gqo(nmol*nm*9),aqo(nmol*nm*27)) lqo=0.0d0 gqo=0.0d0 aqo=0.0d0 c start geometry reading ("dynamic"): c ######################################################## open(10,file='FILE.X') 31111 istep=istep+1 c current time [s] t=dt*dble(istep)*1.0d-12 if(mod(istep,iprint).eq.0)then write(6,6003)istep,int(t*1.0d15) 6003 format(2i20,' fs') endif c read current geometry: read(10,*,end=31112,err=31112) read(10,*,end=31112,err=31112)n call rg(n,x,y,z,atomic,i12,0,maxval) c if(lcct)then if(istep.eq.1)write(6,*)'Tensors by transfer' call trtena(x,y,z,ALPHA1,G1,AA1,P,A,n,mxb,R,IND,RS,INDBIG, 1 IBONDSB,NBOB,PS,AS,ALPHAS,AAS,GS,IOX,luse,asname,psname,aaname, 1 gname,alphaname,rname,NATNAME,po,gpo,apo,bigpol,biggpol,bigapol, 1 lpol,lpolg,lpola,polname,gpolname,apolname,lcidr,lave, 1 nm,lq,gq,aq,nmol,lqo,gqo,aqo,coo,opol,gopol,aopol,natx) c project translations and rotations from alpha if(lproject)call protr(x,y,z,ALPHA1,P,n,istep,luse) if(istep.eq.1)then CALL WRITETTT(n,n,ALPHA1,AA1,G1,'FILE.1.TTT') write(6,6010) 6010 format('First point common origin Tensors to FILE.1.TTT') endif endif c dipole derivatives from internal tensors: if(ltei)call dodd(istep,nmol,npair,n,natm,al,x,y,z,lnnames, 1iii,K0,IIM,imk,u,AM,APTM,AI,APTI,A,P) c c dipole derivatives from charges if(lexca)call doddq(istep,n,pchgft,P,A,x,y,z) c c Alpha from atomic polarizabilities: if(lexp)call doaap(istep,npair,n,natm,al,luse,AA1,G1, 1ALPHA1,x,y,z,Xwork,nmol,AA,G,ALPHA,pol,natx) c c Alpha from internal tensors: if(ltti)call doai(istep,n,nmol,natm,al,x,y,z,lnnames,K0,IIM,imk, 1POLM,POLGM,POLAM,po,gpo,apo,opol,gopol,aopol,ALPHAM,GM,AAM, 1ALPHA1,AA1,G1,Xwork,luse,ALPHA,AA,G,LNAMES,ALPHAI,GI,AAI,npair, 1ALP29,ALPHM,perx,pery,perz,px2,py2,pz2,lp27,lp28,lp29,ldress, 1lq,gq,aq,lqo,gqo,aqo,coo,nm,natx) c numerical averaging of molecular tensor invariants: if(lave)then call avepol(nmol, 1 n,opol,gopol,aopol,alave,coo,perx,pery,perz,lper) if(mod(istep,irep).eq.0) 1 call wrave(pkind,alave,nmol,istep,po,gpo,apo,2) endif c opol(n,3,3) gopol(n,3,3) aopol(n,3,3) molecular polarizabilities c in local coordinates c pkind ... dimension of averaging components c n ... number of atoms c nmol (=NO for cct) .. number of molecules c nm .. number of normal modes c lqo,gqo,aqo normal mode polarizability derivatives, local c coo(n,3) molecular positions if(lcct)then if(istep.eq.1)write(6,*)'Alpha transferred ' DO 201 mm=1,nmol DO 199 ia=1,natm(mm) DO 199 IX=1,3 L=al(mm,ia) i2=3*(ia-1)+IX IIND=3*(L-1)+IX if(luse(L))then DO 200 I=1,3 DO 200 J=1,3 DO 206 K=1,3 206 AA( I+(J-1)*3+(K-1)*9,mm,i2)=AA1( IIND,I,J,K) G (I+(J-1)*3 ,mm,i2)=G1 ( IIND,I,J) 200 ALPHA(I+(J-1)*3, mm,i2)=ALPHA1(IIND,I,J) endif 199 continue 201 continue endif c correct for periodic jumps: if(lper)call cpj(n,lcell,ac,bc,cc,x,y,z,xm,ym,zm,xmm,ymm,zmm, 1perx,pery,perz,px2,py2,pz2) c c update dipoles: if(ltei.or.lexca.or.LABS.or.LVCD)then c dipole time derivatives du/dt=du/dx . du/dt: c constants changing vx,ax into atomic units fac=1.0d0/bohr/dtau cf2=1.0d0/bohr/dtau/dtau do 209 mm=1,nmol do 2010 ia=1,6 2010 uu((mm-1)*6+ia)=0.0d0 do 209 ia=1,natm(mm) L=al(mm,ia) vx=(x(L)-xm(L))*fac vy=(y(L)-ym(L))*fac vz=(z(L)-zm(L))*fac ax=(x(L)+xmm(L)-xm(L)*2.0d0)*cf2 ay=(y(L)+ymm(L)-ym(L)*2.0d0)*cf2 az=(z(L)+zmm(L)-zm(L)*2.0d0)*cf2 c uu-index: 1 - electric moment x c 2 - electric moment y c 3 - electric moment z c 4 - magnetic moment x c 5 - magnetic moment y c 6 - magnetic moment z uu((mm-1)*6+1)=uu((mm-1)*6+1)+vx*P(L,1,1)+vy*P(L,2,1)+vz*P(L,3,1) uu((mm-1)*6+2)=uu((mm-1)*6+2)+vx*P(L,1,2)+vy*P(L,2,2)+vz*P(L,3,2) uu((mm-1)*6+3)=uu((mm-1)*6+3)+vx*P(L,1,3)+vy*P(L,2,3)+vz*P(L,3,3) uu((mm-1)*6+4)=uu((mm-1)*6+4)+ax*A(L,1,1)+ay*A(L,2,1)+az*A(L,3,1) uu((mm-1)*6+5)=uu((mm-1)*6+5)+ax*A(L,1,2)+ay*A(L,2,2)+az*A(L,3,2) 209 uu((mm-1)*6+6)=uu((mm-1)*6+6)+ax*A(L,1,3)+ay*A(L,2,3)+az*A(L,3,3) if(lconv)then call update2c(wmin,dw,nw,wsu,uu,t,sc,nmol) else call update2(wmin,dw,nw,wsu,uu,t,sc,nw0,nmol) endif endif c c update polarizabilities: if(ltti.or.lexp.or.LROA.or.LRAM)then c polarizability time derivatives da/dt=da/dx . dx/dt: fac=1.0d0/dtau/bohr do 2002 mm=1,nmol do 2001 ix=1,9 at(ix,mm)=0.0d0 gt(ix,mm)=0.0d0 aat(ix,mm)=0.0d0 aat(ix+9,mm)=0.0d0 2001 aat(ix+18,mm)=0.0d0 do 2002 ia=1,natm(mm) i1=1+3*(ia-1) i2=2+3*(ia-1) i3=3+3*(ia-1) L=al(mm,ia) vx=(x(L)-xm(L))*fac vy=(y(L)-ym(L))*fac vz=(z(L)-zm(L))*fac do 2002 ii=1,9 gt(ii,mm)=gt(ii,mm) +G( ii,mm,i1)*vx 1 + G( ii,mm,i2)*vy 2 + G( ii,mm,i3)*vz at(ii,mm)=at(ii,mm) +ALPHA(ii,mm,i1)*vx 1 + ALPHA(ii,mm,i2)*vy 2 + ALPHA(ii,mm,i3)*vz aat(ii ,mm)=aat(ii ,mm)+AA( ii ,mm,i1)*vx 1 + AA( ii ,mm,i2)*vy 2 + AA( ii ,mm,i3)*vz aat(ii+ 9,mm)=aat(ii+ 9,mm)+AA( ii+ 9,mm,i1)*vx 1 + AA( ii+ 9,mm,i2)*vy 2 + AA( ii+ 9,mm,i3)*vz 2002 aat(ii+18,mm)=aat(ii+18,mm)+AA( ii+18,mm,i1)*vx 1 + AA( ii+18,mm,i2)*vy 2 + AA( ii+18,mm,i3)*vz c if(lconv)then call update3c(wmin,dw,nw,wsat, 1 wcat,wsgt,wcgt,wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol) else call update3(wmin,dw,nw,wsat, 1 wcat,wsgt,wcgt,wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol,nw0) endif endif c c remember last point do 1 ia=1,n xmm(ia)=xm(ia) ymm(ia)=ym(ia) zmm(ia)=zm(ia) xm(ia)=x(ia) ym(ia)=y(ia) 1 zm(ia)=z(ia) c if(istep.gt.0.and.mod(istep,iprint).eq.0)then c write the normalized power spectra if(ltei.or.lexca.or.LABS.or.LVCD)then c abs/vcd: do 711 iw=1,nw rams(1,iw)=0.0d0 rams(2,iw)=0.0d0 do 711 mm=1,nmol w=wmin+dble(iw-1)*dw fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol 1 /dble(istep) ii=(iw-1)*12+(mm-1)*12*nw0 rams(1,iw)=rams(1,iw)+fac* 1 (wsu(1+ii)**2+wsu(3+ii)**2+wsu(5+ii)**2 2 +wsu(2+ii)**2+wsu(4+ii)**2+wsu(6+ii)**2) 711 rams(2,iw)=rams(2,iw)+fac* 1 (wsu(1+ii)*wsu(8+ii)+wsu(3+ii)*wsu(10+ii)+wsu(5+ii)*wsu(12+ii) 2 -wsu(2+ii)*wsu(7+ii)+wsu(9+ii)*wsu(10+ii)+wsu(6+ii)*wsu(11+ii)) c do 1238 iw=1,nw 1238 spec(iw)=rams(1,iw) if(lconv)then call wrspec(istep,lani,'WD.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WDraw.PRN',nw,wmin,dw,spec,tp,sol) c convolute it with a smoothing function: call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WD.PRN',nw,wmin,dw,spec,tp,sol) endif c do 1239 iw=1,nw 1239 spec(iw)=rams(2,iw) if(lconv)then call wrspec(istep,lani,'WR.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WRraw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WR.PRN',nw,wmin,dw,spec,tp,sol) endif endif c Raman: if(ltti.or.lexp.or.LROA.or.LRAM)then do 1233 iw=1,nw w=wmin+dble(iw-1)*dw rams(1,iw)=0.0d0 rams(2,iw)=0.0d0 rams(3,iw)=0.0d0 rams(4,iw)=0.0d0 do 1233 mm=1,nmol c temperature factor:*1d5(arbitrary factor) if(dabs(w).lt.1.0d-40)then tf=0.0d0 else tf=1.0d0/w/(1.0d0-exp(-w*1.44d0/kelvin))*1.0d5 endif if(lexp)then tf=tf*1.0d-10 else tf=tf*1.0d6 endif sal0=0.0d0 sal1=0.0d0 sag0=0.0d0 sag1=0.0d0 sa1=0.0d0 do 121 jx=1,3 ii=iw+(mm-1)*nw0 sa1=sa1+wsat(1+(jx-1)*3,ii)* 1 (wsaat(2+(3-1)*3+(jx-1)*9,ii)-wsaat(3+(2-1)*3+(jx-1)*9,ii)) 2 +wcat(1+(jx-1)*3,ii)* 2 (wcaat(2+(3-1)*3+(jx-1)*9,ii)-wcaat(3+(2-1)*3+(jx-1)*9,ii)) 3 +wsat(2+(jx-1)*3,ii)* 3 (wsaat(3+(1-1)*3+(jx-1)*9,ii)-wsaat(1+(3-1)*3+(jx-1)*9,ii)) 4 +wcat(2+(jx-1)*3,ii)* 4 (wcaat(3+(1-1)*3+(jx-1)*9,ii)-wcaat(1+(3-1)*3+(jx-1)*9,ii)) 5 +wsat(3+(jx-1)*3,ii)* 5 (wsaat(1+(2-1)*3+(jx-1)*9,ii)-wsaat(2+(1-1)*3+(jx-1)*9,ii)) 6 +wcat(3+(jx-1)*3,ii)* 6 (wcaat(1+(2-1)*3+(jx-1)*9,ii)-wcaat(2+(1-1)*3+(jx-1)*9,ii)) do 121 ix=1,3 sal0=sal0+wsat(ix+(ix-1)*3,ii)*wsat(jx+(jx-1)*3,ii) 1 +wcat(ix+(ix-1)*3,ii)*wcat(jx+(jx-1)*3,ii) sal1=sal1+wsat(ix+(jx-1)*3,ii)*wsat(ix+(jx-1)*3,ii) 1 +wcat(ix+(jx-1)*3,ii)*wcat(ix+(jx-1)*3,ii) c coordinates sag0=sag0+wsat(ix+(ix-1)*3,ii)*wsgt(jx+(jx-1)*3,ii) 1 +wcat(ix+(ix-1)*3,ii)*wcgt(jx+(jx-1)*3,ii) sag1=sag1+wsat(ix+(jx-1)*3,ii)*wsgt(ix+(jx-1)*3,ii) 1 +wcat(ix+(jx-1)*3,ii)*wcgt(ix+(jx-1)*3,ii) 121 continue sa1=sa1/3.0d0 c IL+IR Backscattering: s180=7.0d0*sal1+sal0 c IL+IR Backscattering-xpol: ydy=3.0d0*sal1-sal0 c IL+IR Backscattering-ypol: ydx=4.0d0*sal1+2.0d0*sal0 c IR-IL Backscattering: d180=8.0d0*(3.0d0*sag1-sag0+sa1) ramconst=(w/sol/1000.0d0)**2/dble(istep) rams(1,iw)=rams(1,iw)+s180*tf*ramconst rams(2,iw)=rams(2,iw)+ydy *tf*ramconst rams(3,iw)=rams(3,iw)+ydx *tf*ramconst 1233 rams(4,iw)=rams(4,iw)+d180*tf*ramconst do 1234 iw=1,nw 1234 spec(iw)=rams(1,iw) if(lconv)then call wrspec(istep,lani,'WRam180.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WRam180raw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WRam180.PRN',nw,wmin,dw,spec,tp,sol) endif c do 1235 iw=1,nw c1235 spec(iw)=rams(2,iw) c call wrspec(istep,lani,'WRamyraw.PRN',nw,wmin,dw,spec,tp,sol) c call conv(res,nw,stem,spec,wmin,dw,tp,sol) c call wrspec(istep,lani,'WRamy.PRN',nw,wmin,dw,spec,tp,sol) c do 1236 iw=1,nw c1236 spec(iw)=rams(3,iw) c call wrspec(istep,lani,'WRamxraw.PRN',nw,wmin,dw,spec,tp,sol) c call conv(res,nw,stem,spec,wmin,dw,tp,sol) c call wrspec(istep,lani,'WRamx.PRN',nw,wmin,dw,spec,tp,sol) c do 1237 iw=1,nw c1237 spec(iw)=rams(4,iw) do 1237 iw=1,nw 1237 spec(iw)=rams(4,iw) if(lconv)then call wrspec(istep,lani,'WRoa180.PRN',nw,wmin,dw,spec,tp,sol) else call wrspec(istep,lani,'WRoa180raw.PRN',nw,wmin,dw,spec,tp,sol) call conv(res,nw,stem,spec,wmin,dw,tp,sol) call wrspec(istep,lani,'WRoa180.PRN',nw,wmin,dw,spec,tp,sol) endif endif endif c ######################################################## if(lincr)then if(mod(istep,iprint).eq.0) 1 call ifcw(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) if(lave)call rpca(pkind,1,alave,istep) endif goto 31111 31112 close(10) if(lave)call wrave(pkind,alave,nmol,istep,po,gpo,apo,1) write(6,6002)istep 6002 format(/,i10,' geometries in FILE.X') end subroutine conv(res,nw,stem,spec,wmin,dw,tp,sol) implicit none real*8 res,stem(*),spec(*),wmin,dw,fac,tp,sol,w,wp integer*4 nw,iw,iwp if(res.gt.0)then do 72 iw=1,nw 72 stem(iw)=spec(iw) do 73 iw=1,nw w=(wmin+dble(iw-1)*dw)/tp/sol spec(iw)=0.0d0 do 73 iwp=1,nw wp=(wmin+dble(iwp-1)*dw)/tp/sol fac=((wp-w)/res)**2 73 if(fac.lt.20.0d0)spec(iw)=spec(iw)+stem(iwp)*exp(-fac) endif return end subroutine wrspec(istep,lani,sn,nw,wmin,dw,spec,tp,sol) implicit none character*(*) sn real*8 wmin,dw,spec(*),w,tp,sol integer*4 nw,iw,is,istep logical*4 lani character*20 ts if(lani)then write(ts,90)istep 90 format(i20) do 2 is=1,20 2 if(ts(is:is).eq.' ')ts(is:is)='0' open(8,file=ts(1:20)//sn) else open(8,file=sn) endif do 1 iw=1,nw w=wmin+dble(iw-1)*dw 1 write(8,805)w/tp/sol,spec(iw) 805 format(f12.2,' ',g15.6) close(8) return end subroutine update2(wmin,dw,nw,wsu,u,t,sc,nw0,nmol) implicit none integer*4 nw,iw,nmol,mm,ii,nw0,ix real*8 wmin,w,t,ss,cc,sc,dw,tp,sol,wsu(*),wau,u(*) c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw do 1 iw=1,nw w=w+dw if(dabs(w).gt.1.0d-10)then wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 11 mm=1,nmol ii=(iw-1)*12+(mm-1)*12*nw0 do 11 ix=1,6 wsu( 2*ix-1+ii)=wsu(2*ix-1+ii)*sc+u((mm-1)*6+ix)*ss 11 wsu( 2*ix +ii)=wsu(2*ix +ii)*sc+u((mm-1)*6+ix)*cc endif 1 continue return end subroutine update2c(wmin,dw,nw,wsu,u,t,sc,nmol) c use convolution implicit none integer*4 nw,iw,nmol,mm,ii,nw0,ix,nc,iwp parameter (nw0=2000) real*8 wmin,w,t,ss,cc,sc,dw,tp,sol,wsu(*),wau,u(*),fc,ff common/convc/ff(2*nw0+1),nc c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw do 1 iw=1,nw w=w+dw if(dabs(w).gt.1.0d-10)then wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 11 iwp=max(1,iw-nc+1),min(iw+nc-1,nw) fc=ff(nc+iwp-iw) do 11 mm=1,nmol ii=(iwp-1)*12+(mm-1)*12*nw0 do 11 ix=1,6 wsu(2*ix-1+ii)=wsu(2*ix-1+ii)*sc+u((mm-1)*6+ix)*ss*fc 11 wsu(2*ix +ii)=wsu(2*ix +ii)*sc+u((mm-1)*6+ix)*cc*fc endif 1 continue return end subroutine update3c(wmin,dw,nw,wsat,wcat,wsgt,wcgt, 1wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol) c alpha,A,G -> Raman, ROA c use convolution implicit none integer*4 nw,iw,ix,nn0,nmol,mm,ii,nw0 parameter (nw0=2000) real*8 wmin,w,at(9,nmol),gt(9,nmol),aat(27,nmol), 1t,ss,cc,sc,dw,tp,sol, 2wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0),wcgt(9,nn0), 2wsaat(27,nn0),wcaat(27,nn0),wau integer*4 nc,iwp real*8 fc,ff common/convc/ff(2*nw0+1),nc c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw c do iwp=1,2*nc-1 c write(6,*)iwp,ff(iwp) c enddo c stop do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 1 iwp=max(1,iw-nc+1),min(iw+nc-1,nw) fc=ff(nc+iwp-iw) do 1 mm=1,nmol ii=iwp+(mm-1)*nw0 do 1 ix=1,9 wsat(ix,ii)=wsat(ix,ii)*sc+at(ix,mm)*ss*fc wcat(ix,ii)=wcat(ix,ii)*sc+at(ix,mm)*cc*fc wsgt(ix,ii)=wsgt(ix,ii)*sc+gt(ix,mm)*ss*fc wcgt(ix,ii)=wcgt(ix,ii)*sc+gt(ix,mm)*cc*fc wsaat(ix,ii)=wsaat(ix,ii)*sc+aat(ix,mm)*ss*fc wcaat(ix,ii)=wcaat(ix,ii)*sc+aat(ix,mm)*cc*fc wsaat(ix+9,ii)=wsaat(ix+9,ii)*sc+aat(ix+9,mm)*ss*fc wcaat(ix+9,ii)=wcaat(ix+9,ii)*sc+aat(ix+9,mm)*cc*fc wsaat(ix+18,ii)=wsaat(ix+18,ii)*sc+aat(ix+18,mm)*ss*fc 1 wcaat(ix+18,ii)=wcaat(ix+18,ii)*sc+aat(ix+18,mm)*cc*fc return end subroutine update3(wmin,dw,nw,wsat,wcat,wsgt,wcgt, 1wsaat,wcaat,at,gt,aat,t,sc,nn0,nmol,nw0) c alpha,A,G -> Raman, ROA implicit none integer*4 nw,iw,ix,nn0,nmol,mm,ii,nw0 real*8 wmin,w,at(9,nmol),gt(9,nmol),aat(27,nmol), 1t,ss,cc,sc,dw,tp,sol, 2wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0),wcgt(9,nn0), 2wsaat(27,nn0),wcaat(27,nn0),wau c tp=2.0d0*4.0d0*datan(1.0d0) sol=3.0d10 w=wmin-dw do 1 iw=1,nw w=w+dw wau=w/tp/sol/219470.0d0 ss=-cos(w*t)/wau cc= sin(w*t)/wau do 1 mm=1,nmol ii=iw+(mm-1)*nw0 do 1 ix=1,9 wsat(ix,ii)=wsat(ix,ii)*sc+at(ix,mm)*ss wcat(ix,ii)=wcat(ix,ii)*sc+at(ix,mm)*cc wsgt(ix,ii)=wsgt(ix,ii)*sc+gt(ix,mm)*ss wcgt(ix,ii)=wcgt(ix,ii)*sc+gt(ix,mm)*cc wsaat(ix,ii)=wsaat(ix,ii)*sc+aat(ix,mm)*ss wcaat(ix,ii)=wcaat(ix,ii)*sc+aat(ix,mm)*cc wsaat(ix+9,ii)=wsaat(ix+9,ii)*sc+aat(ix+9,mm)*ss wcaat(ix+9,ii)=wcaat(ix+9,ii)*sc+aat(ix+9,mm)*cc wsaat(ix+18,ii)=wsaat(ix+18,ii)*sc+aat(ix+18,mm)*ss 1 wcaat(ix+18,ii)=wcaat(ix+18,ii)*sc+aat(ix+18,mm)*cc return end subroutine ifc(wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) implicit none integer*4 ix,nn0,ii real*8 wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0) do 2 ii=1,12*nn0 2 wsu(ii)=0.0d0 do 1 ii=1,nn0 do 1 ix=1,9 wsat(ix,ii)=0.0d0 wcat(ix,ii)=0.0d0 wsgt(ix,ii)=0.0d0 wcgt(ix,ii)=0.0d0 wcaat(ix,ii)=0.0d0 wsaat(ix,ii)=0.0d0 wcaat(ix+9,ii)=0.0d0 wsaat(ix+9,ii)=0.0d0 wcaat(ix+18,ii)=0.0d0 1 wsaat(ix+18,ii)=0.0d0 write(6,*)'coefficients initiated' return end subroutine ifcr(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) implicit none integer*4 ix,nn0,ii,istep real*8 wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0) logical lex inquire(file='RESTART.SCR',exist=lex) if(.not.lex)then istep=0 call ifc(wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) write(6,*)'RESTART.SCR not found' else open(9,file='RESTART.SCR',form='unformatted') read(9)istep read(9)(wsu(ii),ii=1,12*nn0) read(9)((wsat(ix,ii),ix=1,9),ii=1,nn0) read(9)((wcat(ix,ii),ix=1,9),ii=1,nn0) read(9)((wsgt(ix,ii),ix=1,9),ii=1,nn0) read(9)((wcgt(ix,ii),ix=1,9),ii=1,nn0) read(9)((wsaat(ix,ii),ix=1,27),ii=1,nn0) read(9)((wcaat(ix,ii),ix=1,27),ii=1,nn0) close(9) write(6,*)'RESTART.SCR read' endif return end subroutine ifcw(istep,wsu,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nn0) implicit none integer*4 ix,nn0,ii,istep real*8 wsu(12*nn0),wsat(9,nn0),wcat(9,nn0),wsgt(9,nn0), 1wcgt(9,nn0),wsaat(27,nn0),wcaat(27,nn0) open(9,file='RESTART.SCR',form='unformatted') write(9)istep write(9)(wsu(ii),ii=1,12*nn0) write(9)((wsat(ix,ii),ix=1,9),ii=1,nn0) write(9)((wcat(ix,ii),ix=1,9),ii=1,nn0) write(9)((wsgt(ix,ii),ix=1,9),ii=1,nn0) write(9)((wcgt(ix,ii),ix=1,9),ii=1,nn0) write(9)((wsaat(ix,ii),ix=1,27),ii=1,nn0) write(9)((wcaat(ix,ii),ix=1,27),ii=1,nn0) close(9) write(6,*)'RESTART.SCR updated' return end SUBROUTINE READTTT(N,ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*N,3,3),G(3*N,3,3),A(3*N,3,3,3) READ(15,*) READ(15,*)NAT READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(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) RETURN END SUBROUTINE READTENI(N,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N,3,3),A(N,3,3) READ (15,*)NAT DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) 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,*) RETURN END subroutine readatms(io,NAT,iii,n) integer*4 n,NAT,iii(4,n),ia,io,ii do 1 ia=1,NAT 1 read(io,*)(iii(ii,ia),ii=1,1),(iii(ii,ia),ii=1,4) return end subroutine dtm(ii,NAT,u,n,x,y,z) implicit none integer*4 n,NAT,ii(4,n),ia,i1,i2,j1,j2,ix,jx,ixx real*8 u(3,3,n),x(*),v1(3),v2(3),v3(3),y(*),z(*),dxx c if(NAT.eq.1)then do 2 ix=1,3 do 2 jx=1,3 2 u(jx,ix,1)=0.0d0 do 3 jx=1,3 3 u(jx,jx,1)=1.0d0 else if(NAT.eq.2)then do 4 ia=1,NAT i1=ii(1,ia) i2=ii(2,ia) v1(1)=x(i2)-x(i1) v1(2)=y(i2)-y(i1) v1(3)=z(i2)-z(i1) v2(1)=0.0d0 v2(2)=0.0d0 v2(3)=0.0d0 dxx=dabs(v1(1)) ixx=1 if(dabs(v1(2)).lt.dxx)then ixx=2 dxx=dabs(v1(2)) endif if(dabs(v1(3)).lt.dxx)ixx=3 v2(ixx)=1.0d0 call make3v(v1,v2,v3) do 4 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 4 u(3,ix,ia)=v3(ix) else do 1 ia=1,NAT i1=ii(1,ia) i2=ii(2,ia) j1=ii(3,ia) j2=ii(4,ia) if(i1.gt.0.and.i2.gt.0.and.j1.gt.0.and.j2.gt.0)then v1(1)=x(i2)-x(i1) v2(1)=x(j2)-x(j1) v1(2)=y(i2)-y(i1) v2(2)=y(j2)-y(j1) v1(3)=z(i2)-z(i1) v2(3)=z(j2)-z(j1) call make3v(v1,v2,v3) do 13 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 13 u(3,ix,ia)=v3(ix) else do 14 ix=1,3 do 14 jx=1,3 14 u(jx,ix,ia)=0.0d0 u(1,1,ia)=1.0d0 u(2,2,ia)=1.0d0 u(3,3,ia)=1.0d0 endif 1 continue endif endif return end SUBROUTINE TTTT(ALPHA,A,G,NAT,ICO,xx,y,z,X) C Transforms A and G to local origins (ICO=1) c or common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3),xx(*), 1y(*),z(*),X(3,NAT) C C transfer into atomic coordinates: DO 5 J=1,NAT X(1,J)=xx(J)/0.529177D0 X(2,J)= y(J)/0.529177D0 5 X(3,J)= z(J)/0.529177D0 C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID) G(IIND,IA,IB)=G(IIND,IA,IB)+SUM IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0 2 CONTINUE c IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins' c IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' c IF(ICO.EQ.3)WRITE(6,*)' G set to zero' C DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM 1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0 3 CONTINUE c IF(ICO.EQ.1)WRITE(6,*)' A transformed into atomic origins' c IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' c IF(ICO.EQ.3)WRITE(6,*)' A set to zero' RETURN END subroutine report(s) character*(*) s write(6,*)s stop end FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF(I.EQ.1)then if(J.EQ.2.AND.K.EQ.3)then EPS= 1.0d0 else IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 endif else IF(I.EQ.2)then IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 else IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 endif endif RETURN END subroutine norm(v) implicit none real*8 v(*),n n=dsqrt(v(1)**2+v(2)**2+v(3)**2) if(n.gt.0.0d0)then v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n endif return end subroutine rg(n,x,y,z,zim,i12,ic,maxval) implicit none integer*4 n,ia,maxval,ic,ii real*8 x(*),y(*),z(*) integer*4 zim(*) integer*4 i12(maxval,n) if(ic.eq.1)then do 2 ia=1,n 2 read(10,*)zim(ia),x(ia),y(ia),z(ia),(i12(ii,ia),ii=1,7) else do 1 ia=1,n 1 read(10,*)zim(ia),x(ia),y(ia),z(ia) endif return end subroutine cctn_tinker2(NAT,atomic,x,y,z, 1R,ND,NDS,RS,INDBIG,IBONDSB,IBONDS,NBO,NBOB, 2PS,AS,ALPHAS,AAS,GS,AAT,GT, 3asname,psname,aaname,gname,alphaname,rname,NATNAME, 4NDNAME,IOX,luse,po,apo,gpo,lpol,lpolg,lpola, 1polname,gpolname,apolname,fr,natx) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind,atomic(*) dimension R(3,NAT), 1ND(NAT),NDS(NAT), 1INDBIG(IOX,NAT),IBONDSB(8,NAT),IBONDS(8,NAT),NBO(NAT), 1NBOB(NAT),PS(NAT,3,3),AS(NAT,3,3),ALPHAS(3*NAT,3,3), 3AAS(3*NAT,3,3,3),GS(3*NAT,3,3),AAT(3*NAT,3,3,3),GT(3*NAT,3,3), 1RS(3,NAT),cot(3), 1asname(IOX,natx,3,3),psname(IOX,natx,3,3), 1aaname(IOX,3*natx,3,3,3),gname(IOX,3*natx,3,3), 1alphaname(IOX,3*natx,3,3),rname(IOX,3,natx),NATNAME(IOX), 1polname(IOX,3,3),gpolname(IOX,3,3),apolname(IOX,3,3,3), 1NDNAME(IOX,natx) LOGICAL LWR,LABS,LVCD,LROA,LSTRICT,LINV, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM,lpol,lpolg,lpola logical lxfile CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname DIMENSION x(*),y(*),z(*),po(3,3),gpo(3,3),apo(3,3,3) common/cctnopts/NO,NATS, 1LWR,LABS,LVCD,IWG,LROA,LRAM,LSTRICT,LINV,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,n1,n2,okind logical luse(*) fr=0.0d0 apo=0.0d0 gpo=0.0d0 po=0.0d0 bohr=0.529177d0 c open(3,file='CCT.OUT') NATi=NAT c c NAT ... number of atoms in the tinekr molecule c NATi ... number of atoms in the tinker mol for transfer c NATS ... number of atoms in the small moelcule c Open file with option and overlaps and read it in: c For selected atoms only: CALL readopt(NAT,INDBIG,IOX,NO,NATi, 1LWR,LABS,LVCD,IWG,LROA,LRAM,LSTRICT,LINV,LAPT, 2LALF,LRDO,LNAMES,NAME,LPOLYMER,LHALTONERROR,n1,n2,okind, 3lpol,lpolg,lpola) c c if tensors constructed from named fragments, load them all: if(LNAMES)then DO 11 I=1,NO NE=80 do 10 NE=67,1,-1 10 if(NAME(I)(NE:NE).ne.' ')goto 12 12 basicname(1:NE)=NAME(I)(1:NE) xname=basicname(1:NE)//'.x' write(6,*)xname(1:NE+2) inquire(file=xname,exist=lxfile) if(.not.lxfile)then write(6,*)xname call report('Coordinates not found') endif open(2,file=xname) read(2,*) read(2,*)NATS if(NATS.gt.NAT)call report('NATS>NAT') rewind 2 CALL READRXS(NAT,NATS,RS,NDS,NBO,IBONDS) NATNAME(I)=NATS do 13 ia=1,NATS NDNAME(I,ia)=NDS(ia) do 13 ix=1,3 13 rname(I,ix,ia)=RS(ix,ia) 11 WRITE(6,*)NATS, ' atoms in the fragment' else OPEN(2,FILE='SMALL.X',STATUS='OLD') read(2,*) read(2,*)NATS if(NATS.gt.NAT)call report('NATS>NAT') rewind 2 CALL READRXS(NAT,NATS,RS,NDS,NBO,IBONDS) write(6,*)NATS,' atoms in SMALL.X' CLOSE(2) endif c names c i12 bond table for the small molecule,rewrite to big for all overlaps c BIG: i -----------ib NBOB c SMALL: INDBIG(1,i)-----------iss NBO do 2011 i=1,NAT 2011 NBOB(i)=0 do 201 i=1,NAT if(luse(i))then do 2012 io=1,NO if(NBOB(i).eq.0)then if(INDBIG(io,i).ne.0)then do 202 j=1,NBO(INDBIG(io,i)) iss=IBONDS(j,INDBIG(io,i)) ibb=0 do 203 ii=1,NAT if(INDBIG(io,ii).eq.iss)then ibb=ii goto 204 endif 203 continue 204 if(ibb.gt.0)then NBOB(i)=NBOB(i)+1 IBONDSB(NBOB(i),i)=ibb endif 202 continue endif c write(6,*)i,(IBONDSB(ii,i),ii=1,NBOB(i)) endif 2012 continue endif 201 continue do 101 i=1,NAT ND(i)=atomic(i) R(1,i)=x(i) R(2,i)=y(i) 101 R(3,i)=z(i) c c check consistency in atomic types: do 22 i=1,NO if(LNAMES)then NATS=NATNAME(i) do 23 ia=1,NATS 23 NDS(ia)=ndname(i,ia) endif DO 22 IAB=1,NATi IF(INDBIG(I,IAB).NE.0)THEN IF(ND(IAB).NE.NDS(ABS(INDBIG(I,IAB))))THEN WRITE(6,3300)I,IAB,ND(IAB),NDS(ABS(INDBIG(I,IAB))) 3300 FORMAT('Overlap ',I2,', atom ',I3,': ',2I4) CALL REPORT(' Different atomic types !') ENDIF IF(IABS(INDBIG(I,IAB)).gt.NATS)then write(3,*)'overlap ',i write(3,*)' NATS ',NATS write(3,*)'IND ',IAB write(3,*)'INDB ',INDBIG(I,IAB) write(6,*)'overlap ',i write(6,*)' NATS ',NATS write(6,*)'IND ',IAB write(6,*)'INDB ',INDBIG(I,IAB) c call report('Overflow NATS') ENDIF ENDIF 22 CONTINUE c c c Special section for named fragments: c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn if(lnames)then do 14 i=1,NO NATS=natname(i) NE=80 do 15 NE=80,1,-1 15 if(NAME(I)(NE:NE).ne.' ')goto 16 16 basicname=NAME(I)(1:NE) fname=basicname(1:NE)//'.fc' write(6,*)fname(1:NE+3) c do 19 ia=1,NATS do 19 ix=1,3 19 RS(ix,ia)=rname(i,ix,ia) c IF(LABS.OR.LVCD)THEN nname=basicname(1:NE)//'.ten' write(6,*)nname(1:NE+4) OPEN(15,FILE=nname,STATUS='OLD') CALL READTEN1(NAT,PS,AS,LAPT,RS,NATS) CLOSE(15) do 18 ia=1,NATS do 18 ix=1,3 do 18 jx=1,3 asname(i,ia,ix,jx)=AS(ia,ix,jx) 18 psname(i,ia,ix,jx)=PS(ia,ix,jx) ENDIF C IF(LROA.or.LRAM)THEN tname=basicname(1:NE)//'.ttt' write(6,*)tname(1:NE+4) OPEN(15,FILE=tname,STATUS='OLD') CALL READTTT(NAT,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(6,*)' A and G obtained from alpha as DO' CALL TTTT1(3*NAT,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*NAT,ALPHAS,AAS,GS,NATS,1,RS) ENDIF do 20 iax=1,3*NATS do 20 ix=1,3 do 20 jx=1,3 alphaname(i,iax,ix,jx)=ALPHAS(iax,ix,jx) gname(i,iax,ix,jx)=GS(iax,ix,jx) do 20 kx=1,3 20 aaname(i,iax,ix,jx,kx)=AAS(iax,ix,jx,kx) C Now small tensors are in atomic origins ENDIF c read alpha G A polarizability: if(lpol)call readpolal(po,basicname(1:NE)//'.pol',fr) if(lpolg)call readpolg(gpo,basicname(1:NE)//'.g.pol') if(lpola)call readpola(apo,basicname(1:NE)//'.a.pol') c transform G, A to small molecule center if(lpola.or.lpolg)then do 301 ix=1,3 cot(ix)=0.0d0 do 3011 ia=1,NATS 3011 cot(ix)=cot(ix)+rname(i,ix,ia) 301 cot(ix)=cot(ix)/bohr/dble(NATS) call TTT1(po,apo,gpo,1,cot,1) endif c c transcript alpha,G,A to named strings: do 302 ix=1,3 do 302 jx=1,3 polname(I,ix,jx)=po(ix,jx) gpolname(I,ix,jx)=gpo(ix,jx) do 302 kx=1,3 302 apolname(I,ix,jx,kx)=apo(ix,jx,kx) C 14 continue c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn else IF(LABS.OR.LVCD)THEN OPEN(15,FILE='SMALL.TEN',STATUS='OLD') CALL READTEN1(NAT,PS,AS,LAPT,RS,NATS) write(6,*)'SMALL.TEN read' CLOSE(15) ENDIF C IF(LROA.OR.LRAM)THEN OPEN(15,FILE='SMALL.TTT',STATUS='OLD') CALL READTTT(NAT,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO for small' CALL TTTT1(3*NAT,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT1(3*NAT,ALPHAS,AAS,GS,NATS,1,RS) ENDIF C Now small tensors are in atomic origins C IF(LRDO)THEN WRITE(3,*)' Alpha from SMALL_ALPHA.TTT for small' OPEN(15,FILE='SMALL_ALPHA.TTT',STATUS='OLD') CALL READTTT(NAT,ALPHAS,AAT,GT,NATS) CLOSE(15) ENDIF write(6,*)'Raman tensors read' endif if(lpol)call readpolal(po,'SMALL.POL',fr) if(lpolg)call readpolg(gpo,'SMALL.G.POL') if(lpola)call readpola(apo,'SMALL.A.POL') c transform G, A to small molecule center if(lpola.or.lpolg)then do 313 ix=1,3 cot(ix)=0.0d0 do 3131 ia=1,NATS 3131 cot(ix)=cot(ix)+RS(ix,ia) 313 cot(ix)=cot(ix)/bohr/dble(NATS) call TTT2(po,apo,gpo,1,cot,1) c call WRITEPOL(po,'POL.TTT',fr) c call WRITEPOL(gpo,'POL.G.TTT',fr) c call WRITEPOLA(apo,fr) endif ENDIF c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn c lnames c write(6,*)'CCTN options read in' close(3) return end subroutine trtena(x,y,z,ALPHA,G,AA,PB,AB,n,mxb,R,IND,RS,INDBIG, 1IBONDSB,NBOB,PS,AS,ALPHAS,AAS,GS,IOX,luse,asname,psname,aaname, 1gname,alphaname,rname,NATNAME,pol,gpo,apo,bigpol,biggpol,bigapol, 1lpol,lpolg,lpola,polname,gpolname,apolname,lcidr,lave, 1nm,lq,gq,aq,nmol,lqo,gqo,aqo,coo,opol,gopol,aopol,natx) IMPLICIT none integer*4 IOX,okind,nm,nmol,n,INDBIG(IOX,n),IBONDSB(8,n), 1NBOB(n),NATNAME(*),i,J,K,L,IX,IERR,IIND,IWG,MXB,N1,N2,NATS,NO, 1natx,IND(*) real*8 lq(*),gq(*),aq(*),R(3,n),RS(3,n),PS(n,3,3),AS(n,3,3), 1ALPHAS(3*n,3,3),AAS(3*n,3,3,3),GS(3*n,3,3),PB(n,3,3),AB(n,3,3), 1ALPHA(3*n,3,3),AA(3*n,3,3,3),G(3*n,3,3),asname(IOX,natx,3,3), 1psname(IOX,natx,3,3), 1aaname(IOX,3*natx,3,3,3),gname(IOX,3*natx,3,3), 1alphaname(IOX,3*natx,3,3),rname(IOX,3,natx),bigapol(3,3,3), 1pol(3,3),gpo(3,3),apo(3,3,3),bigpol(3,3),biggpol(3,3), 1polname(IOX,3,3),gpolname(IOX,3,3),apolname(IOX,3,3,3), 1opol(n,3,3),gopol(n,3,3),aopol(n,3,3,3),coo(n,3),x(*),y(*),z(*), 1lqo(*),gqo(*),aqo(*) LOGICAL LWR,LABS,LVCD,LROA,LSTRICT,LINV,LAPT,LALF,LRDO,LNAMES, 1LPOLYMER,LHALTONERROR,LRAM,luse(*),lpol,lpolg,lpola,lcidr,lave common/cctnopts/NO,NATS, 1LWR,LABS,LVCD,IWG,LROA,LRAM,LSTRICT,LINV,LAPT, 2LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,n1,n2,okind c c transcript current geometry: do 101 i=1,n R(1,i)=x(i) R(2,i)=y(i) 101 R(3,i)=z(i) c c transcript current geometry: IF(LABS.OR.LVCD)THEN DO 2 I=1,n DO 2 J=1,3 DO 2 K=1,3 PB(I,J,K)=0.0d0 2 AB(I,J,K)=0.0d0 ENDIF C IF(LROA.OR.LRAM)THEN DO 3 L=1,n DO 3 IX=1,3 IIND=3*(L-1)+IX DO 3 I=1,3 DO 3 J=1,3 ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 DO 3 K=1,3 3 AA(IIND,I,J,K)=0.0d0 ENDIF C CALL IMPROVE(R,RS,n,NATS,IND,INDBIG, 1NO,IBONDSB,NBOB,PB,AB,LABS,LVCD,AS,PS,IWG,LROA,LRAM, 1ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LNAMES, 3LHALTONERROR,n1,n2,okind,IERR,LWR,luse,IOX,mxb, 1asname,psname,aaname,gname,alphaname,rname,NATNAME, 1pol,gpo,apo,bigpol,biggpol,bigapol,lpol,lpolg,lpola, 1polname,gpolname,apolname,lcidr,lave, 1nm,lq,gq,aq,nmol,lqo,gqo,aqo,coo, 1opol,gopol,aopol,natx) c if(IERR.ne.0)then write(6,*)'MD point skipped' return endif c c c transform back: CALL TRTENVCD(n,PB,AB,R,n) IF(LROA.or.LRAM)THEN CALL TTTT1(3*n,ALPHA,AA,G,n,2,R) ENDIF return END SUBROUTINE IMPROVE(RB,RS,NAT,NATS,IND,INDBIG, 1NO,IBONDS ,NBOB,PB,AB,LABS,LVCD,AS,PS,IWG,LROA,LRAM, 2ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LNAMES, 3LHALTONERROR,n1,n2,okind,IERR,LWR,luse,IOX,mxb, 4aname,pname,aaname,gname,alphaname,rname,NATNAME, 1pol,gpol,apol,bigpol,biggpol,bigapol,lpol,lpolg,lpola, 1polname,gpolname,apolname,lcidr,lave, 1nm,lq,gq,aq,nmol,lqo,gqo,aqo,coo, 1opol,gopol,aopol,natx) implicit none INTEGER*4 mxb,IOX,okind,nmol,natx,IND(*), 1INDBIG(IOX,NAT),ITU,natname(IOX),iao,ic,id,ig,IOP,KX, 4IBONDS(8,NAT),NBOB(*),IOL(IOX),ISUM,I,IA,IAS,IEND,IERR,II, 1IIND,IIOII,INDS,IO,nu, 2IOIJ,IP,IPASS,IPO,ISTART,IWG,IX,IXX,J,JJ,IERC,noo, 3JX,K,KK,N1,N2,NAT,NATIO,NATS,NF,NO, 4ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,NAT),nats_o(IOX), 1nm real*8 RB(3,NAT),RS(3,NATS),TIJ(3),OTIJ(3), 1PB(NAT,3,3),AB(NAT,3,3),PS(NAT,3,3),AS(NAT,3,3),WF(IOX), 1ALPHA(3*NAT,3,3),ALPHAS(3*NAT,3,3),AA(3*NAT,3,3,3), 1AAS(3*NAT,3,3,3),G(3*NAT,3,3),GS(3*NAT,3,3),DISTM,SUM,PNEW,ANEW, 1ALPHANEW,GNEW,AANEW,DIST,POLD,AOLD,AD1,AD2,AD0,bohr,EPS, 6RM(IOX,3,3),aname(IOX,NAT,3,3),pname(IOX,NAT,3,3), 1aaname(IOX,3*natx,3,3,3),gname(IOX,3*natx,3,3), 2alphaname(IOX,3*natx,3,3),rname(IOX,3,natx), 1A(3,3),RMS,alst(3,3,3),glst(3,3,3),aast(3,3,3,3),pol(3,3), 1gpol(3,3),apol(3,3,3),bigpol(3,3),biggpol(3,3),bigapol(3,3,3), 1polname(IOX,3,3),gpolname(IOX,3,3),apolname(IOX,3,3,3), 1lq(*),gq(*),aq(*),opol(NAT,3,3), 1gopol(NAT,3,3),aopol(NAT,3,3,3),lqo(nmol*nm*9),gqo(nmol*nm*9), 1aqo(nmol*nm*27),coo(NAT,3) LOGICAL LABS,LVCD,LROA,LRAM,LSTRICT,LINV, 1LNAMES,LHALTONERROR,LWR,lpol,lpolg,lpola logical luse(*),lcidr,lave common/rmssom/RMS real*8,allocatable::XS(:,:),XB(:,:),zsum(:) allocate(XS(mxb,3),XB(mxb,3),zsum(IOX)) c IPASS=0 bohr=0.529177d0 c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp c group polarizabilities: opol=0.0d0 gopol=0.0d0 aopol=0.0d0 c polarizability derivatives in normal modes: lqo=0.0d0 gqo=0.0d0 aqo=0.0d0 if(nmol.ne.NO)call report('nmol <> NO') if(lpol)then if(IOX.gt.NAT)call report('IOX > n') c c total polarizability biggpol=0.0d0 bigpol=0.0d0 bigapol=0.0d0 c loop over overlaps: DO 551 IO=1,NO c make the transformation matrix for whole overlap: call CRMSo(NAT,INDBIG,IOX,NAT,IO,LNAMES,natname,RS,RB,IND, 2 rname,NATS,LINV,mxb,A,XS,XB,natx) c c calculate overlap center: do 53 ix=1,3 53 coo(IO,ix)=0.0d0 noo=0 nu=0 zsum(IO)=0.0d0 do 52 iao=1,NAT if(INDBIG(IO,iao).GT.0)then zsum(IO)=zsum(IO)+1.0d0 noo=noo+1 do 54 ix=1,3 54 coo(IO,ix)=coo(IO,ix)+RB(ix,iao) c c look at other overlaps, and find out if this atoms is part of them: ic=0 do 552 IOP=1,NO 552 if(IOP.ne.IO.and.INDBIG(IOP,iao).ne.0)ic=ic+1 if(ic.eq.0)nu=nu+1 endif 52 continue if(noo.gt.0)then do 55 ix=1,3 55 coo(IO,ix)=coo(IO,ix)/zsum(IO)/bohr endif if(lwr)write(6,6009)IO,noo,NATS,nu 6009 format('Overlap',i5,':',i5,' atoms transferred of',i5,',', 1 i5,' unique') c c if named fragments, fill alpha, G, A for each overlap: if(lnames)then call ctd(IO,polname,pol,IOX) if(lpolg)call ctd(IO,gpolname,gpol,IOX) if(lpola)call ctdd(IO,apolname,apol,IOX) endif c c rotate polarizabilites: c write(6,*)'IO',IO c call wrm(' A:',A) c call wrm(' pol:',pol) if(lpol)call rpp(NAT,IO,opol,pol,A) if(lpolg)call rpp(NAT,IO,gopol,gpol,A) if(lpola)call rppp(NAT,IO,aopol,apol,A) c write(6,*)'opol:' c write(6,444)((opol(IO,IX,JX),IX=1,3),JX=1,3) c44 format(3g12.4) c rotate polarizability derivatibes (normal mode) lq,gq,aq, c put them to lqo,gqo,aqo if(lave)call rpd(IO,nm,A,lq,gq,aq,lqo,gqo,aqo) c c construct alpha,G,A of all system: do 5519 IX=1,3 do 5519 JX=1,3 c c plain sum: bigpol(IX,JX) =bigpol(IX,JX) +opol(IO,IX,JX) biggpol(IX,JX)=biggpol(IX,JX)+gopol(IO,IX,JX) do 55191 KX=1,3 55191 bigapol(IX,JX,KX)=bigapol(IX,JX,KX)+aopol(IO,IX,JX,KX) c c polarizability part for origin dependence of G: DO 618 IG=1,3 DO 618 ID=1,3 618 biggpol(IX,JX)=biggpol(IX,JX) 1 -0.5d0*EPS(JX,IG,ID)*coo(IO,IG)*opol(IO,IX,ID) c c polarizability part for origin dependence of A: do 5519 KX=1,3 if(JX.eq.IX)bigapol(IX,JX,KX)=bigapol(IX,JX,KX) 1 -coo(IO,1)*opol(IO,KX,1) 2 -coo(IO,2)*opol(IO,KX,2) 3 -coo(IO,3)*opol(IO,KX,3) 5519 bigapol(IX,JX,KX)=bigapol(IX,JX,KX) 1 +1.5d0*(coo(IO,JX)*opol(IO,KX,IX)+coo(IO,IX)*opol(IO,KX,JX)) 551 continue c correction of polarizability for mutual interaction of groups: call mkcorrect(NAT,NO,coo,opol,gopol,aopol,lpola,lpolg, 1 bigpol,biggpol,bigapol) if(lcidr)then c call WRITEPOL(bigpol,'POL.TTT',fr) c if(lpolg)call WRITEPOL(biggpol,'POL.G.TTT',fr) c if(lpola)call WRITEPOLA(bigapol,fr) endif if(LWR)then write(6,6301) 6301 format('Isotropic Polarizabilities:') do 555 IO=1,NO 555 write(6,6302)IO,(opol(IO,1,1)+opol(IO,2,2)+opol(IO,3,3))/3.0d0 6302 format(i3,f10.3) write(6,6303)NO 6303 format('Overlap centers:',/,I4) do 556 IO=1,NO 556 write(6,6304)(coo(IO,ix)*bohr,ix=1,3) 6304 format('6 ',3f12.4) endif endif c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp c if(lpol)then C C Atom IA istart=1 iend=NAT if(n1.gt.0)then istart=n1 iend=n2 endif DO 20 IA= istart,iend if(luse(IA))then c IF(okind.ne.2.or.(IA.eq.istart))THEN C DO 22 IX=1,3 22 TIJ(IX)=RB(IX,IA) C C Look at all overlaps io, find that which fits best the i-j pair: DISTM=1.0d5 NF=0 DO 21 IO=1,NO C C Check, if io contains i: IF(INDBIG(IO,IA).GT.0)THEN NF=NF+1 IOL(NF)=IO C C Calculate center of the overlap and distance from i-j center c i.........TIJ.............j c ----------OTIJ--------------------- DO 23 IX=1,3 23 OTIJ(IX)=0.0D0 NATIO=0 DO 24 II=1,NAT IIOII=INDBIG(IO,II) C If Lstrict, take only transformed atoms IF(LSTRICT.AND.IIOII.LT.0)IIOII=0 IIOII=ABS(IIOII) IF(IIOII.NE.0)THEN NATIO=NATIO+1 DO 25 IX=1,3 25 OTIJ(IX)=OTIJ(IX)+RB(IX,II) ENDIF 24 CONTINUE c c If required, calculate RMS overlap and use this criterium instead c of the distances,remember rotation matrix (RM) if(okind.gt.0)then call CRMS(IPASS,NAT,INDBIG,IOX,IO, 1 IA,IA,NBOB,okind,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS,IERR, 1 A,ITU,mxb,XS,XB,natx) if(IERR.ne.0)return do 62 IX=1,3 do 62 JX=1,3 62 RM(NF,IX,JX)=A(IX,JX) ipass_o(NF)=IPASS WF(NF)=DSQRT(RMS) ierc_o(NF)=IERC itu_o(NF)=ITU do 58 I=1,ITU 58 ind_o(NF,I)=IND(I) do 59 I=1,3 do 59 J=1,3 59 RM(NF,I,J)=A(I,J) nats_o(NF)=NATS IF(RMS.LT.DISTM)THEN DISTM=RMS IOIJ=IO ENDIF else c Use distances: DO 26 IX=1,3 26 OTIJ(IX)=OTIJ(IX)/DBLE(NATIO) DIST=(OTIJ(1)-TIJ(1))**2+(OTIJ(2)-TIJ(2))**2 1 +(OTIJ(3)-TIJ(3))**2 IF(DIST.LT.DISTM)THEN DISTM=DIST IOIJ=IO ENDIF WF(NF)=DSQRT(DIST) endif ENDIF 21 CONTINUE if(LWR)then write(6,3009)IA 3009 format(' Atom',I6,$) if(NF.eq.0)then write(6,3008) 3008 format(' No overlap') else write(6,3003)NF,IOIJ 3003 format(i5,' overlaps found; best is # ',i4) endif endif c C Calculate the weighting factors: SUM=0.0d0 DO 43 I=1,NF 43 IF(WF(I).GT.1.0d-9)SUM=SUM+1.0d0/WF(I) DO 42 I=1,NF IF(WF(I).LE.1.0d-9)THEN WF(1)=1.0d0 NF=1 GOTO 44 ENDIF WF(I)=1.0d0/WF(I)/SUM 42 CONTINUE 44 CONTINUE c c cccccccccccccccccccccccccccccccccc ENDIF C C Loop over possible overlaps: ISUM=0 DO 202 IP=1,NF IPO=IOL(IP) IF(IWG.EQ.0)THEN IF(IPO.NE.IOIJ)GOTO 202 WF(IP)=1.0d0 ENDIF ISUM=ISUM+1 C C Find the atomic groups to be placed on together for IA: c and the rotation matrix if(okind.eq.0)then call CRMS(IPASS,NAT,INDBIG,IOX,IPO, 1 IA,IA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS,IERR, 1 A,ITU,mxb,XS,XB,natx) if(IERR.ne.0)return else IPASS=ipass_o(IP) IERC=ierc_o(IP) ITU=itu_o(IP) do 56 I=1,ITU 56 IND(I)=ind_o(IP,I) do 57 I=1,3 do 57 J=1,3 57 A(I,J)=RM(IP,I,J) NATS=nats_o(IP) endif if(IERC.ne.0)goto 202 if(LWR)then write(6,3004)IPO,WF(IP),ITU 3004 format('Overlap ',I2,': weighting factor: ',F7.4,';', 1 I4,' atom bed:',/, 2 ' Big :',$) do 60 I=1,ITU 60 write(6,3007)IND(I) 3007 format(I3,$) write(6,*) write(6,3005) 3005 format( ' Small:',$) do 61 I=1,ITU 61 write(6,3007)INDBIG(IPO,IND(I)) write(6,*) endif C IF(LINV)THEN DO 50 I=1,3 DO 50 J=1,3 50 A(I,J)=-A(I,J) ENDIF C C C Transform the dipole derivatives IF(LABS.OR.LVCD)THEN IAS=INDBIG(IPO,IA) DO 40 IX=1,3 DO 40 JX=1,3 POLD=PB(IA,IX,JX) AOLD=AB(IA,IX,JX) C Zero out at the start of the summation IF(ISUM.EQ.1)THEN PB(IA,IX,JX)=0.0d0 AB(IA,IX,JX)=0.0d0 ENDIF PNEW=0.0d0 ANEW=0.0d0 DO 41 II=1,3 AD1=A(IX,II) DO 41 JJ=1,3 AD2=AD1*A(JX,JJ) if(LNAMES)then PNEW=PNEW+pname(IPO,IAS,II,JJ)*AD2 ANEW=ANEW+aname(IPO,IAS,II,JJ)*AD2 else PNEW=PNEW+PS(IAS,II,JJ)*AD2 ANEW=ANEW+AS(IAS,II,JJ)*AD2 endif 41 continue PB(IA,IX,JX)=PNEW*WF(IP)+PB(IA,IX,JX) AB(IA,IX,JX)=ANEW*WF(IP)+AB(IA,IX,JX) IF (.NOT.LVCD)AB(IA,IX,JX)=AOLD IF (.NOT.LVCD)ANEW=AOLD IF (.NOT.LABS)PB(IA,IX,JX)=POLD IF (.NOT.LABS)PNEW=POLD 40 continue ENDIF C C Transform the ROA and Raman polarizability derivatives IF(LROA.or.LRAM)THEN IAS=INDBIG(IPO,IA) if(lwr)then call wrt2(IA,NAT,ALPHA,'Old Alpha:') call wrt2(IA,NAT,G,'Old G:') call wrt3(IA,NAT,AA,'Old A:') endif if(LNAMES)then do 599 IXX=1,3 INDS=3*(IAS-1)+IXX do 599 II=1,3 do 599 JJ=1,3 do 5993 KK=1,3 5993 aast(IXX,II,JJ,KK)=aaname(IPO,INDS,II,JJ,KK) alst(IXX,II,JJ)=alphaname(IPO,INDS,II,JJ) 599 glst(IXX,II,JJ)=gname(IPO,INDS,II,JJ) else do 5991 IXX=1,3 INDS=3*(IAS-1)+IXX do 5991 II=1,3 do 5991 JJ=1,3 do 5994 KK=1,3 5994 aast(IXX,II,JJ,KK)=AAS(INDS,II,JJ,KK) alst(IXX,II,JJ)=ALPHAS(INDS,II,JJ) 5991 glst(IXX,II,JJ)=GS(INDS,II,JJ) if(lwr)then call wrt2(IAS,NAT,ALPHAS,'Small Alpha:') call wrt2(IAS,NAT,GS,'Small G:') call wrt3(IAS,NAT,AAS,'Small A:') endif endif IF(LRAM)THEN DO 46 IX=1,3 IIND=3*(IA-1)+IX DO 46 I=1,3 DO 46 J=1,3 C Zero out at the start of the summation IF(ISUM.EQ.1)ALPHA(IIND,I,J)=0.0d0 ALPHANEW=0.0d0 DO 47 IXX=1,3 AD0=A(IX,IXX) DO 47 II=1,3 AD1=AD0*A(I,II) DO 47 JJ=1,3 47 ALPHANEW=ALPHANEW+alst(IXX,II,JJ)*AD1*A(J,JJ) 46 ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+ALPHANEW*WF(IP) ENDIF IF(LROA)THEN c Rotation of G: DO 461 IX=1,3 IIND=3*(IA-1)+IX DO 461 I=1,3 DO 461 J=1,3 C Zero out at the start of the summation IF(ISUM.EQ.1)G(IIND,I,J)=0.0d0 GNEW=0.0d0 DO 471 IXX=1,3 AD0=A(IX,IXX) DO 471 II=1,3 AD1=AD0*A(I,II) DO 471 JJ=1,3 471 GNEW=GNEW+glst(IXX,II,JJ)*AD1*A(J,JJ) 461 G(IIND,I,J)=G(IIND,I,J)+GNEW*WF(IP) c c Rotation of A: DO 48 IX=1,3 IIND=3*(IA-1)+IX DO 48 I=1,3 DO 48 J=1,3 DO 48 K=1,3 C Zero out at the start of the summation IF(ISUM.EQ.1)AA(IIND,I,J,K)=0.0d0 AANEW=0.0d0 DO 49 IXX=1,3 AD0=A(IX,IXX) INDS=3*(IAS-1)+IXX DO 49 II=1,3 AD1=AD0*A(I,II) DO 49 JJ=1,3 AD2=AD1*A(J,JJ) DO 49 KK=1,3 49 AANEW=AANEW+aast(IXX,II,JJ,KK)*AD2*A(K,KK) 48 AA(IIND,I,J,K)=AA(IIND,I,J,K)+AANEW*WF(IP) ENDIF if(lwr)then call wrt2(IA,NAT,ALPHA,'New Alpha:') call wrt2(IA,NAT,G,'New G:') call wrt3(IA,NAT,AA,'New A:') endif ENDIF 202 CONTINUE c loop over overlaps, IP endif c use IA 20 CONTINUE c IA NF C RETURN C END c SUBROUTINE READRXS(N0,N,R,ND,NB,IBONDS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),NB(N0) integer*4 BT(7),IBONDS(8,N0) READ(2,*) read(2,*)N do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(BT(k),k=1,7) NB(i)=0 do 1 k=1,7 if(BT(k).ne.0)then NB(i)=NB(i)+1 IBONDS(NB(i),i)=BT(k) endif 1 continue RETURN END SUBROUTINE DOMATRIX(IERR,A,XB,XS,mxb,ITU) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3),A(3,3), 1XS(mxb,3),XB(mxb,3) A=0.0d0 IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.1d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(mxb,ANG,ITU,A,XS,XB) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(mxb,ANG,ITU,A,XS,XB) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(mxb,ANG,ITU,A,XS,XB) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(mxb,ANG,ITU,A,XS,XB) FTOL=1.0d-7 ITMAX=5000 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)then RETURN endif IF(ITER.EQ.ITMAX)THEN write(6,*)'mxb',mxb write(6,*)'ITU',ITU write(6,*)1,XS(1,1),XS(1,2),XS(1,3),XB(1,1),XB(1,2),XB(1,3) write(6,*)mxb,XS(mxb,1),XS(mxb,2),XS(mxb,3), 1 XB(mxb,1),XB(mxb,2),XB(mxb,3) write(6,*)'A',A WRITE(6,*)' Rotation has not converged !',ITER IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(mxb,PR,ITU,A,XS,XB) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(mxb,PRR,ITU,A,XS,XB) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(mxb,PRR,ITU,A,XS,XB) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(mxb,PR,ITU,A,XS,XB) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(mxb,U,n,A,XS,XB) implicit none INTEGER*4 mxb,n,I REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,U(3),FU, 1A(3,3),XS(mxb,3),XB(mxb,3),RMS common/rmssom/RMS S1=SIN(U(1)) S2=SIN(U(2)) S3=SIN(U(3)) C1=COS(U(1)) C2=COS(U(2)) C3=COS(U(3)) C u1 = theta C u2 = phi C u3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 RMS=0.0d0 DO 1 I=1,n TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 RMS=RMS+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=RMS RETURN END C c SUBROUTINE SREPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(3,3000) WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(3,*)RS WRITE(6,*)RS WRITE(3,3001) WRITE(6,3001) 3001 FORMAT(80(1H*),/,/, 1'PROGRAM CONTINUES BECAUSE LHALTONERROR DISABLED') END C SUBROUTINE READTEN1(N0,P,A,LAPT,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3),X(3,NAT) real*8,allocatable::PV(:,:,:),R(:,:),AI(:,:,:),AJ(:,:,:) LOGICAL LAPT allocate(PV(NAT,3,3),R(3,NAT),AI(NAT,3,3),AJ(NAT,3,3)) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR READ (15,*) NOAT NAT=NOAT WRITE(3,*)NAT,' atoms' WRITE(*,*)NAT,' atoms' C C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) WRITE(6,*)' Dipole derivatives read - in ' C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) AJ(L,J,I)=0.0d0 PV(L,J,I)=P(L,J,I) 2203 AI(L,J,I)=SUM WRITE(6,*)' APT part of AAT constructed' WRITE(3,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (AI(L,J,I),I=1,3) DO 230 L=1,NOAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NOAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) WRITE(6,*)' VCD parameters read in ...' WRITE(3,*)' VCD parameters read in ...' C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) WRITE(6,*)'DOG used' WRITE(3,*)'DOG used' C C Total axial tensor = electronic + nuclear: DO 4 L=1,NAT DO 4 I=1,3 DO 4 J=1,3 4 A(L,I,J)=AI(L,I,J)+AJ(L,I,J) C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(3,*)' AAT shifted to atomic origin' WRITE(6,*)' AAT shifted to atomic origin' RETURN END C SUBROUTINE TRTENVCD(N0,P,A,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3),X(3,NAT) real*8,allocatable::R(:,:) allocate(R(3,NAT)) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM RETURN END C SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA DIMENSION VCD(NAT,3,3),VEL(NAT,3,3),DIP(NAT,3,3), 1C(3,NAT) DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END C SUBROUTINE TTTT1(N,ALPHA,A,G,NAT,ICO,R) C Transforms A and G to local origins (ICO=1) or c common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT) real*8,allocatable::X(:,:) allocate(X(3,NAT)) C C transfer into atomic coordinates: DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I,J)/0.529177D0 C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID) G(IIND,IA,IB)=G(IIND,IA,IB)+SUM IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0 2 CONTINUE c IF(ICO.EQ.1)WRITE(3,*)' G transformed into atomic origins' c IF(ICO.EQ.2)WRITE(3,*)' G transformed into common origin' c IF(ICO.EQ.3)WRITE(3,*)' G set to zero' C DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM 1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0 3 CONTINUE c IF(ICO.EQ.1)WRITE(3,*)' A transformed into atomic origins' c IF(ICO.EQ.2)WRITE(3,*)' A transformed into common origin' c IF(ICO.EQ.3)WRITE(3,*)' A set to zero' RETURN END c SUBROUTINE readopt(n,INDBIG,IOX,NO, 1NAT,LWR,LABS,LVCD,IWG,LROA,LRAM,LSTRICT,LINV, 2LAPT,LALF,LRDO,LNAMES,NAME,LPOLYMER,LHALTONERROR, 2n1,n2,okind,lpol,lpolg,lpola) implicit none INTEGER*4 n,INDBIG,IOX,NO,nt, 1NAT,n1,n2,okind,I,IA,IAB,IC,IEND,II,IO,ISMALL,IWG,JA,NDD DIMENSION INDBIG(IOX,n) 1,NDD(20) CHARACTER*7 STR,STRC CHARACTER*80 NAME(IOX) LOGICAL LWR,LABS,LVCD,LROA,LSTRICT,LINV,LAPT, 1LALF,LRDO,LPOLYMER,LNAMES,LHALTONERROR,LRAM,lpol,loldformat, 1lpolg,lpola,lnew integer*4,allocatable::bb(:,:) c n1=0 n2=0 okind=0 LRDO=.FALSE. lpol=.FALSE. lpola=.FALSE. lpolg=.FALSE. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LHALTONERROR=.TRUE. lnew=.false. loldformat=.true. LROA=.FALSE. LRAM=.FALSE. LSTRICT=.TRUE. LINV=.FALSE. LNAMES=.FALSE. IWG=1 OPEN(2,FILE='CCT.INP',STATUS='OLD') 9 READ(2,2000)STR 2000 FORMAT(A7) WRITE(6,2000)STR IF(STR(1:4).EQ.'N1N2') READ(2,*)n1,n2 IF(STR(1:3).EQ.'OKI') READ(2,*)okind IF(STR(1:3).EQ.'IWG') READ(2,*)IWG IF(STR(1:6).EQ.'LNAMES') READ(2,*)LNAMES IF(STR(1:4).EQ.'LINV') READ(2,*)LINV IF(STR(1:4).EQ.'GPOL') READ(2,*)lpolg IF(STR(1:4).EQ.'APOL') READ(2,*)lpola IF(STR(1:5).EQ.'ALPOL') READ(2,*)lpol IF(STR(1:4).EQ.'LAPT') READ(2,*)LAPT IF(STR.EQ.'LSTRICT') READ(2,*)LSTRICT IF(STR(1:3).EQ.'LWR') READ(2,*)LWR IF(STR(1:4).EQ.'LALF') READ(2,*)LALF IF(STR(1:4).EQ.'LRDO') READ(2,*)LRDO IF(STR(1:4).EQ.'LRAM') READ(2,*)LRAM IF(STR(1:4).EQ.'LROA') READ(2,*)LROA IF(STR(1:4).EQ.'LABS') READ(2,*)LABS IF(STR(1:4).EQ.'LVCD') READ(2,*)LVCD IF(STR(1:4).EQ.'LNEW') READ(2,*)lnew IF(STR(1:7).EQ.'LHALTON')READ(2,*)LHALTONERROR IF(STR(1:7).EQ.'OLDFORM')READ(2,*)loldformat BACKSPACE 2 READ(2,2000)STRC if(STR.ne.STRC)goto 9 NO=0 c c older laborius definition of overlap: IF(STR.NE.'POLYMER')then call report('Only POLYMER allowed') else LPOLYMER=.true. READ(2,*)NO WRITE(3,*)NO,' overlaps found, atoms assigned automatically ...' IF(NO.GT.IOX)CALL REPORT(' NO > IOX ! Change the dimensions.') DO 2 I=1,NO READ(2,*)IC IF(IC.NE.I)THEN WRITE(6,*)IC,I,' nat = ',NAT CALL REPORT('Error in reading ') ENDIF IF(LNAMES)READ(2,2900)NAME(I) IF(LNAMES)write(6,2900)NAME(I) 2900 FORMAT(a80) c write(6,*)NAT c newest format if(lnew)then DO 52 IAB=1,NAT 52 INDBIG(I,IAB)=0 read(2,*)nt allocate(bb(2,nt)) backspace 2 read(2,*)nt,((bb(IO,IAB),IO=1,2),IAB=1,nt) do 53 IAB=1,nt 53 INDBIG(I,bb(1,IAB))=bb(2,IAB) deallocate(bb) else if(loldformat)then READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT) 2001 FORMAT(20I3) DO 5 IAB=1,NAT 5 INDBIG(I,IAB)=0 READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT) else READ(2,*)(INDBIG(I,IAB),IAB=1,NAT) DO 51 IAB=1,NAT 51 INDBIG(I,IAB)=0 READ(2,*)(INDBIG(I,IAB),IAB=1,NAT) endif endif 2 CONTINUE c c Control Output: DO 4 IO=1,NO WRITE(3,2002)IO 2002 FORMAT(' Overlap number ',I5) IA=0 7 IA=IA+1 IEND=IA+19 IF(IEND.GT.NAT)IEND=NAT IF(IA.GT.NAT)GOTO 6 DO 8 II=1,(IEND-IA+1) 8 NDD(II)=II+IA-1 WRITE(3,2004)(NDD(II),II=1,IEND-IA+1) 2004 FORMAT(20I3) IA=IEND GOTO 7 6 WRITE(3,2001)(INDBIG(IO,IAB),IAB=1,NAT) 4 CONTINUE c c check, that atoms are not defined twice in the overlaps: do io=1,no do ia=1,nat ismall=indbig(io,ia) if(ismall.ne.0)then do ja=ia+1,nat if(ismall.eq.indbig(io,ja))then write(3,*)'overlap ',io write(3,*)'big ',ia,ja write(3,*)'small ',ismall call report('umbiqous overlap') endif enddo endif enddo enddo ENDIF CLOSE(2) write(6,6011)n1,n2,okind,IWG,LNAMES,LINV,LPOL,LAPT, 1LSTRICT,LWR,LALF,LRDO,LRAM,LROA,LABS,LVCD, 1LHALTONERROR,loldformat 6011 format(/,' CTT Options',/,' ***********',/, 1' N1 N2 ',2i6,/, 1' OKIND ',i6,/, I' IWG ',i6,/, I' LNAMES ',l6,/, I' LINV ',l6,/, I' LPOL ',l6,/, I' LAPT ',l6,/, I' LSTRICT ',l6,/, I' LWR ',l6,/, I' LALFA ',l6,/, I' LRDO ',l6,/, I' RAM ROA ABS VCD ',4l6,/, I' HALT ON ERROR ',l6,/, I' OLD FORMAT ',l6,/,/) RETURN END C subroutine CRMS(IPASS,NAT,INDBIG,IOX,IPO, 1IA,JA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS,IERR, 1A,ITU,mxb,XS,XB,natx) c find covalently bond atoms to IA..JA, oevrlap and make the c rotation matrix c IC=0 ... take atoms two bonds away only if directly c bonded do not suffice c IC=1 ... include two bonds away by default c IC=2 ... include all transfered atoms implicit none integer*4 IA,IPASS,II,NAT,INDBIG(IOX,NAT),IBA,IBONDS(8,NAT), 1NBOB(*),JB,IIA,NATS,natname(*),IX,IND(*),iu,IERC,I,IC,IERR, 2IPO,J,JA,ITU,IOX,mxb,natx real*8 RS(3,NAT),RB(3,NAT),TB(3),TS(3),rname(IOX,3,natx), 1A(3,3),XS(mxb,3),XB(mxb,3) logical LNAMES,LHALTONERROR,LSTRICT,LINV IERC=1 IPASS=0 666 ITU=0 IPASS=IPASS+1 DO 27 II=1,NAT if(INDBIG(IPO,II).NE.0)then C Take IA and JA right away in the list: IF(II.EQ.IA.OR.II.EQ.JA.or.IC.eq.2)GOTO 277 C Take also atoms connected to IA or JA: DO 29 IBA=1,NBOB(II) 29 IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA)GOTO 277 C For second pass, take also atoms two bonds away: IF(IPASS.EQ.2.or.IC.eq.1)THEN DO 31 IBA=1,NBOB(II) IIA=IBONDS(IBA,II) DO 31 JB=1,NBOB(IIA) 31 IF(IBONDS(JB,IIA).EQ.IA.OR.IBONDS(JB,IIA).EQ.JA)GOTO 277 ENDIF c c This atom is of no use, try another one: GOTO 27 277 IF(INDBIG(IPO,II).LT.0.AND.LSTRICT)GOTO 27 C Include atom II into the transformation set: ITU=ITU+1 IF(ITU.GT.mxb)CALL REPORT('ITU > mxb !') if(LNAMES)NATS=natname(IPO) DO 28 IX=1,3 if(LNAMES)then XS(ITU,IX)=rname(IPO,IX,ABS(INDBIG(IPO,II))) else XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II))) endif 28 XB(ITU,IX)=RB(IX,II) IND(ITU)=II endif 27 CONTINUE c c IF(ITU.LT.3)WRITE(3,*)' ITU < 3 !' IF(ITU.LT.3.AND.IPASS.EQ.1.and.IC.eq.0)GOTO 666 if(ITU.LT.3)then write(6,*)'Only following atoms found in the big molecule:' write(6,*)(IND(iu) ,iu=1,ITU) write(6,*)'Pair',IA,JA write(6,*)'pass ic',IPASS,IC if(LHALTONERROR)then CALL REPORT(' ITU < 3 !') else CALL SREPORT(' ITU < 3 !') return endif endif c c C Calculate center of the coordinates that will be rotated: DO 33 I=1,3 TB(I) =0.0d0 33 TS(I)=0.0d0 DO 34 II=1,ITU DO 34 J=1,3 TB(J)=TB(J)+XB(II,J) 34 TS(J)=TS(J)+XS(II,J) DO 38 J=1,3 TB(J)=TB(J)/DBLE(ITU) 38 TS(J)=TS(J)/DBLE(ITU) C Subtract mass-center: DO 35 II=1,ITU DO 35 J=1,3 XS(II,J)=XS(II,J)-TS(J) 35 XB(II,J)=XB(II,J)-TB(J) C IF(LINV)THEN DO 39 II=1,ITU DO 39 J=1,3 39 XS(II,J)=-XS(II,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR,A,XB,XS,mxb,ITU) IF(IERR.EQ.1.AND.IPASS.EQ.1)THEN write(6,*)'domat error,try other pass' GOTO 666 ENDIF IF(IERR.EQ.1)then if(LHALTONERROR)then CALL REPORT('Program cannot continue') else CALL SREPORT('Program cannot continue') return endif endif IERC=0 return end c ==================================================================== subroutine protr(x,y,z,ALPHA,PB,N,ic,luse) c project translations/rotations from the electric tensors IMPLICIT none integer*4 I,J,K,IA,IX,IIA,ic,IIX,N,II,nc,jx,kx real*8 x(*),y(*),z(*),PB(N,3,3),ALPHA(3*N,3,3),r(3,N), 1t(3*N),PT(N,3,3),AT(3*N,3,3) real*8,allocatable::A(:,:) logical luse(*) allocate(A(3*N,3*N)) c c N is number of atoms here nc=0 do 1 i=1,N if(luse(i))then nc=nc+1 r(1,nc)=x(i) r(2,nc)=y(i) r(3,nc)=z(i) do 2 ix=1,3 do 2 jx=1,3 PT(nc,ix,jx)=PB(i,ix,jx) do 2 kx=1,3 2 AT(ix+(nc-1)*3,jx,kx)=ALPHA(ix+(i-1)*3,jx,kx) endif 1 continue CALL DOMA(r,A,nc) if(ic.eq.1)call rules(AT,nc,N) c DO 6 J=1,3 DO 11 I=1,3*nc t(I)=0.0d0 DO 11 IIA=1,nc DO 11 IIX=1,3 11 t(I)=t(I)+PT(IIA,IIX,J)*A(I,IIX+3*(IIA-1)) do 6 IA=1,nc do 6 IX=1,3 6 PT(IA,IX,J)=t(IX+3*(IA-1)) c DO 3 J=1,3 DO 3 K=1,3 DO 33 I=1,3*nc t(I)=0.0d0 DO 33 II=1,3*nc 33 t(I)=t(I)+AT(II,J,K)*A(I,II) do 3 I=1,3*nc 3 AT(I,J,K)=t(I) if(ic.eq.1)write(6,*)'projection ok' if(ic.eq.1)call rules(AT,nc,N) nc=0 do 4 i=1,N if(luse(i))then nc=nc+1 do 5 ix=1,3 do 5 jx=1,3 PB(i,ix,jx)=PB(nc,ix,jx) do 5 kx=1,3 5 ALPHA(ix+(i-1)*3,jx,kx)=AT(ix+(nc-1)*3,jx,kx) endif 4 continue RETURN END C ============================================================ SUBROUTINE DOMA(r,A,nat) implicit none integer*4 n,nat,N6,ICOL,I,J,II real*8 r(3,nat),A(3*nat,3*nat),S real*8,allocatable::TEM(:,:) allocate(TEM(3*nat,3*nat)) n=3*nat DO 12 I=1,n DO 12 J=1,n 12 A(I,J)=0.0d0 DO 1 I=1,nat A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-r(2,I) A(3*I-1,4)= r(1,I) A(3*I-1,5)=-r(3,I) A(3*I ,5)= r(2,I) A(3*I-2,6)= r(3,I) 1 A(3*I ,6)=-r(1,I) N6=6 CALL ORT(A,n,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1.0d-10)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END C ============================================================ SUBROUTINE ORT(A,NDIM,N6,NAT) implicit none integer*4 NDIM,N6,NAT,NAT3,I,J,IC,ICP REAL*8 A(NDIM,NDIM),C NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.0.0d0)C=1.0d0/DSQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 C=0.0d0 DO 5 J=1,NAT3 5 C=C+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-C*A(J,ICP) C Normalize : C=0.0d0 DO 7 J=1,NAT3 7 C=C+A(J,IC)**2 IF(C.GT.0.0d0)C=1.0d0/DSQRT(C) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*C RETURN END SUBROUTINE rules(ALPHA,NAT,N0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*N0,3,3),ap(3,3) character*2 xyz(3) data xyz/' X',' Y',' Z'/ c write(6,6077)' Sum rules-alpha:' 6077 format(a17) do 11 ii=1,3 do 21 ix=1,3 do 21 jx=1,3 ap(ix,jx)=0.0d0 do 21 ia=2,NAT 21 ap(ix,jx)=ap(ix,jx)+ALPHA((ia-1)*3+ii,ix,jx) write(6,6017)xyz(ii) 6017 format(a2) write(6,600)(( -ap(ix,jx),jx=1,3),ix=1,3) 600 format(9F8.3) 11 write(6,600)((ALPHA(ii,ix,jx),jx=1,3),ix=1,3) return end SUBROUTINE readpolal(a,fn,w) character*(*)fn real*8 a(3,3),w integer*4 i,j character*80 s80 open(90,file=fn) w=0.085645d0 read(90,7021)s80 7021 format(a80) read(s80(22:len(s80)),*,end=99,err=99)w 99 read(90,*) read(90,*)((a(i,j),i=1,j),j=1,3) close(90) do 1 j=1,3 do 1 i=1,j 1 a(j,i)=a(i,j) write(6,*)fn//' read' return end SUBROUTINE readpolg(a,fn) character*(*)fn real*8 a(3,3) integer*4 i,j open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(j,i),i=1,3),j=1,3) close(90) write(6,*)fn//' read' return end c SUBROUTINE readpola(a,fn) character*(*)fn integer*4 i,j,k real*8 a(3,3,3) open(90,file=fn) read(90,*) read(90,*) do 1 k=1,3 read(90,*)((a(i,j,k),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(j,i,k)=a(i,j,k) close(90) write(6,*)fn//' read' return end c subroutine dc(d,lg,dw,wmin,nw) c predefine convolution function implicit none real*8 dd,ff,w,dw,wmin,pi,spi,d character*1 lg integer*4 nw0,i,nc,nw parameter (nw0=2000) common/convc/ff(2*nw0+1),nc pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) if(nw.gt.nw0)call report('nw too big in dc') w=wmin-dw if(lg.eq.'L')then do 1 nc=1,nw w=w+dw dd=((w-wmin)/d)**2 ff(nc)=1.0d0/d/(dd+1.0d0)/pi 1 if(dd.gt.1.0d3)goto 2 else do 3 nc=1,nw w=w+dw dd=((w-wmin)/d)**2 ff(nc)=exp(-dd)/d/spi 3 if(dd.gt.20.0d0)goto 2 endif 2 do 4 i=1,nc-1 4 ff(nc+i)=ff(i+1) ff(nc)=ff(1) do 5 i=1,nc-1 5 ff(nc-i)=ff(nc+i) return end subroutine trabc(ia,xm,ym,zm,a,ss) integer*4 ia real*8 xm(*),ym(*),zm(*),a(*),ss xm(ia)=xm(ia)+a(1)*ss ym(ia)=ym(ia)+a(2)*ss zm(ia)=zm(ia)+a(3)*ss return end subroutine make3v(v1,v2,v3) c makes system of orthonormal vectors from v1 and v2 implicit none real*8 v1(*),v2(*),v3(*) v3(1)=v1(2)*v2(3)-v1(3)*v2(2) v3(2)=v1(3)*v2(1)-v1(1)*v2(3) v3(3)=v1(1)*v2(2)-v1(2)*v2(1) call norm(v1) call norm(v3) v2(1)=v3(2)*v1(3)-v3(3)*v1(2) v2(2)=v3(3)*v1(1)-v3(1)*v1(3) v2(3)=v3(1)*v1(2)-v3(2)*v1(1) call norm(v2) return end subroutine readopt1(K0,n,NSA,iprint,nmol,nnames,isr,natm, 1imk,al,dt,kelvin,bohr,pi,dbb,perx,pery,perz,tp,sol,dw,wmin, 1wmax,res,ac,bc,cc,lper,lcell,lp27,lp28,lp29,lproject,lnnames, 1lpair,LABS,LVCD,LRAM,LROA,ltti,ldress,ltei,lani,lconv,lg,snam, 1lincr,lcct,lexp,lexca,px2,py2,pz2,sc,dtau,nw,nw0,lcidr,lave,irep, 1natx) implicit none integer*4 K0,n,NSA,iprint,nmol,nnames,ia,isr(*),natm(*), 1imk(*),al(nmol,n),i,mm,nw,irep,nw0,natx real*8 dt,kelvin,bohr,pi,dbb,perx,pery,perz,tp,sol,dw,wmin,wmax, 1res,ac(*),bc(*),cc(*),px2,py2,pz2,sc,dtau logical lper,lcell,lp27,lp28,lp29,lproject,lnnames,lpair,LABS, 1LVCD,LRAM,LROA,ltti,ldress,ltei,lani,lconv,ldyn,lincr,lcct,lexp, 1lexca,lcidr,lave character*1 lg character*10 s10,snam(*) NSA=0 c write tensor average each irep geometry irep=10 iprint=100 nmol=0 natx=0 dt=0.00009672d0 kelvin=298.15d0 bohr=0.529177d0 pi=4.0d0*datan(1.0d0) c au^2 -> debye^2 dbb=2.541765d0**2 perx=0.0d0 pery=0.0d0 perz=0.0d0 tp=2.0d0*pi sol=3.0d10 c dw,wmin,wmax step and frequency limits in cm-1 dw=2.0d0 wmin=50.0d0 wmax=4000.0d0 res=5.0d0 lper=.false. lcell=.false. lp27=.false. lp28=.false. lp29=.false. lproject=.false. lnnames=.false. lpair=.false. LABS=.false. LVCD=.false. LRAM=.false. LROA=.false. ltti=.false. ldress=.false. ltei=.false. lani=.false. lconv=.false. lincr=.false. c transfer Cartesian tensors via CCT: lcct=.false. c read atomic polarizabilities from POLAR.TXT: lexp=.false. lexca=.false. c Rayleigh parameters: lcidr=.false. c Average polarizabilities lave=.false. lg='G' inquire(file='DYNAMIM.OPT',exist=ldyn) if(.not.ldyn)call report('DYNAMIM.OPT not found') open(10,file='DYNAMIM.OPT') 311 read(10,71,end=3111,err=3111)s10 71 format(a10) if(s10(1:4).eq.'NRES')then read(10,*)NSA read(10,*)(isr(ia),ia=1,NSA) endif if(s10(1:4).eq.'TSTE')read(10,*)dt if(s10(1:4).eq.'IPRI')read(10,*)iprint if(s10(1:2).eq.'DW')read(10,*)dw if(s10(1:3).eq.'RES')read(10,*)res if(s10(1:3).eq.'CCT')read(10,*)lcct if(s10(1:4).eq.'LEXP')read(10,*)lexp if(s10(1:4).eq.'LEXC')read(10,*)lexca if(s10(1:4).eq.'WMIN')read(10,*)wmin if(s10(1:4).eq.'WMAX')read(10,*)wmax if(s10(1:4).eq.'TTTI')read(10,*)ltti if(s10(1:4).eq.'TEMK')read(10,*)kelvin if(s10(1:4).eq.'DRES')read(10,*)ldress if(s10(1:4).eq.'TENI')read(10,*)ltei if(s10(1:4).eq.'ANIM')read(10,*)lani if(s10(1:4).eq.'CONV')read(10,*)lconv if(s10(1:4).eq.'INCR')read(10,*)lincr if(s10(1:5).eq.'LCIDR')read(10,*)lcidr if(s10(1:4).eq.'IREP')read(10,*)irep if(s10(1:4).eq.'LAVE')read(10,*)lave if(s10(1:4).eq.'CONT')read(10,*)lg if(s10(1:4).eq.'PROJ')read(10,*)lproject if(s10(1:4).eq.'PERI')read(10,*)lper,perx,pery,perz if(s10(1:4).eq.'CELL')then read(10,*)lcell do 72 i=1,3 72 read(10,*)ac(i),bc(i),cc(i) endif if(s10(1:3).eq.'P27')read(10,*)lp27 if(s10(1:3).eq.'P28')read(10,*)lp28 if(s10(1:3).eq.'P29')read(10,*)lp29 if(s10(1:4).eq.'PAIR')read(10,*)lpair if(s10(1:4).eq.'NMOL')then natx=0 read(10,*)nmol do 191 mm=1,nmol read(10,*)natm(mm),(al(mm,ia),ia=1,natm(mm)) 191 if(natm(mm).gt.natx)natx=natm(mm) endif c c named tensor files: names and molecular assignment imk must be given: if(s10(1:6).eq.'NNAMES')then read(10,*)lnnames if(lnnames)then read(10,*)nnames if(nnames.gt.K0)call report('too many molecular kinds') do 1911 i=1,nnames 1911 read(10,71)snam(i) read(10,*)nmol do 160 mm=1,nmol 160 read(10,*)imk(mm) endif endif goto 311 c NSA - number of atoms for spectral response (followed by their list) c TSTE - time step in picoseconds c lper ... correct for periodic conditions given by perx pery perz (A) c lcell ... non-orthogonal crystall cell c lp27 ... calculate with the 3^3 large periodic box c lp28 ... damp towards walls of the smaller box c lp29 ... real superbox c lincr ... restart accumulation of spectra from the disk 3111 close(10) px2=perx/2.0d0 py2=pery/2.0d0 pz2=perz/2.0d0 nw=int((wmax-wmin)/dw) if(nw.gt.nw0)then write(6,6000)nw,nw0 6000 format(i6,' > ',i6,':too many spectral points') stop endif sc=1.0d0 dtau=dt*1.0d-12/2.418D-17 return end subroutine readopt3(nmol,lpair) integer*4 nmol logical ldyn,lpair character*10 s10 lpair=.false. nmol=0 inquire(file='DYNAMIM.OPT',exist=ldyn) if(.not.ldyn)call report('DYNAMIM.OPT not found') open(10,file='DYNAMIM.OPT') 311 read(10,71,end=3111,err=3111)s10 71 format(a10) if(s10(1:4).eq.'NMOL')read(10,*)nmol if(s10(1:4).eq.'PAIR')read(10,*)lpair goto 311 3111 close(10) return end SUBROUTINE WRITETTT(MX,NAT,ALPHA,A,G,s) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3) character*(*)s OPEN(2,FILE=s) WRITE(2,2000)NAT 2000 FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 WRITE(2,2001)L,IX,(ALPHA(IIND,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3F18.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx jy jz') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 WRITE(2,2001)L,IX,(G(IIND,I,J),J=1,3) WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3) CLOSE(2) RETURN END subroutine wrt2(IA,MX,ALPHA,s) implicit none integer*4 MX,IIND,IA,IX,I,J real*8 ALPHA(3*MX,3,3) character*(*) s write(6,*)s do 611 IX=1,3 IIND=3*(IA-1)+IX 611 write(6,3007)((ALPHA(IIND,I,J),I=1,3),J=1,3) 3007 format(9f8.3) return end subroutine wrt3(IA,MX,ALPHA,s) implicit none integer*4 MX,IIND,IA,IX,I,J,K real*8 ALPHA(3*MX,3,3,3) character*(*) s write(6,*)s do 611 IX=1,3 IIND=3*(IA-1)+IX 611 write(6,3007)(((ALPHA(IIND,I,J,K),I=1,3),J=1,3),K=1,3) 3007 format(9f8.3) return end SUBROUTINE TTT2(ALPHA,A,G,ICO,c,iwr) c Tensors, not derivatives C Transforms A and G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C c supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,IA,IB,IG,ID,IC,iwr real*8 ALPHA(3,3),G(3,3),A(3,3,3),c(3),SIGN,SUM,EPS C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*c(IG)*ALPHA(IA,ID) G(IA,IB)=G(IA,IB)+SUM IF(ICO.EQ.3)G(IA,IB)=0.0d0 2 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(6,*)' G transformed into overlap origin' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' endif c A .. last index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM= 1c(1)*ALPHA(IA,1)+c(2)*ALPHA(IA,2)+c(3)*ALPHA(IA,3) A(IC,IB,IA)=A(IC,IB,IA)+ 1SIGN*(SUM-1.5d0*(c(IB)*ALPHA(IA,IC)+c(IC)*ALPHA(IA,IB))) IF(ICO.EQ.3)A(IC,IB,IA)=0.0d0 3 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(6,*)' A transformed into overlap origin' IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A set to zero' endif RETURN END c =============================== subroutine readcharge(n,g) implicit none integer*4 n,i real*8 g(*) open(99,file='CHARGE.LST',status='old') read(99,*) do 19 i=1,n 19 read(99,*)g(i) close(99) return end subroutine readatpolars(n,pol) implicit none integer*4 n,ia real*8 pol(*) open(99,file='POLAR.LST',status='old') read(99,*) do 18 ia=1,n 18 read(99,*)pol(ia) close(99) return end SUBROUTINE TTT1(ALPHA,A,G,ICO,c,iwr) c Tensors, not derivatives C Transforms A and G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C c supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,IA,IB,IG,ID,IC,iwr real*8 ALPHA(3,3),G(3,3),A(3,3,3),c(3),SIGN,SUM,EPS C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*c(IG)*ALPHA(IA,ID) G(IA,IB)=G(IA,IB)+SUM IF(ICO.EQ.3)G(IA,IB)=0.0d0 2 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(6,*)' G transformed into overlap origin' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' endif c A .. last index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM= 1c(1)*ALPHA(IA,1)+c(2)*ALPHA(IA,2)+c(3)*ALPHA(IA,3) A(IC,IB,IA)=A(IC,IB,IA)+ 1SIGN*(SUM-1.5d0*(c(IB)*ALPHA(IA,IC)+c(IC)*ALPHA(IA,IB))) IF(ICO.EQ.3)A(IC,IB,IA)=0.0d0 3 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(6,*)' A transformed into overlap origin' IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A set to zero' endif RETURN END SUBROUTINE CRMSo(NAT,INDBIG,IOX,MX,IPO,LNAMES,natname,RS,RB,IND, 2rname,NATS,LINV,MXB,A,XS,XB,natx) c make the rotation matrix for all atoms in the overlap IPO implicit none integer*4 IOX,MXB,MX,natx integer*4 II,NAT,INDBIG(IOX,MX),NATS,natname(*),IX,IND(*),I,IERR, 2IPO,J,ITU real*8 RS(3,MX),RB(3,MX),TB(3),TS(3),rname(IOX,3,natx), 1A(3,3),XS(MXB,3),XB(MXB,3) logical LNAMES,LINV c if(LNAMES)NATS=natname(IPO) ITU=0 DO 27 II=1,NAT if(INDBIG(IPO,II).NE.0.and.ITU.lt.MXB)then C Include atom II into the transformation set: ITU=ITU+1 IF(ITU.EQ.MXB)write(6,*)'Warning: MXB reached in CRMSo' DO 28 IX=1,3 if(LNAMES)then XS(ITU,IX)=rname(IPO,IX,ABS(INDBIG(IPO,II))) else XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II))) endif 28 XB(ITU,IX)=RB(IX,II) IND(ITU)=II endif 27 CONTINUE c if(ITU.LT.3) CALL REPORT(' ITU < 3 !') c C Calculate center of the coordinates that will be rotated: DO 33 I=1,3 TB(I) =0.0d0 33 TS(I)=0.0d0 DO 34 II=1,ITU DO 34 J=1,3 TB(J)=TB(J)+XB(II,J) 34 TS(J)=TS(J)+XS(II,J) DO 38 J=1,3 TB(J)=TB(J)/DBLE(ITU) 38 TS(J)=TS(J)/DBLE(ITU) C Subtract mass-centers: DO 35 II=1,ITU DO 35 J=1,3 XS(II,J)=XS(II,J)-TS(J) 35 XB(II,J)=XB(II,J)-TB(J) C IF(LINV)THEN DO 39 II=1,ITU DO 39 J=1,3 39 XS(II,J)=-XS(II,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR,A,XB,XS,mxb,ITU) IF(IERR.EQ.1)CALL REPORT('CRMs DoMatrix error') return end subroutine ct(ta,t,o) real*8 ta(3,3),t(3,3),o(3,3) integer*4 i,j do 1 i=1,3 do 1 j=1,3 1 ta(i,j)=t(i,1)*o(1,j)+t(i,2)*o(2,j)+t(i,3)*o(3,j) return end c subroutine ctd(IIND,ALPHA,tem,M) integer*4 IIND,M,i,j real*8 ALPHA(M,3,3),tem(3,3) do 1 i=1,3 do 1 j=1,3 1 tem(i,j)=ALPHA(IIND,i,j) return end c subroutine ctdd(IIND,A,tem,M) integer*4 IIND,M,i,j,k real*8 A(M,3,3,3),tem(3,3,3) do 1 i=1,3 do 1 j=1,3 do 1 k=1,3 1 tem(i,j,k)=A(IIND,i,j,k) return end subroutine mkttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3),rt2,rt5 integer*4 ix,jx,kx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 ttt(ix,jx,kx)=-15.0d0*rt(ix)*rt(jx)*rt(kx)/rt5/rt2 if(ix.eq.jx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(kx)/rt5 if(ix.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(jx)/rt5 613 if(jx.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(ix)/rt5 return end c subroutine mktt(rt,tt) implicit none real*8 rt(3),tt(3,3),rt2,rt5 integer*4 ix,jx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 611 ix=1,3 do 611 jx=1,3 611 tt(ix,jx)=3.0d0*rt(ix)*rt(jx)/rt5 do 610 ix=1,3 610 tt(ix,ix)=tt(ix,ix)-rt2/rt5 return end subroutine tadd(bigpol,tc) implicit none c B=B+C real*8 bigpol(3,3),tc(3,3) integer*4 ix,jx do 615 ix=1,3 do 615 jx=1,3 615 bigpol(ix,jx)=bigpol(ix,jx)+tc(ix,jx) return end subroutine ctrl3(iprint,dt,dtau,n,NSA,dw,wmin,wmax,nw,res,lper, 1perx,pery,perz,lnnames,lincr,lp27,lp28,lp29,lcell,LRAM,LROA, 1LABS,LVCD,lani,lconv,ltei,ltti,lcct,lexp,lexca,lproject,ldress,lg) implicit none integer*4 iprint,n,NSA,nw real*8 dt,dtau,dw,wmin,wmax,res,perx,pery,perz logical lper,lnnames,lincr,lp27,lp28,lp29,lcell, 1LRAM,LROA,LABS,LVCD, 1lani,lconv,ltei, 1ltti,lcct,lexp,lexca,lproject,ldress character*1 lg write(6,6013)iprint,dt,dtau,n,NSA,dw,wmin,wmax,nw,res,lper,perx, 1pery,perz,lnnames,lincr, 1lp27,lp28,lp29,lcell, 1LRAM,LROA,LABS,LVCD, 1lani,lconv,ltei, 1ltti,lcct,lexp,lexca,lproject,ldress,lg 6013 format(/,' Parameters',/,' **********',/, 1' printing interval : ',i10,' steps',/, 1' dt : ',f10.6,' ps (',f10.4,' au)',/, 1' n : ',i10,' atoms',/, 1' NSA : ',i10,' atoms',/, 1' DwA : ',f10.2,' cm-1',/, 1' wmin : ',f10.2,' cm-1',/, 1' wmax : ',f10.2,' cm-1',/, 1' Num. of freqs. : ',i10,/, 1' Convol. broadening: ',f10.2,' cm-1'/, 1' Per. box : ',l1,3f10.2,' A',/, 1' lnnames : ',l1,/, 1' lincr : ',l1,/, 1' lp27 lp28 lp29 lcell: ',4l3,/, 1' RAM ROA ABS VCD : ',4l3,/, 1' animated : ',l1,/, 1' convolution : ',l1,/, 1' internal ten : ',l1,/, 1' internal ttt : ',l1,/, 1' Cartesian transfer: ',l1,/, 1' at pol POLAR.TXT : ',l1,/, 1' charges CHARGE.LST: ',l1,/, 1' project pol : ',l1,/, 1' dress pol : ',l1,/, 1' convolution type : ',a1,/,/) return end subroutine inttt(lnnames,snam,NAT1,n,ALPHAI,AAI,GI,iii, 1K0,IIM,ALPHAM,GM,AAM,POLI,POLGI,POLAI,POLM,POLGM,POLAM, 1nnames,fr) implicit none logical lnnames integer*4 mm,istart,K0,iend,L,I,n,NAT1,iii(4,n),IIM(K0,4,n), 1IIND,IX,nnames,IP,JP,J,K character*10 snam(K0) real*8 ALPHAI(3*n,3,3),GI(3*n,3,3),AAI(3*n,3,3,3),fr, 1ALPHAM(9,K0,3*n),GM(9,K0,3*n),POLI(3,3),POLGI(3,3), 1POLAI(3,3,3),POLM(K0,9),POLGM(K0,9),POLAM(K0,27),AAM(27,K0,3*n) if(lnnames)then do 186 mm=1,nnames do 1851 istart=1,len(snam(mm)) 1851 if(snam(mm)(istart:istart).ne.' ')goto 1852 1852 do 1853 iend=len(snam(mm)),1,-1 1853 if(snam(mm)(iend:iend).ne.' ')goto 1854 1854 open(15, 1 FILE=snam(mm)(istart:iend)//'.int.ttt',status='old') write(6,7006)mm,' Raman',(snam(mm)(i:i),i=istart,iend) 7006 format(' Kind ',i4,', ',a6,' tensors from ',40a1) CALL READTTT(n,ALPHAI,AAI,GI,NAT1) if(NAT1.gt.2)call readatms(15,NAT1,iii,n) close(15) do 18541 L=1,NAT1 do 18542 I=1,4 18542 IIM(mm,I,L)=iii(I,L) do 18541 IX=1,3 IIND=3*(L-1)+IX do 18541 I=1,3 do 18541 J=1,3 ALPHAM(I+(J-1)*3,mm,IIND)=ALPHAI(IIND,I,J) GM(I+(J-1)*3,mm,IIND)=GI(IIND,I,J) do 18541 K=1,3 18541 AAM(I+(J-1)*3+(K-1)*9,mm,IIND)=AAI(IIND,I,J,K) CALL readpolal(POLI,snam(mm)(istart:iend)//'.int.pol',fr) CALL readpolg( POLGI,snam(mm)(istart:iend)//'g.int.pol') CALL readpola( POLAI,snam(mm)(istart:iend)//'a.int.pol') do 186 IP=1,3 do 186 JP=1,3 POLAM(mm,IP+(JP-1)*3 )=POLAI(IP,JP,1) POLAM(mm,IP+(JP-1)*3+ 9)=POLAI(IP,JP,2) POLAM(mm,IP+(JP-1)*3+18)=POLAI(IP,JP,3) POLGM(mm,IP+(JP-1)*3)=POLGI(IP,JP) 186 POLM(mm,IP+(JP-1)*3)=POLI(IP,JP) else write(6,7007) 7007 format(' All molecules: INT.TTT tensors') open(15,FILE='INT.TTT',status='old') CALL READTTT(n,ALPHAI,AAI,GI,NAT1) if(NAT1.gt.2)call readatms(15,NAT1,iii,n) close(15) CALL readpolal(POLI,'INT.POL',fr) CALL readpolg(POLGI,'INTG.POL') CALL readpola(POLAI,'INTA.POL') endif return end subroutine intte(lnnames,nnames,snam,n,APTI,AI,NAT1,iii, 1K0,IIM,APTM,AM) implicit none integer*4 nnames,n,NAT1,iii(4,n),K0,IIM(K0,4,n),mm,istart,iend, 1I,L,J logical lnnames character*10 snam(K0) real*8 APTI(n,3,3),AI(n,3,3),APTM(n,9,K0),AM(n,9,K0) if(lnnames)then do 189 mm=1,nnames do 1891 istart=1,len(snam(mm)) 1891 if(snam(mm)(istart:istart).ne.' ')goto 1892 1892 do 1893 iend=len(snam(mm)),1,-1 1893 if(snam(mm)(iend:iend).ne.' ')goto 1894 1894 open(15, 1 FILE=snam(mm)(istart:iend)//'.int.ten',status='old') write(6,7006)mm,'dipole',(snam(mm)(i:i),i=istart,iend) 7006 format(' Kind ',i4,', ',a6,' tensors from ',40a1) CALL READTENI(n,APTI,AI,NAT1) write(6,*)'NAT1',NAT1,' n:',n if(NAT1.gt.2)call readatms(15,NAT1,iii,n) CLOSE(15) do 189 L=1,NAT1 do 18842 I=1,4 18842 IIM(mm,I,L)=iii(I,L) do 189 I=1,3 do 189 J=1,3 APTM(L,I+(J-1)*3,mm)=APTI(L,I,J) 189 AM( L,I+(J-1)*3,mm)=AI(L,I,J) else OPEN(15,FILE='INT.TEN',status='old') CALL READTENI(n,APTI,AI,NAT1) if(NAT1.gt.2)call readatms(15,NAT1,iii,n) CLOSE(15) write(6,7008) 7008 format(' All molecules: INT.TEN tensors') endif return end subroutine wrave(pkind,a1,nmol,istep,pol,gpol,apol,icw) c ic=1 write summary c ic=2 progress errors implicit none integer*4 pkind,ic,I,istep,nmol,J,ityp,id,ii,K,icw,JK real*8 a1(pkind),pol(3,3),sp,ab,aa, 1s1,s2,s3,gpol(3,3),apol(3,3,3),Gaa,Gab,sG,s1g,s3g,s2g,s1a,s3a, 1xyxy,zxzx,yzyz,xxyy,zzxx,yyzz,xxyxz,yyzyx,zzxzy,xxxyz, 1yyyzx,zzzxy,xxxx,xyxxz,yzyyx,zxzzy,S0,S31,S32,EXCA,gpisvejc, 1AMU,BOHR,T0,T3,zxxyx,yzzxz,S0e,S31e,S32e,yyyy,zzzz,xzxz, 1xyyzy,su(3,3),st(3),sw(3,3,3),T31,T32,T33,err,rerr real*8,allocatable::alave(:) character*1 xyz(3) data xyz/'X','Y','Z'/ character*50,allocatable:: ss(:) c rewrite sums and average: allocate(alave(pkind)) do 800 i=1,pkind 800 alave(i)=a1(i)/dble(nmol*istep) EXCA=5320.0d0 AMU=1822.0d0 BOHR=0.529177d0 gpisvejc=(AMU*BOHR**5)*2.0d0*4.0d0*atan(1.0d0)/EXCA c average molecular polarizabilities laboratory frame c equilibrium and derivs ityp=15 allocate(ss(ityp)) ss( 1)=' alpha_AA x alpha_AA ' ss( 2)=' alpha_AA x G_AA ' ss( 3)=' alpha_AB ' ss( 4)=' alpha_AB x alpha_AB ' ss( 5)=' alpha_AA x alpha_BB ' ss( 6)=' G_AB ' ss( 7)=' alpha_AB x G_AB ' ss( 8)=' alpha_AA x G_BB ' ss( 9)=' alpha_AA x A_B,AC ' ss(10)=' alpha_AA x A_A,BC ' ss(11)=' alpha_AB x A_A,AC ' ss(12)=' ai_AB x aj_AB ' ss(13)=' 2 aj_AB Gi_AB +ai_AA(G_AA-G_BB) ' ss(14)=' Rij_C ai_AA aj_AB ' ss(15)=' [ai_AA (Aj_B,AC+Aj_A,BC)-2 aj_AB Ai_A,AC] / 3 ' aa= pol(1,1)+ pol(2,2)+ pol(3,3) Gaa=gpol(1,1)+gpol(2,2)+gpol(3,3) sp=aa/3.0d0 sG=Gaa/3.0d0 ab =pol(1,1)* pol(1,1)+pol(1,2)* pol(1,2)+pol(1,3)* pol(1,3) 1 +pol(2,1)* pol(2,1)+pol(2,2)* pol(2,2)+pol(2,3)* pol(2,3) 1 +pol(3,1)* pol(3,1)+pol(3,2)* pol(3,2)+pol(3,3)* pol(3,3) Gab=pol(1,1)*Gpol(1,1)+pol(1,2)*Gpol(1,2)+pol(1,3)*Gpol(1,3) 1 +pol(2,1)*Gpol(2,1)+pol(2,2)*Gpol(2,2)+pol(2,3)*Gpol(2,3) 1 +pol(3,1)*Gpol(3,1)+pol(3,2)*Gpol(3,2)+pol(3,3)*Gpol(3,3) c = (3ab ab - aa bb)/30 c = (2aa bb - ab ab)/15 c = ( aa bb + 2 ab ab)/15 s1 =(3.0d0 * ab- aa *aa)/30.0d0 s1g=(3.0d0 *Gab- aa*Gaa)/30.0d0 s2 =(2.0d0 * aa* aa- ab)/15.0d0 s2g=(2.0d0 * aa*Gaa- Gab)/15.0d0 s3 =(2.0d0* aa* aa +4.0d0* ab)/30.0d0 s3g=( aa*Gaa +2.0d0*Gab)/15.0d0 c EPS_abc a_cd A_bda / 15: s1a=0.0d0 do 14 id=1,3 14 s1a=s1a+(pol(3,id)*apol(2,id,1) -pol(2,id)*apol(3,id,1) 1 +pol(1,id)*apol(3,id,2) -pol(3,id)*apol(1,id,2) 2 +pol(2,id)*apol(1,id,3) -pol(1,id)*apol(2,id,3))/15.0d0 s3a=-s1a/2.0d0 c AA B,AC = EPS_abg gd a,bd /15 c AA A,BC = 0 c AB A,AC = -EPS_abg gd a,bd /30 if(icw.eq.2)open(60,file='PAVE.TXT') if(icw.eq.2)write(6,6004) 6004 format(/,' Equilibrium polarizability invariants:') c 1 1 .. 3 a_AA x a_AA c 2 4 .. 6 a_AA x G_AA c 3 7 .. 15 a_AB c 4 16 .. 24 a_AB x a_AB c 5 25 .. 33 a_AA x a_BB c 6 34 .. 42 G_AB c 7 43 .. 51 a_AB x G_AB c 8 52 .. 60 a_AA x G_BB c 9 61 .. 87 a_AA x A_B,AC c 10 88 ..114 a_AA x A_A,BC c 11 115 ..141 a_AB x A_A,AC c Two-molecule averages c 12 142 ..150 ai_AB x aj_AB c 13 151 ..159 2ai_AB x Gj_AB + ai_AA x (Gj_AA - Gj_BB) c 14 160 ..186 2 R_ijC ai_AA x aj_AB c 15 187 ..213 (ai_AA Aj_B,AC + ai_AA Aj_A,BC - 2 aj_AB Ai_A,AC )/3 do 1 ic=1,ityp if(icw.eq.1)then write(6 ,6001)ss(ic) ,' Ideal:' write(60,6001)ss(ic) ,' Ideal:' 6001 format(/,A50,A7) write(6,6002) write(60,6002) 6002 format(1x,2(9x,'X',9x,'Y',9X,'Z',3x)) else write(6 ,60011)ss(ic) ,' Error:' 60011 format(A50,A7,$) endif if(ic.eq.1.or.ic.eq.2)then if(ic.eq.1)then do 201 J=1,3 201 st(J)=s3 ii=0 else do 202 J=1,3 202 st(J)=s3g ii=3 endif call ferr(err,rerr,alave,st,ii) if(icw.eq.1)then write(60,6007)(alave(J+ii),J=1,3),(st(J),J=1,3) write(6 ,6007)(alave(J+ii),J=1,3),(st(J),J=1,3) 6007 format(2x,2(1x,3e11.3)) endif endif if((ic.ge.3.and.ic.le.8).or.(ic.ge.12.and.ic.le.13))then c analytical elements: call ssu(nmol,su,ic,sp,s1,s2,s3,sG,s1g,s2g,s3g) ii=6+(ic-3)*9 if(ic.eq.12.or.ic.eq.13)ii=141+(ic-12)*9 call serr(err,rerr,alave,su,ii) if(icw.eq.1)then do 11 J=1,3 write(6, 631)xyz(J),(alave(ii+I+3*(J-1)),I=1,3),(su(I,J),I=1,3) 631 format(1x,a1,1x,3e11.3,1x,3e11.3) 11 write(60,631)xyz(J),(alave(ii+I+3*(J-1)),I=1,3),(su(I,J),I=1,3) endif endif if((ic.ge.9.and.ic.le.11).or.(ic.ge.14.and.ic.le.15))then c analytical elements: call ssw(s1a,s3a,sw,ic) ii=60+27*(ic-9) if(ic.eq.14.or.ic.eq.15)ii=159+27*(ic-14) call terr(err,rerr,alave,sw,ii) if(icw.eq.1)then do 12 K=1,3 write(6,6014)K 6014 format(' K = ',i1,':') do 12 J=1,3 JK=ii+3*(J-1)+9*(K-1) write(6 ,631)xyz(I),(alave(JK+I),I=1,3),(sw(I,J,K),I=1,3) 12 write(60,631)xyz(I),(alave(JK+I),I=1,3),(sw(I,J,K),I=1,3) endif endif if(icw.eq.1)write(60,60071) 60071 format(' Error/relative:',$) write(6,60072)err,rerr 60072 format(2e11.3) 1 continue c SCP backscattering: xx xx + yx yx c to better average, use "xx xx"=(xx xx + yy yy + zz zz)/3 c "yx yx"=(xy xy + xz xz + yz yz)/3 S0=(alave(1)+alave(2)+alave(3)+alave(15+1+3*(2-1)) 1+alave(15+1+3*(3-1))+alave(15+2+3*(3-1)))/3.0d0 S0e=s1+s3 c SCP backscattering ROA: 2 xy xy -xx yy + xx xx c +(w/3)(xx y,xz + xx x,yz - yx x,xz) xyxy=alave(42+1+3*(2-1)) yzyz=alave(42+2+3*(3-1)) zxzx=alave(42+3+3*(1-1)) xyxy=(xyxy+zxzx+yzyz)/3.0d0 xxyy=alave(51+1+3*(2-1)) zzxx=alave(51+3+3*(1-1)) yyzz=alave(51+2+3*(3-1)) xxyy=(xxyy+zzxx+yyzz)/3.0d0 xxxx=(alave(4)+alave(5)+alave(6))/3.0d0 S31=2.0d0*xyxy-xxyy+xxxx S31e=2.0d0/15.0d0*(3.0d0*Gab-aa*Gaa) xxyxz=alave(60+1+3*(2-1)+9*(3-1)) yyzyx=alave(60+2+3*(3-1)+9*(1-1)) zzxzy=alave(60+3+3*(1-1)+9*(2-1)) xxyxz=(xxyxz+yyzyx+zzxzy)/3.0d0 xxxyz=alave(87+1+3*(2-1)+9*(3-1)) yyyzx=alave(87+2+3*(3-1)+9*(1-1)) zzzxy=alave(87+3+3*(1-1)+9*(2-1)) xxxyz=(xxxyz+yyyzx+zzzxy)/3.0d0 xyxxz=alave(114+1+3*(2-1)+9*(3-1)) yzyyx=alave(114+2+3*(3-1)+9*(1-1)) zxzzy=alave(114+3+3*(1-1)+9*(2-1)) xyxxz=(xyxxz+yzyyx+zxzzy)/3.0d0 S32=(xxyxz+xxxyz-2.0d0*xyxxz)/3.0d0 S32e=2.0d0*s1a/3.0d0 S0 =S0 *180.0d0*AMU*BOHR**4 S0e =S0e*180.0d0*AMU*BOHR**4 S31 =S31 *gpisvejc*360.0d0 S31e=S31e*gpisvejc*360.0d0 S32 =S32 *gpisvejc*360.0d0 S32e=S32e*gpisvejc*360.0d0 write(6,6003)S0,S0e,S31,S31e,S32,S32e,S31+S32,S31e+S32e 6003 format(' RamSCP180 :',e11.3,', ideal:',e11.3,/,/, 1 ' ROASCP180G:',e11.3,', ',e11.3,/, 1 ' A:',e11.3,', ',e11.3,/, 1 ' ',11(1h-) ,'--------',11(1h-),/, 1 ' T:',e11.3,', ',e11.3,/) c 1 1 .. 3 a_AA x a_AA c 2 4 .. 6 a_AA x G_AA c 3 7 .. 15 a_AB c 4 16 .. 24 a_AB x a_AB c 5 25 .. 33 a_AA x a_BB c 6 34 .. 42 G_AB c 7 43 .. 51 a_AB x G_AB c 8 52 .. 60 a_AA x G_BB c 9 61 .. 87 a_AA x A_B,AC c 10 88 ..114 a_AA x A_A,BC c 11 115 ..141 a_AB x A_A,AC c 12 142 ..150 ai_AB x aj_AB c 13 151 ..159 2ai_AB x Gj_AB + ai_AA x (Gj_AA - Gj_BB) c 14 160 ..186 R_Cj ai_AA x aj_AB c 15 187 ..213 (ai_AA Aj_B,AC + ai_AA Aj_A,BC - 2 aj_AB Ai_A,AC )/3 c Two molecules: c Raman: xxxx=alave(141+1+3*(1-1)) yyyy=alave(141+2+3*(2-1)) zzzz=alave(141+3+3*(3-1)) xxxx=(xxxx+yyyy+zzzz)/3.0d0 xyxy=alave(141+1+3*(2-1)) xzxz=alave(141+1+3*(3-1)) yzyz=alave(141+2+3*(3-1)) xyxy=(xyxy+xzxz+yzyz)/3.0d0 T0=(XXXX+xyxy-sp**2*dble(nmol-1))*180.0d0*AMU*BOHR**4 c ROA: c 2 Zij ai_XX x aj_XY: zxxyx=alave(159+1+3*(2-1)+9*(3-1)) xyyzy=alave(159+2+3*(3-1)+9*(1-1)) yzzxz=alave(159+3+3*(1-1)+9*(2-1)) zxxyx=(zxxyx+xyyzy+yzzxz)/3.0d0 T31=2.0d0*zxxyx*gpisvejc*360.0d0 c 2ai_XY x Gj_XY + ai_XX x (Gj_XX - Gj_YY): xyxy=alave(150+1+3*(2-1)) yzyz=alave(150+2+3*(3-1)) zxzx=alave(150+3+3*(1-1)) xyxy=(xyxy+yzyz+zxzx)/3.0d0 T32=xyxy*gpisvejc*360.0d0 c (ai_XX Aj_Y,XZ + ai_XX Aj_X,YZ - 2 aj_XY Ai_X,XZ )/3 xxyxz=alave(186+1+3*(2-1)+9*(3-1)) yyzyx=alave(186+2+3*(3-1)+9*(1-1)) zzxzy=alave(186+3+3*(1-1)+9*(2-1)) xxyxz=(xxyxz+yyzyx+zzxzy)/3.0d0 T33=xxyxz/3.0d0*gpisvejc*360.0d0 T3 =T31+T32+T33 write(6,6008)T0,T31,T32,T33,T3 6008 format(' Two-molecule contribution:',/, 1 ' T-RamSCP180 :',e11.3,/,/, 1 ' T-ROASCP180a:',e11.3,/, 1 ' T-ROASCP180G:',e11.3,/, 1 ' T-ROASCP180A:',e11.3,/, 1 ' ',11(1h-) ,'--------',11(1h-),/, 1 ' T-ROASCP180T:',e11.3,/) close(60) return end subroutine inpder(ndim,im,l,g,a,e) c load normal mode derivatives of alpha, Gp and A tensors) implicit none integer*4 nm,im,II,i,j,k,i9,ndim real*8 l(*),g(*),a(*),ecm,e(*) logical lex inquire(file='SMALL.Q.TTT',exist=lex) if(lex)then open(22,file='SMALL.Q.TTT') read(22,*) read(22,*)nm write(6,*)nm,' modes in SMALL.Q.TTT' if(nm.gt.ndim)call report('too many normal modes') read(22,*) read(22,*) im=0 DO 1 II=1,nm read(22,*)ecm,ecm if(dabs(ecm).gt.1.0d-1)then im=im+1 e(im)=ecm DO 11 i=1,3 11 read(22,*)l(i+9*(im-1)),(l(i+3*(j-1)+9*(im-1)),j=1,3) endif 1 continue read(22,*) read(22,*) im=0 DO 2 II=1,nm read(22,*)ecm,ecm if(dabs(ecm).gt.1.0d-1)then im=im+1 DO 21 i=1,3 21 read(22,*)g(i+9*(im-1)),(g(i+3*(j-1)+9*(im-1)),j=1,3) c j index magnetic endif 2 continue read(22,*) read(22,*) im=0 DO 3 II=1,nm read(22,*)ecm,ecm if(dabs(ecm).gt.1.0d-1)then im=im+1 i9=27*(im-1) DO 31 I=1,3 DO 31 J=1,3 31 read(22,*)A(I+3*(J-1)+i9),A(I+3*(J-1)+i9), 1 (A(I+3*(J-1)+9*(K-1)+i9),K=1,3) endif 3 continue close(22) write(6,*)im,' energies positive' else call report('SMALL.Q.TTT not found') endif return end subroutine setactive(NSA,luse,isr,n) implicit none integer*4 NSA,ia,isr(*),n,i,ic logical luse(*) if(NSA.ne.0)then do 315 ia=1,NSA 315 luse(isr(ia))=.true. else do 316 ia=1,n 316 luse(ia)=.true. endif write(6,*)'Active atoms:' ic=0 do 1 i=1,n 1 if(luse(i))ic=ic+1 if(ic.eq.n)then write(6,*)'(all)' else do 2 i=1,n if(luse(i))then write(6,6012)'X' else write(6,6012)'O' 6012 format(a1,$) endif 2 if(mod(i,60).eq.0)write(6,*) endif write(6,*) return end function DE(I,J) implicit none integer*4 I,J real*8 DE if(I.EQ.J)then DE=1.0d0 else DE=0.0d0 endif return end subroutine mkcorrect(n,NO,coo,opol,gopol,aopol,lpola,lpolg, 1 bigpol,biggpol,bigapol) c correction of polarizability for mutual interaction of groups implicit none integer*4 IO,IOP,NO,n,i,j,ix,jx,kx,id,ie real*8 rt(3),coo(n,3),rt2,tt(3,3),ttt(3,3,3),tem(3,3),ta(3,3), 1opol(n,3,3),gopol(n,3,3),aopol(n,3,3,3),tc(3,3),EPS, 1bigpol(3,3),ata(3,3),tai(3,3),taj(3,3),pi,sf,bohr,atai(3,3), 1ataj(3,3),tv(3),su,bigapol(3,3,3),biggpol(3,3) logical lpola,lpolg bohr=0.529177d0 do 608 IO=1,NO do 608 IOP=1,NO if(IOP.ne.IO)then do 609 ix=1,3 609 rt(ix)=coo(IO,ix)-coo(IOP,ix) rt2=rt(1)**2+rt(2)**2+rt(3)**2 if(rt2.gt.0.0d0)then c c T-tensor call mktt(rt,tt) c c TTT-tensor call mkttt(rt,ttt) c c aTa to a: call ctd(IO,opol,tem,n) call ct(ta,tt,tem) call ctd(IOP,opol,tem,n) call ct(tc,tem,ta) do 612 ix=1,3 do 612 jx=1,3 612 bigpol(ix,jx)=bigpol(ix,jx)+tc(ix,jx) c ATa-aTA to a: if(lpola)then do 614 i=1,3 do 614 j=1,3 tc(i,j)=0.0d0 do 614 ix=1,3 do 614 jx=1,3 do 614 kx=1,3 614 tc(i,j)=tc(i,j)+ 1 (aopol(IO ,ix,jx,i)*ttt(ix,jx,kx)*opol(IOP,kx,j)- 1 aopol(IOP,ix,jx,j)*ttt(ix,jx,kx)*opol(IO ,kx,i))/3.0d0 call tadd(bigpol,tc) endif c c GTG to alpha: if(lpolg)then pi=4.0d0*atan(1.0d0) sf=(2.0d0*pi/(532.0d0*10.0d0/bohr))**2 do 616 i=1,3 do 616 j=1,3 tc(i,j)=0.0d0 do 616 ix=1,3 do 616 jx=1,3 616 tc(i,j)=tc(i,j)+ 1 gopol(IO ,i,jx)*tt(ix,jx)*gopol(IOP,j,jx)*sf call tadd(bigpol,tc) endif c c contribution of alpha...alpha coupling to G: call ctd(IOP,opol,tem,n) call ct(taj,tt,tem) call ctd(IO,opol,tem,n) call ct(ata,tem,taj) do 620 ix=1,3 do 620 jx=1,3 tc(ix,jx)=0.0d0 do 620 i=1,3 do 620 j=1,3 620 tc(ix,jx)=tc(ix,jx)-EPS(jx,i,j)/2.0d0*coo(IO,i)*ata(j,ix) call tadd(biggpol,tc) c c contribution of alpha...A coupling to G if(lpola)then do 624 ix=1,3 do 624 jx=1,3 tai(ix,jx)=0.0d0 taj(ix,jx)=0.0d0 do 624 i=1,3 do 624 j=1,3 tai(ix,jx)=tai(ix,jx)+ttt(ix,i,j)*aopol(IO ,i,j,jx) 624 taj(ix,jx)=taj(ix,jx)+ttt(ix,i,j)*aopol(IOP,i,j,jx) do 625 ix=1,3 do 625 jx=1,3 atai(ix,jx)=0.0d0 ataj(ix,jx)=0.0d0 do 625 i=1,3 atai(ix,jx)=atai(ix,jx)+opol(IOP,ix,i)*tai(i,jx) 625 ataj(ix,jx)=ataj(ix,jx)+opol(IO ,ix,i)*taj(i,jx) do 622 ix=1,3 do 622 jx=1,3 tc(ix,jx)=0.0d0 do 622 ie=1,3 do 622 id=1,3 622 tc(ix,jx)=tc(ix,jx)- 1 EPS(jx,ie,id)*coo(IO,ie)*(atai(ix,id)-ataj(id,ix))/6.0d0 call tadd(biggpol,tc) endif c c contribution of Alpha..G to G: call ctd(IOP,opol,tem,n) call ct(taj,tt,tem) do 623 ix=1,3 do 623 jx=1,3 623 tc(jx,ix)=-(gopol(IO,jx,1)*taj(1,ix)+ 1 gopol(IO,jx,2)*taj(2,ix)+gopol(IO,jx,3)*taj(3,ix)) call tadd(biggpol,tc) c contribution of alpha..alpha to A if(lpola)then call ctd(IOP,opol,tem,n) call ct(taj,tt,tem) call ctd(IO,opol,tem,n) tv(1)=coo(IO,1)*tem(1,1)+coo(IO,2)*tem(2,1)+coo(IO,3)*tem(3,1) tv(2)=coo(IO,1)*tem(1,2)+coo(IO,2)*tem(2,2)+coo(IO,3)*tem(3,2) tv(3)=coo(IO,1)*tem(1,3)+coo(IO,2)*tem(2,3)+coo(IO,3)*tem(3,3) do 626 ix=1,3 do 626 jx=1,3 do 626 kx=1,3 su=1.5d0*( 1 (coo(IO,ix)*tem(jx,1)+coo(IO,jx)*tem(ix,1))*taj(1,kx)+ 1 (coo(IO,ix)*tem(jx,2)+coo(IO,jx)*tem(ix,2))*taj(2,kx)+ 1 (coo(IO,ix)*tem(jx,3)+coo(IO,jx)*tem(ix,3))*taj(3,kx)) if(ix.eq.jx) 1 su=su-tv(1)*taj(1,kx)-tv(2)*taj(2,kx)-tv(3)*taj(3,kx) 626 bigapol(ix,jx,kx)=bigapol(ix,jx,kx)+su endif c contribution of alpha...A to A c A_iAab Tijae a_jeg if(lpola)then call ctd(IOP,opol,tem,n) call ct(taj,tt,tem) do 627 ix=1,3 do 627 jx=1,3 do 627 kx=1,3 su=aopol(IO,ix,jx,1)*taj(1,kx) 1 + aopol(IO,ix,jx,2)*taj(2,kx)+aopol(IO,ix,jx,3)*taj(3,kx) 627 bigapol(ix,jx,kx)=bigapol(ix,jx,kx)+su endif endif endif 608 continue return end subroutine avepol(NO,NAT,opol,gopol,aopol, 1alave,coo,perx,pery,perz,lper) implicit none integer*4 IX,JX,KX,NO,NAT,IO,JO,ii real*8 opol(NAT,3,3),gopol(NAT,3,3),aopol(NAT,3,3,3), 1pi(3,3),gi(3,3),ai(3,3,3),alave(*), 1coo(NAT,3),dk,perx,pery, 1perz,pj(3,3),gj(3,3),aj(3,3,3),ddp logical lper c Polarizabilities: DO 1 IO=1,NO DO 621 IX=1,3 DO 621 JX=1,3 pi(IX,JX)=opol(IO,IX,JX) gi(IX,JX)=gopol(IO,IX,JX) do 621 KX=1,3 621 ai(IX,JX,KX)=aopol(IO,IX,JX,KX) c last index-dipole c Single-molecule averages c 1 1 .. 3 a_AA x a_AA c 2 4 .. 6 a_AA x G_AA c 3 7 .. 15 a_AB c 4 16 .. 24 a_AB x a_AB c 5 25 .. 33 a_AA x a_BB c 6 34 .. 42 G_AB c 7 43 .. 51 a_AB x G_AB c 8 52 .. 60 a_AA x G_BB c 9 61 .. 87 a_AA x A_B,AC c 10 88 ..114 a_AA x A_A,BC c 11 115 ..141 a_AB x A_A,AC do 101 IX=1,3 c alphaAA alphaAA 1.. 3 alave(IX )=alave(IX )+pi(IX,IX)*pi(IX,IX) c alphaAA GAA 4.. 6 alave(IX+ 3)=alave(IX+ 3)+pi(IX,IX)*gi(IX,IX) do 101 JX=1,3 ii=IX+3*(JX-1) c alphaAB 7.. 15 alave(ii+ 6)=alave(ii+ 6)+pi(IX,JX) c alphaAB alphaAB 16.. 24 alave(ii+ 15)=alave(ii+ 15)+pi(IX,JX)*pi(IX,JX) c alphaAA alphaBB 25.. 33 alave(ii+ 24)=alave(ii+ 24)+pi(IX,IX)*pi(JX,JX) c GAB 34.. 42 alave(ii+ 33)=alave(ii+ 33)+gi(IX,JX) c alphaAB GAB 43.. 51 alave(ii+ 42)=alave(ii+ 42)+pi(IX,JX)*gi(IX,JX) c alphaAA GBB 52.. 60 alave(ii+ 51)=alave(ii+ 51)+pi(IX,IX)*gi(JX,JX) do 101 KX=1,3 ii=IX+3*(JX-1)+9*(KX-1) c alphaAA A_B,AC 61.. 87 alave(ii+ 60)=alave(ii+ 60)+pi(IX,IX)*ai(IX,KX,JX) c alphaAA A_A,BC 88..114 alave(ii+ 87)=alave(ii+ 87)+pi(IX,IX)*ai(JX,KX,IX) c alphaAB A_A,AC 115..141 101 alave(ii+114)=alave(ii+114)+pi(IX,JX)*ai(IX,KX,IX) c A - last index-dipole c Two-molecule averages c 12 142 ..150 ai_AB x aj_AB c 13 151 ..159 2ai_AB x Gj_AB + ai_AA x (Gj_AA - Gj_BB) c 14 160 ..186 2 R_ijC ai_AA x aj_AB c 15 187 ..213 (ai_AA Aj_B,AC + ai_AA Aj_A,BC - 2 aj_AB Ai_A,AC )/3 DO 1 JO=1,NO if(JO.NE.IO)then DO 625 IX=1,3 DO 625 JX=1,3 pj(IX,JX)=opol(JO,IX,JX) gj(IX,JX)=gopol(JO,IX,JX) do 625 KX=1,3 625 aj(IX,JX,KX)=aopol(JO,IX,JX,KX) DO 626 IX=1,3 DO 626 JX=1,3 ii=IX+3*(JX-1) alave(ii+141)=alave(ii+141) +pi(IX,JX)* pj(IX,JX) alave(ii+150)=alave(ii+150)+2.0d0*pi(IX,JX)* gj(IX,JX) 1 + pi(IX,IX)*(gj(IX,IX)-gj(JX,JX)) DO 626 KX=1,3 ii=IX+3*(JX-1)+9*(KX-1) dk=2.0d0*ddp(coo,IO,JO,NAT,KX,lper,perx,pery,perz) alave(ii+159)=alave(ii+159)+dk*pi(IX,IX)*pj(IX,JX) 626 alave(ii+186)=alave(ii+186)+(pi(IX,IX)*aj(IX,KX,JX) 1 +pi(IX,IX)*aj(JX,KX,IX)-2.0d0*pj(IX,JX)*ai(IX,KX,IX))/3.0d0 c A - last index-dipole endif 1 continue return end function ddp(coo,IO,JO,n,K,lper,perx,pery,perz) implicit none integer*4 IO,JO,n,K logical lper real*8 coo(n,3),perx,pery,perz,d,ddp,p d=coo(IO,K)-coo(JO,K) if(lper)then if(K.eq.1)p=perx if(K.eq.2)p=pery if(K.eq.3)p=perz if(dabs(d+p).lt.dabs(d))d=d+p if(dabs(d-p).lt.dabs(d))d=d-p endif ddp=d return end subroutine rpd(IO,nm,A,lq,gq,aq,lqo,gqo,aqo) implicit none integer*4 nm,i,ix,jx,kx,IO,ii real*8 tl(3,3),tlo(1,3,3),A(3,3),tg(3,3),tgo(1,3,3), 1ta(3,3,3),tao(1,3,3,3),lq(*),gq(*),aq(*),lqo(*),gqo(*),aqo(*) c tlo,tgo,tao .. first fake dimension, to utilize rpp,rppp do 1 i=1,nm do 2 ix=1,3 do 2 jx=1,3 tl(ix,jx)=lq(ix+3*(jx-1)+9*(i-1)) tg(ix,jx)=gq(ix+3*(jx-1)+9*(i-1)) do 2 kx=1,3 2 ta(ix,jx,kx)=aq(ix+3*(jx-1)+9*(kx-1)+27*(i-1)) call rpp( 1,1,tlo,tl,A) call rpp( 1,1,tgo,tg,A) call rppp(1,1,tao,ta,A) do 1 ix=1,3 do 1 jx=1,3 ii=ix+3*(jx-1)+9*(i-1)+9*nm*(IO-1) lqo(ii)=tlo(1,ix,jx) gqo(ii)=tgo(1,ix,jx) do 1 kx=1,3 ii=ix+3*(jx-1)+9*(kx-1)+27*(i-1)+27*nm*(IO-1) 1 aqo(ii)=tao(1,ix,jx,kx) return end subroutine dressp(mm,nmol,rm,lp29,natm,al,npair,n,ALP29,ALPHA, 1ALPHM,opol,perx,pery,perz,px2,py2,pz2,istep,lp27,lp28,luse,natx) implicit none integer*4 mm,nmol,natm(*),npair,n,al(npair,n),L,ia,IIND,I,J, 1istep,ibx,ibxp,iby,ibyp,ibz,ibzp,ix,nib,mp,natx real*8 qx,qy,qz,rm(*),ALP29(9,3*n),ALPHA(9,nmol,3*natx), 1ALPHM(9,3*n),p1(9),p2(9),opol(n,3,3),perx,pery,perz,pi,px2,py2, 1DM,dx,dxp,dy,dyp,dz,dzp,fac,fac1,x1,x2,y1,y2,z1,z2,pz2, 1rdx,rdy,rdz,rd2,rd5,s1,s2,s3,smd,spd logical lp29,lp28,lp27,luse(*) pi=4.0d0*atan(1.0d0) c calculate cube center qx=0.0d0 qy=0.0d0 qz=0.0d0 do 1003 mm=1,nmol qx=qx+rm(1+3*(mm-1)) qy=qy+rm(2+3*(mm-1)) 1003 qz=qz+rm(3+3*(mm-1)) qx=qx/dble(nmol) qy=qy/dble(nmol) qz=qz/dble(nmol) c **************************************** IF(lp29)THEN c true superbox c c molecules in the box1: do 11517 mm=1,nmol DM=dble(natm(mm)) c save original polarizability for each molecule do 406 ia=1,natm(mm) L=al(mm,ia) DO 406 IX=1,3 IIND=3*(L-1)+IX DO 406 I=1,3 DO 406 J=1,3 406 ALP29(I+(J-1)*3,IIND)=ALPHA(I+(J-1)*3,mm,IIND) c temporary polarizability where all is collected do 410 ia=1,natm(mm) L=al(mm,ia) DO 410 IX=1,3 IIND=3*(L-1)+IX DO 410 I=1,3 DO 410 J=1,3 410 ALPHM(I+(J-1)*3,IIND)=0.0d0 c save static polarizability to faster array: do 409 I=1,3 do 409 J=1,3 409 p1(I+3*(J-1))=opol(mm,I,J) c box1 do 11516 ibx=1,3 dx=dble(ibx-2)*perx do 11516 iby=1,3 dy=dble(iby-2)*pery do 11516 ibz=1,3 dz=dble(ibz-2)*perz x1=rm(1+3*(mm-1))+dx y1=rm(2+3*(mm-1))+dy z1=rm(3+3*(mm-1))+dz c damp factor towards box walls: fac1=1.0d0 fac1=fac1*(1.0d0+cos(pi*(qx-x1)/px2/3.0d0))/2.0d0 fac1=fac1*(1.0d0+cos(pi*(qy-y1)/py2/3.0d0))/2.0d0 fac1=fac1*(1.0d0+cos(pi*(qz-z1)/pz2/3.0d0))/2.0d0 c damp original polarizability, add all to molecule mm: do 407 ia=1,natm(mm) L=al(mm,ia) DO 407 IX=1,3 IIND=3*(L-1)+IX DO 407 I=1,3 DO 407 J=1,3 407 ALPHM(I+(J-1)*3,IIND)=ALPHM(I+(J-1)*3,IIND) 1 +ALP29(I+(J-1)*3,IIND)*fac1 c molecules in the box2: do 11516 mp=1,nmol c save polarizability to faster array: do 408 I=1,3 do 408 J=1,3 408 p2(I+3*(J-1))=opol(mp,I,J) c box2 do 11516 ibxp=1,3 dxp=dble(ibxp-2)*perx do 11516 ibyp=1,3 dyp=dble(ibyp-2)*pery do 11516 ibzp=1,3 dzp=dble(ibzp-2)*perz x2=rm(1+3*(mp-1))+dxp y2=rm(2+3*(mp-1))+dyp z2=rm(3+3*(mp-1))+dzp rdx=x1-x2 rdy=y1-y2 rdz=z1-z2 rd2=rdx**2+rdy**2+rdz**2 rd5=rd2**2*dsqrt(rd2) c damp pair contributions towards box walls: fac=fac1 fac=fac*(1.0d0+cos(pi*(qx-x2)/px2/3.0d0))/2.0d0 fac=fac*(1.0d0+cos(pi*(qy-y2)/py2/3.0d0))/2.0d0 fac=fac*(1.0d0+cos(pi*(qz-z2)/pz2/3.0d0))/2.0d0 if(rd2.gt.1.0d-1)then c add all interactions to molecule mm: do 15162 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 402 IX=1,3 IIND=3*(L-1)+IX DO 402 I=1,3 DO 402 J=1,3 s1=p2(J )*ALP29(1+(I-1)*3,IIND) 1 + p2(J+3)*ALP29(2+(I-1)*3,IIND) 2 + p2(J+6)*ALP29(3+(I-1)*3,IIND) s2=p2(J)*rdx+p2(J+3)*rdy+p2(J+6)*rdz s3=ALP29(1+(I-1)*3,IIND)*rdx 1 + ALP29(2+(I-1)*3,IIND)*rdy 2 + ALP29(3+(I-1)*3,IIND)*rdz 402 ALPHM(I+(J-1)*3,IIND)=ALPHM(I+(J-1)*3,IIND)+ 1 2.0d0*fac*(3.0d0*s2*s3-s1*rd2)/rd5 c c corrections on molecular positions: DO 404 I=1,3 DO 404 J=1,3 s1=p2(I)*p1(J)+p2(I+3)*p1(J+3)+p2(I+6)*p1(J+6) s2=p2(I)*rdx+p2(I+3)*rdy+p2(I+6)*rdz s3=p1(J)*rdx+p1(J+3)*rdy+p1(J+6)*rdz smd=s3*p2(I) spd=s2*p1(J) ALPHM(I+(J-1)*3,3*(L-1)+1)=ALPHM(I+(J-1)*3,3*(L-1)+1) 1 +6.0d0*fac*(smd+spd+s1*rdx-5.0d0*s2*s3*rdx/rd2) 1 /DM/rd5 smd=s3*p2(I+3) spd=s2*p1(J+3) ALPHM(I+(J-1)*3,3*(L-1)+2)=ALPHM(I+(J-1)*3,3*(L-1)+2) 1 +6.0d0*fac*(smd+spd+s1*rdy-5.0d0*s2*s3*rdy/rd2) 1 /DM/rd5 smd=s3*p2(I+6) spd=s2*p1(J+6) 404 ALPHM(I+(J-1)*3,3*(L-1)+3)=ALPHM(I+(J-1)*3,3*(L-1)+3) 1 +6.0d0*fac*(smd+spd+s1*rdz-5.0d0*s2*s3*rdz/rd2) 1 /DM/rd5 endif 15162 continue endif 11516 continue c record resultant derivatives: do 11517 ia=1,natm(mm) L=al(mm,ia) DO 11517 IX=1,3 IIND=3*(L-1)+IX DO 11517 I=1,3 DO 11517 J=1,3 11517 ALPHA(I+(J-1)*3,mm,IIND)=ALPHM(I+(J-1)*3,IIND) c **************************************** ELSE c set superbox if required if(lp27)then nib=3 else nib=1 endif c c molecules in the box: do 1517 mm=1,nmol DM=dble(natm(mm)) c save original polarizability for each molecule do 411 ia=1,natm(mm) L=al(mm,ia) DO 411 IX=1,3 IIND=3*(L-1)+IX DO 411 I=1,3 DO 411 J=1,3 411 ALP29(I+(J-1)*3,IIND)=ALPHA(I+(J-1)*3,mm,IIND) c temporary polarizability where all is collected do 412 ia=1,natm(mm) L=al(mm,ia) DO 412 IX=1,3 IIND=3*(L-1)+IX DO 412 I=1,3 DO 412 J=1,3 412 ALPHM(I+(J-1)*3,IIND)=0.0d0 c save static polarizability to faster array: do 413 I=1,3 do 413 J=1,3 413 p1(I+3*(J-1))=opol(mm,I,J) x1=rm(1+3*(mm-1)) y1=rm(2+3*(mm-1)) z1=rm(3+3*(mm-1)) c molecules in the box and the boxes around (if lp27): do 1516 ibx=1,nib dx=dble(ibx-nib+1)*perx do 1516 iby=1,nib dy=dble(iby-nib+1)*pery do 1516 ibz=1,nib dz=dble(ibz-nib+1)*perz do 1516 mp=1,nmol if(mm.ne.mp.or. 1 (lp27.and.(ibx.ne.2.or.iby.ne.2.or.ibz.ne.2)))then c save polarizability to faster array: do 414 I=1,3 do 414 J=1,3 414 p2(I+3*(J-1))=opol(mp,I,J) x2=rm(1+3*(mp-1)) y2=rm(2+3*(mp-1)) z2=rm(3+3*(mp-1)) if(nib.eq.3)then x2=x2+dx y2=y2+dy z2=z2+dz endif rdx=x1-x2 rdy=y1-y2 rdz=z1-z2 rd2=rdx**2+rdy**2+rdz**2 rd5=rd2**2*dsqrt(rd2) fac=1.0d0 c damp towards box walls: if(lp28)then fac=fac*(1.0d0+cos(pi*(qx-x2)/px2))/2.0d0 fac=fac*(1.0d0+cos(pi*(qy-y2)/py2))/2.0d0 fac=fac*(1.0d0+cos(pi*(qz-z2)/pz2))/2.0d0 endif if(dabs(rd5).lt.1.0d-10)then write(6,*)mm,mp call report(' zero distance between molecules!') endif do 15161 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 400 IX=1,3 IIND=3*(L-1)+IX DO 400 I=1,3 DO 400 J=1,3 s1=p2(J )*ALP29(1+(I-1)*3,IIND) 1 + p2(J+3)*ALP29(2+(I-1)*3,IIND) 2 + p2(J+6)*ALP29(3+(I-1)*3,IIND) s2=p2(J)*rdx+p2(J+3)*rdy+p2(J+6)*rdz s3=ALP29(1+(I-1)*3,IIND)*rdx 1 + ALP29(2+(I-1)*3,IIND)*rdy 2 + ALP29(3+(I-1)*3,IIND)*rdz 400 ALPHM(I+(J-1)*3,IIND)=ALPHM(I+(J-1)*3,IIND)+ 1 2.0d0*fac*(3.0d0*s2*s3-s1*rd2)/rd5 c c corrections on molecular positions: DO 403 I=1,3 DO 403 J=1,3 s1=p2(I)*p1(J)+p2(I+3)*p1(J+3)+p2(I+6)*p1(J+6) s2=p2(I)*rdx+p2(I+3)*rdy+p2(I+6)*rdz s3=p1(J)*rdx+p1(J+3)*rdy+p1(J+6)*rdz ALPHM(I+(J-1)*3,3*(L-1)+1)=ALPHM(I+(J-1)*3,3*(L-1)+1) 1 +6.0d0*fac*(s3*p2(I )+s2*p1(J )+(s1-5.0d0*s2*s3/rd2)*rdx) 1 /DM/rd5 ALPHM(I+(J-1)*3,3*(L-1)+2)=ALPHM(I+(J-1)*3,3*(L-1)+2) 1 +6.0d0*fac*(s3*p2(I+3)+s2*p1(J+3)+(s1-5.0d0*s2*s3/rd2)*rdy) 1 /DM/rd5 403 ALPHM(I+(J-1)*3,3*(L-1)+3)=ALPHM(I+(J-1)*3,3*(L-1)+3) 1 +6.0d0*fac*(s3*p2(I+6)+s2*p1(J+6)+(s1-5.0d0*s2*s3/rd2)*rdz) 1 /DM/rd5 endif 15161 continue endif 1516 continue c colect back to slow array: do 1517 ia=1,natm(mm) L=al(mm,ia) DO 1517 IX=1,3 IIND=3*(L-1)+IX DO 1517 I=1,3 DO 1517 J=1,3 1517 ALPHA(I+(J-1)*3,mm,IIND)=ALP29(I+(J-1)*3,IIND)+ 1 ALPHM(I+(J-1)*3,IIND) ENDIF c **************************************** if(istep.eq.1)write(6,*)'Alpha dressed' return end subroutine fillpairs(nmol,natm,natmp,npair,n,al,alp) implicit none integer*4 nmol,natm(*),ip,m1,m2,natp,ia,npair,n,alp(npair,n), 1al(npair,n),natmp(*) ip=0 do 181 m1=1,nmol do 181 m2=m1+1,nmol ip=ip+1 natp=0 do 182 ia=1,natm(m1) natp=natp+1 182 alp(ip,natp)=al(m1,ia) do 183 ia=1,natm(m2) natp=natp+1 183 alp(ip,natp)=al(m2,ia) 181 natmp(ip)=natp do 184 ip=1,npair natm(ip)=natmp(ip) do 184 ia=1,natm(ip) 184 al(ip,ia)=alp(ip,ia) return end c do dipole derivatives from internal tensors subroutine dodd(istep,nmol,npair,n,natm,al,x,y,z,lnnames, 1iii,K0,IIM,imk,u,AM,APTM,AI,APTI,A,P) implicit none integer*4 istep,mm,nmol,ia,natm(*),npair,n,al(npair,n),I,J, 1iii(4,n),K0,IIM(K0,4,n),imk(*),L real*8 x(*),y(*),z(*),u(3,3,n),ut(3,3),AM(n,9,K0),APTM(n,9,K0), 1AI(n,3,3),APTI(n,3,3),at1(3,3),pt1(3,3),FBOHR,bohr, 1A(n,3,3),P(n,3,3) logical lnnames real*8,allocatable::xx(:),yy(:),zz(:) bohr=0.529177d0 FBOHR=0.25d0/bohr allocate(xx(n),yy(n),zz(n)) if(istep.eq.1)write(6,*)'DD from internal tensors' do 154 mm=1,nmol do 1561 ia=1,natm(mm) xx(ia)=x(al(mm,ia)) yy(ia)=y(al(mm,ia)) 1561 zz(ia)=z(al(mm,ia)) if(lnnames)then do 158 ia=1,natm(mm) do 158 I=1,4 158 iii(I,ia)=IIM(imk(mm),I,ia) endif c define the transformation matrices call dtm(iii,natm(mm),u,n,xx,yy,zz) c do 154 ia=1,natm(mm) L=al(mm,ia) c c rewrite transformation to faster arrays: DO 14121 I=1,3 DO 14121 J=1,3 14121 ut(I,J)=u(I,J,ia) if(lnnames)then DO 141 I=1,3 DO 141 J=1,3 at1(I,J)= AM( ia,I+(J-1)*3,imk(mm)) 141 pt1(I,J)=APTM(ia,I+(J-1)*3,imk(mm)) else DO 1413 I=1,3 DO 1413 J=1,3 at1(I,J)= AI(ia,I,J) 1413 pt1(I,J)=APTI(ia,I,J) endif c call tas2(at1,ut) call tas2(pt1,ut) c c Transform to common origin: DO 2201 J=1,3 at1(J,1)=at1(J,1)+(y(ia)*pt1(J,3)-z(ia)*pt1(J,2))*FBOHR at1(J,2)=at1(J,2)+(z(ia)*pt1(J,1)-x(ia)*pt1(J,3))*FBOHR 2201 at1(J,3)=at1(J,3)+(x(ia)*pt1(J,2)-y(ia)*pt1(J,1))*FBOHR c c record DO 154 I=1,3 DO 154 J=1,3 A(L,I,J)=at1(I,J) 154 P(L,I,J)=pt1(I,J) return end c dipole derivatives from charges: subroutine doddq(istep,n,pchgft,P,A,x,y,z) implicit none integer*4 istep,n,ia,ix,jx real*8 A(n,3,3),P(n,3,3),x(*),y(*),z(*),fac,pchgft(*) if(istep.eq.1)write(6,*)'DD from charges' do 401 ia=1,n fac=pchgft(ia) do 405 ix=1,3 do 405 jx=1,3 P(ia,ix,jx)=0.0d0 405 A(ia,ix,jx)=0.0d0 P(ia,1,1)=fac*x(ia) P(ia,2,2)=fac*y(ia) P(ia,3,3)=fac*z(ia) A(ia,1,2)= fac*z(ia) A(ia,2,1)=-fac*z(ia) A(ia,2,3)= fac*x(ia) A(ia,3,2)=-fac*x(ia) A(ia,3,1)= fac*y(ia) 401 A(ia,1,3)=-fac*y(ia) return end c Alpha from atomic polarizabilities: subroutine doaap(istep,npair,n,natm,al,luse,AA1,G1, 1ALPHA1,x,y,z,Xwork,nmol,AA,G,ALPHA,pol,natx) implicit none integer*4 istep,npair,natm(*),n,al(npair,n),mm,ia,L,IX, 1IIND,I,J,iap,LP,nmol,K,natx,i2 logical luse(*) real*8 AA1(3*n,3,3,3),G1(3*n,3,3),ALPHA1(3*n,3,3),rki(3), 1x(*),y(*),z(*),rd1,rd67,s1,Xwork(3,n),AA(27,nmol,3*natx), 2G(9,nmol,3*natx),ALPHA(9,nmol,3*natx),pol(*) if(istep.eq.1)write(6,*)'Alpha from atomic polarizabilities' c c molecular polarizabilities: DO 17 mm=1,npair DO 179 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 172 IX=1,3 IIND=3*(L-1)+IX DO 172 I=1,3 DO 172 J=1,3 AA1(IIND,I,J,1)=0.0d0 AA1(IIND,I,J,2)=0.0d0 AA1(IIND,I,J,3)=0.0d0 G1(IIND,I,J)=0.0d0 172 ALPHA1(IIND,I,J)=0.0d0 endif 179 continue DO 192 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 178 iap=1,natm(mm) LP=al(mm,iap) if(luse(LP))then rki(1)=x(LP)-x(L) rki(2)=y(LP)-y(L) rki(3)=z(LP)-z(L) if(LP.ne.L)then rd1=dsqrt(rki(1)**2+rki(2)**2+rki(3)**2+1.0d-9) rd67=6.0d0/rd1**7*pol(L)*pol(LP) DO 173 IX=1,3 IIND=3*(L-1)+IX DO 173 I=1,3 DO 173 J=1,3 s1=-5.0d0*rki(I)*rki(J)*rki(IX) if(I.eq.IX)s1=s1+rd1*rd1*rki(J) if(J.eq.IX)s1=s1+rd1*rd1*rki(I) if(I.eq.J )s1=s1+rd1*rd1*rki(IX) 173 ALPHA1(IIND,I,J)=ALPHA1(IIND,I,J)+s1*rd67 endif endif 178 continue endif 192 continue CALL TTTT(ALPHA1,AA1,G1,n,2,x,y,z,Xwork) c c transcript into the big arrays: DO 193 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 194 IX=1,3 IIND=3*(L-1)+IX i2=3*(ia-1)+IX DO 194 I=1,3 DO 194 J=1,3 DO 195 K=1,3 195 AA(I+(J-1)*3+(K-1)*9,mm,i2)=AA1( IIND,I,J,K) G(I+(J-1)*3, mm,i2)=G1( IIND,I,J) 194 ALPHA(I+(J-1)*3, mm,i2)=ALPHA1(IIND,I,J) endif 193 continue 17 continue return end c Polarizabilities and derivatives from internal tensors: subroutine doai(istep,n,nmol,natm,al,x,y,z,lnnames,K0,IIM,imk, 1POLM,POLGM,POLAM,POLI,POLGI,POLAI,opol,gopol,aopol,ALPHAM,GM,AAM, 1ALPHA1,AA1,G1,Xwork,luse,ALPHA,AA,G,LNAMES,ALPHAI,GI,AAI,npair, 1ALP29,ALPHM,perx,pery,perz,px2,py2,pz2,lp27,lp28,lp29,ldress, 1lq,gq,aq,lqo,gqo,aqo,coo,nm,natx) implicit none integer*4 istep,nmol,natm(*),n,npair,al(npair,n),mm,ia,L,i2, 1K0,IIM(K0,4,n),imk(*),IXP,IP,JP,KP,I,IIND,J,K,ix,nm,im,natx real*8 DM,x(*),y(*),z(*),bohr,ut(3,3),polt(3,3),POLM(K0,9), 1POLGM(K0,9),POLAM(K0,27),POLI(3,3),POLGI(3,3),POLAI(3,3,3), 1pola(3,3,3),opol(n,3,3),gopol(n,3,3),aopol(n,3,3,3),poltg(3,3), 1as1(3,3,3),gs1(3,3,3),au1(3,3,3,3),Xwork(3,n),AAI(3*n,3,3,3), 4ALPHAM(9,K0,3*n),GM(9,K0,3*n),AAM(27,K0, 3*n),GI(3*n,3,3), 2ALPHA1(3*n,3,3),AA1(3*n,3,3,3),G1(3*n,3,3),ALPHAI(3*n,3,3), 3ALPHA(9,nmol,3*natx),AA(27,nmol,3*natx),G(9,nmol,3*natx), 1ALP29(9,3*n), ALPHM(9,3*n),perx,pery,perz,px2,py2,pz2, 1lq(*),gq(*),aq(*),lqo(*),gqo(*),aqo(*),coo(n,3) logical lnames,luse(*),lp29,lp27,lp28,lnnames,ldress real*8,allocatable::rm(:),xx(:),yy(:),zz(:),u(:,:,:) integer*4,allocatable::iii(:,:) allocate(rm(3*n),xx(n),yy(n),zz(n),iii(4,n),u(3,3,n)) bohr=0.529177d0 if(istep.eq.1)write(6,*)'Alpha from internal tensors' c molecular positions in atomic units (x y z are in A): rm=0.0d0 DO 196 mm=1,nmol DM=dble(natm(mm)) DO 198 ia=1,natm(mm) L=al(mm,ia) rm(1+(mm-1)*3)=rm(1+(mm-1)*3)+x(L) rm(2+(mm-1)*3)=rm(2+(mm-1)*3)+y(L) 198 rm(3+(mm-1)*3)=rm(3+(mm-1)*3)+z(L) do 196 ix=1,3 rm(ix+(mm-1)*3)=rm(ix+(mm-1)*3)/DM/bohr 196 coo(mm,ix)=rm(ix+(mm-1)*3) c propagate the internal tensors to all molecules, c order of atoms must be the same do 1515 mm=1,nmol c define the transformation matrices do 156 ia=1,natm(mm) xx(ia)=x(al(mm,ia)) yy(ia)=y(al(mm,ia)) 156 zz(ia)=z(al(mm,ia)) if(lnnames)then do 159 ia=1,natm(mm) do 159 I=1,4 159 iii(I,ia)=IIM(imk(mm),I,ia) endif c if not enough bonds: if(natm(mm).lt.2)then u=0.0d0 DO 1551 ia=1,natm(mm) DO 1551 I=1,3 1551 u(I,I,ia)=1.0d0 else call dtm(iii,natm(mm),u,n,xx,yy,zz) endif do 15 ia=1,natm(mm) L=al(mm,ia) c c transformation matrix for this atom: DO 1412 I=1,3 DO 1412 J=1,3 1412 ut(I,J)=u(I,J,ia) c c Polarizabilities: c polarizability according to the system of first atom: if(ia.eq.1)then if(lnames)then DO 20 I=1,3 DO 20 J=1,3 poltg(I,J)=POLGM(imk(mm),I+(J-1)*3) polt(I,J)=POLM(imk(mm),I+(J-1)*3) DO 20 K=1,3 20 pola(I,J,K)=POLAM(imk(mm),I+(J-1)*3+9*(K-1)) else polt=POLI poltg=POLGI pola=POLAI endif call tas2(polt,ut) call tas2(poltg,ut) call tas3(pola,ut) DO 21 I=1,3 DO 21 J=1,3 opol( mm,I,J)=polt( I,J) gopol(mm,I,J)=poltg(I,J) DO 21 K=1,3 21 aopol(mm,I,J,K)=pola(I,J,K) c normal mode polarizability derivatives: DO 2 im=1,nm DO 1 i=1,3 DO 1 j=1,3 polt( i,j )=lq(i+3*(j-1)+9*(im-1)) poltg(i,j )=gq(i+3*(j-1)+9*(im-1)) DO 1 k=1,3 1 pola(i,j,k)=aq(i+3*(j-1)+9*(k-1)+27*(im-1)) call tas2(polt,ut) call tas2(poltg,ut) call tas3(pola,ut) DO 2 i=1,3 DO 2 j=1,3 lqo(i+3*(j-1)+9*(im-1)+9*nm*(mm-1))=polt( i,j) gqo(i+3*(j-1)+9*(im-1)+9*nm*(mm-1))=poltg(i,j) DO 2 k=1,3 2 aqo(i+3*(j-1)+9*(k-1)+27*(im-1)+27*nm*(mm-1))=pola(i,j,k) endif c c Derivatives: c rewrite ALPHA G A tensors to temporary arrays: if(lnnames)then DO 151 IXP=1,3 IIND=3*(ia-1)+IXP DO 151 IP=1,3 DO 151 JP=1,3 as1(IXP,IP,JP)=ALPHAM(IP+(JP-1)*3,imk(mm),IIND) gs1(IXP,IP,JP)=GM(IP+(JP-1)*3,imk(mm),IIND) DO 151 KP=1,3 151 au1(IXP,IP,JP,KP)=AAM(IP+(JP-1)*3+(KP-1)*3,imk(mm),IIND) else DO 1512 IXP=1,3 IIND=3*(ia-1)+IXP DO 1512 IP=1,3 DO 1512 JP=1,3 as1(IXP,IP,JP )=ALPHAI(IIND,IP,JP) gs1(IXP,IP,JP )=GI( IIND,IP,JP) DO 1512 KP=1,3 1512 au1(IXP,IP,JP,KP)=AAI( IIND,IP,JP,KP) endif c c transform as1,gs1,au1 call tas3(as1,ut) call tas3(gs1,ut) call tas4(au1,ut) c rewrite DO 153 IX=1,3 DO 153 I=1,3 DO 153 J=1,3 ALPHA1(3*(L-1)+IX,I,J )=as1(IX,I,J ) G1( 3*(L-1)+IX,I,J )=gs1(IX,I,J ) DO 153 K =1,3 153 AA1( 3*(L-1)+IX,I,J,K)=au1(IX,I,J,K) 15 continue c ia c transform to common origin: CALL TTTT(ALPHA1,AA1,G1,n,2,x,y,z,Xwork) do 1515 ia=1,natm(mm) L=al(mm,ia) if(luse(L))then DO 300 IX=1,3 IIND=3*(L-1)+IX i2=3*(ia-1)+IX DO 300 I=1,3 DO 300 J=1,3 DO 306 K=1,3 306 AA(I+(J-1)*3+(K-1)*9,mm,i2)=AA1( IIND,I,J,K) G (I+(J-1)*3, mm,i2)=G1 (IIND,I,J) 300 ALPHA(I+(J-1)*3, mm,i2)=ALPHA1(IIND,I,J) endif 1515 continue c mm c c dress the polarizabilities to account for the influence of the c other molecules if(ldress) 1call dressp(mm,nmol,rm,lp29,natm,al,npair,n,ALP29,ALPHA, 1ALPHM,opol,perx,pery,perz,px2,py2,pz2,istep,lp27,lp28,luse,natx) if(istep.eq.1)write(6,*)'Done' return end c TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT subroutine tas2(s,u) integer*4 A,B real*8 s(3,3),u(3,3),t(3,3) DO 21 A=1,3 DO 21 B=1,3 21 t(A,B)=u(1,A)*s(1,B)+u(2,A)*s(2,B)+u(3,A)*s(3,B) DO 22 A=1,3 DO 22 B=1,3 22 s(A,B)=u(1,B)*t(A,1)+u(2,B)*t(A,2)+u(3,B)*t(A,3) return end subroutine tas3(s,u) implicit none integer*4 A,B,C real*8 t(3,3,3),s(3,3,3),u(3,3) DO 157 A=1,3 DO 157 B=1,3 DO 157 C=1,3 157 t(A,B,C)=u(1,A)*s(1,B,C)+u(2,A)*s(2,B,C)+u(3,A)*s(3,B,C) DO 152 A=1,3 DO 152 B=1,3 DO 152 C=1,3 152 s(A,B,C)=u(1,B)*t(A,1,C)+u(2,B)*t(A,2,C)+u(3,B)*t(A,3,C) DO 151 A=1,3 DO 151 B=1,3 DO 151 C=1,3 151 t(A,B,C)=u(1,C)*s(A,B,1)+u(2,C)*s(A,B,2)+u(3,C)*s(A,B,3) s=t return end subroutine tas4(s,u) implicit none integer*4 A,B,C,D real*8 t(3,3,3,3),s(3,3,3,3),u(3,3) DO 157 A=1,3 DO 157 B=1,3 DO 157 C=1,3 DO 157 D=1,3 157 t(A,B,C,D)=u(1,A)*s(1,B,C,D)+u(2,A)*s(2,B,C,D)+u(3,A)*s(3,B,C,D) DO 152 A=1,3 DO 152 B=1,3 DO 152 C=1,3 DO 152 D=1,3 152 s(A,B,C,D)=u(1,B)*t(A,1,C,D)+u(2,B)*t(A,2,C,D)+u(3,B)*t(A,3,C,D) DO 151 A=1,3 DO 151 B=1,3 DO 151 C=1,3 DO 151 D=1,3 151 t(A,B,C,D)=u(1,C)*s(A,B,1,D)+u(2,C)*s(A,B,2,D)+u(3,C)*s(A,B,3,D) DO 153 A=1,3 DO 153 B=1,3 DO 153 C=1,3 DO 153 D=1,3 153 s(A,B,C,D)=u(1,D)*t(A,B,C,1)+u(2,D)*t(A,B,C,2)+u(3,D)*t(A,B,C,3) return end subroutine rpp(IOX,IO,o,g,A) implicit none integer*4 IX,JX,b,IOX,IO real*8 o(IOX,3,3),A(3,3),g(3,3),t(3,3) do 1 IX=1,3 do 1 b=1,3 1 t(IX,b)=A(IX,1)*g(1,b)+A(IX,2)*g(2,b)+A(IX,3)*g(3,b) do 2 IX=1,3 do 2 JX=1,3 2 o(IO,IX,JX)=A(JX,1)*t(IX,1)+A(JX,2)*t(IX,2)+A(JX,3)*t(IX,3) return end subroutine rppp(n,L,o,p,A) implicit none integer*4 i,j,b,n,L,c,k real*8 o(n,3,3,3),A(3,3),p(3,3,3),t(3,3,3),u(3,3,3) do 1 i=1,3 do 1 b=1,3 do 1 c=1,3 1 t( i,b,c)=A(i,1)*p(1,b,c)+A(i,2)*p(2,b,c)+A(i,3)*p(3,b,c) do 2 i=1,3 do 2 j=1,3 do 2 c=1,3 2 u( i,j,c)=A(j,1)*t(i,1,c)+A(j,2)*t(i,2,c)+A(j,3)*t(i,3,c) do 3 i=1,3 do 3 j=1,3 do 3 k=1,3 3 o(L,i,j,k)=A(k,1)*u(i,j,1)+A(k,2)*u(i,j,2)+A(k,3)*u(i,j,3) return end c correct for periodic jumps: subroutine cpj(n,lcell,ac,bc,cc,x,y,z,xm,ym,zm,xmm,ymm,zmm, 1perx,pery,perz,px2,py2,pz2) implicit none logical lcell integer*4 n,ia real*8 ac(3),bc(3),cc(3),aan,abn,acn,x(*),y(*),z(*), 1ap,bp,cp,apm,bpm,cpm,xm(*),ym(*),zm(*),perx,pery,perz,px2,py2,pz2, 2xmm(*),ymm(*),zmm(*) if(lcell)then aan=dsqrt(ac(1)**2+ac(2)**2+ac(3)**2) abn=dsqrt(bc(1)**2+bc(2)**2+bc(3)**2) acn=dsqrt(cc(1)**2+cc(2)**2+cc(3)**2) do 185 ia=1,n ap=(x(ia)*ac(1)+y(ia)*ac(2)+z(ia)*ac(3))/aan bp=(x(ia)*bc(1)+y(ia)*bc(2)+z(ia)*bc(3))/abn cp=(x(ia)*cc(1)+y(ia)*cc(2)+z(ia)*cc(3))/acn apm=(xm(ia)*ac(1)+ym(ia)*ac(2)+zm(ia)*ac(3))/aan bpm=(xm(ia)*bc(1)+ym(ia)*bc(2)+zm(ia)*bc(3))/abn cpm=(xm(ia)*cc(1)+ym(ia)*cc(2)+zm(ia)*cc(3))/acn if(ap-apm.gt.aan)then call trabc(ia,xm ,ym ,zm ,ac,1.0d0) call trabc(ia,xmm,ymm,zmm,ac,1.0d0) endif if(apm-ap.gt.aan)then call trabc(ia,xm ,ym ,zm ,ac,-1.0d0) call trabc(ia,xmm,ymm,zmm,ac,-1.0d0) endif if(bp-bpm.gt.abn)then call trabc(ia,xm , ym, zm,bc,1.0d0) call trabc(ia,xmm,ymm,zmm,bc,1.0d0) endif if(bpm-bp.gt.abn)then call trabc(ia, xm, ym, zm,bc,-1.0d0) call trabc(ia,xmm,ymm,zmm,bc,-1.0d0) endif if(cp-cpm.gt.acn)then call trabc(ia, xm, ym, zm,cc,1.0d0) call trabc(ia,xmm,ymm,zmm,cc,1.0d0) endif if(cpm-cp.gt.acn)then call trabc(ia, xm, ym, zm,cc,-1.0d0) call trabc(ia,xmm,ymm,zmm,cc,-1.0d0) endif 185 continue else do 180 ia=1,n if(x(ia)-xm(ia).gt.px2)then xm(ia)=xm(ia)+perx xmm(ia)=xmm(ia)+perx endif if(xm(ia)-x(ia).gt.px2)then xm(ia)=xm(ia)-perx xmm(ia)=xmm(ia)-perx endif if(y(ia)-ym(ia).gt.py2)then ym(ia)=ym(ia)+pery ymm(ia)=ymm(ia)+pery endif if(ym(ia)-y(ia).gt.py2)then ym(ia)=ym(ia)-pery ymm(ia)=ymm(ia)-pery endif if(z(ia)-zm(ia).gt.pz2)then zm(ia)=zm(ia)+perz zmm(ia)=zmm(ia)+perz endif if(zm(ia)-z(ia).gt.pz2)then zm(ia)=zm(ia)-perz zmm(ia)=zmm(ia)-perz endif 180 continue endif return end c TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT subroutine ferr(err,rerr,a,s,ii) implicit none integer*4 ii,i real*8 err,rerr,a(*),s(*),t,m t=0.0d0 m=0.0d0 do 1 i=1,3 if(dabs(s(i)).gt.m)m=dabs(s(i)) 1 t=t+(a(ii+i)-s(i))**2 err=dsqrt(t/3.0d0) if(m.gt.0.0d0)then rerr=(err/m)*100.0d0 else rerr=err endif return end subroutine serr(err,rerr,a,s,ii) implicit none integer*4 ii,i,j,ij real*8 err,rerr,a(*),s(3,3),t,m t=0.0d0 m=0.0d0 do 1 i=1,3 do 1 j=1,3 ij=3*(j-1)+i if(dabs(s(i,j)).gt.m)m=dabs(s(i,j)) 1 t=t+(a(ii+ij)-s(i,j))**2 err=dsqrt(t/3.0d0) if(m.gt.0.0d0)then rerr=(err/m)*100.0d0 else rerr=err endif return end subroutine terr(err,rerr,a,s,ii) implicit none integer*4 ii,i,j,k,ij real*8 err,rerr,a(*),s(3,3,3),t,m t=0.0d0 m=0.0d0 do 1 i=1,3 do 1 j=1,3 do 1 k=1,3 ij=9*(k-1)+3*(j-1)+i if(dabs(s(i,j,k)).gt.m)m=dabs(s(i,j,k)) 1 t=t+(a(ii+ij)-s(i,j,k))**2 err=dsqrt(t/3.0d0) if(m.gt.0.0d0)then rerr=(err/m)*100.0d0 else rerr=err endif return end subroutine ssu(nmol,st,ic,sp,s1,s2,s3,sG,s1g,s2g,s3g) implicit none integer*4 nmol,ic,I,J real*8 st(3,3),sp,s1,s2,s3,sG,s1g,s2g,s3g,DE if(ic.eq.3)then do 203 I=1,3 do 203 J=1,3 203 st(I,J)=sp*DE(I,J) endif if(ic.eq.4)then do 204 I=1,3 do 204 J=1,3 204 st(I,J)=(s3 -s1 )*DE(I,J)+s1 endif if(ic.eq.5)then do 205 I=1,3 do 205 J=1,3 205 st(I,J)=(s3 -s2 )*DE(I,J)+s2 endif if(ic.eq.6)then do 206 I=1,3 do 206 J=1,3 206 st(I,J)=sG*DE(I,J) endif if(ic.eq.7)then do 207 I=1,3 do 207 J=1,3 207 st(I,J)=(s3g-s1g)*DE(I,J)+s1g endif if(ic.eq.8)then do 208 I=1,3 do 208 J=1,3 208 st(I,J)=(s3g-s2g)*DE(I,J)+s2g endif if(ic.eq.12)then do 212 I=1,3 do 212 J=1,3 212 st(I,J)=sp*sp*dble(nmol-1)*DE(I,J) endif if(ic.eq.13)then do 213 I=1,3 do 213 J=1,3 213 st(I,J)=2.0d0*sp*sG*dble(nmol-1)*DE(I,J) endif return end subroutine ssw(s1a,s3a,sw,ic) implicit none integer*4 ic,I,J,K real*8 s1a,s3a,sw(3,3,3),EPS sw=0.0d0 do 1 I=1,3 do 1 J=1,3 do 1 K=1,3 if(ic.eq. 9)sw(I,J,K)=EPS(I,J,K)*s1a 1 if(ic.eq.11)sw(I,J,K)=EPS(I,J,K)*s3a return end subroutine rpca(pkind,ic,alave,istep) implicit none integer*4 pkind,ic,istep real*8 alave(pkind) logical lex if(ic.eq.0)then inquire(file='ALAVE.SCR',exist=lex) if(lex)then open(9,file='ALAVE.SCR',form='unformatted') read(9)alave read(9)istep write(6,*)'ALAVE.SCR ',istep,' steps read' else istep=0 write(6,*)'ALAVE.SCR not found' endif else open(9,file='ALAVE.SCR',form='unformatted') write(9)alave write(9)istep endif close(9) return end