program rdfu implicit none integer*4 nat,ig,ia,n,nats,nmol,np,irep,ics,nrec real*8 rmin,rmax,rlim,dr,perx,pery,perz,cg(3,3),fg(3,3),cs(3,3), 1ic,ir,apo(3,3,3),po(3,3),gpo(3,3) logical lper,lext,lpol,lrest real*8,allocatable::r(:),g(:),cc(:,:,:),ci(:,:),ue(:),xe(:),ye(:), 1ur(:),xr(:),yr(:),we(:),wr(:),rrec(:),drec(:),arec(:) integer*4,allocatable::icc(:),z(:) c write(6,6000) 6000 format(/, 'Probability distribution - two molecule position',/,/, 1 ' Input: FILE.X ... MD snapshots',/, 1 ' RDFU.OPT ... options',/, 1 ' Output: RDFM.PRN ... radial distribution function',/, 1 ' XY.PRN ... directional angles',/, 1 ' X.PRN ... chirality index',/) call readopt(nats,np,rmin,rmax,rlim,lper,perx,pery,perz,irep,lext, 1lpol,lrest,nrec) allocate(rrec(2*nats*3*nrec),drec(nrec),arec(nrec)) drec=1.0d90 arec=0.0d0 dr=(rmax-rmin)/dble(np-1) ig=0 ic=0.0d0 ir=0.0d0 allocate(g(np),ci(np,3),cc(np,3,3),icc(np)) c g - RDF g=0.0d0 c cg - directional cosines, close cg=0.0d0 c cg - directional cosines, far fg=0.0d0 c cs - directional cosines, single molecule cs=0.0d0 c cc - directional cosines, distance dependence cc=0.0d0 c ci - chirality index, distance dependence ci=0.0d0 icc=0 ics=0 c for extensive cosinus product averaging,allocate even when not used: allocate(ue(81),xe(81),ye(243),ur(81),xr(81),yr(243), 1we(81),wr(81)) ue=0.0d0 xe=0.0d0 ye=0.0d0 ur=0.0d0 xr=0.0d0 yr=0.0d0 we=0.0d0 wr=0.0d0 if(lrest)then write(6,*)'restart requested, averages from AV.SCR' open(9,file='AV.SCR',form='unformatted') read(9)we,ue,xe,ye close(9) goto 999 endif if(lpol)call ldp(nats,po,apo,gpo) open(9,file='FILE.X',status='old') 1 read(9,90,end=99,err=99) 90 format(a80) read(9,*,end=99,err=99)nat ig=ig+1 n=3*nat nmol=nat/nats if(ig.eq.1)then allocate(r(n),z(nat)) write(6,602)nmol 602 format(I6,' molecules') endif do 2 ia=1,nat 2 read(9,*)z(ia),r(3*ia-2),r(3*ia-1),r(3*ia) call incr(g,z,r,nats,nmol,lper,perx,pery,perz,dr,np,rmin, 1rlim,cg,fg,ic,ir,cc,icc,cs,ics,ci,ue,xe,ye,ur,xr,yr,lext, 1we,wr,drec,nrec,rrec,arec) c write if(mod(ig,irep).eq.0)then call wcos(cg,fg,cs,ic,ir,ics,lext,ue,xe,ye, 1 we,ur,xr,yr,wr,nmol,ig) if(lpol)call wrp(we,ue,xe,ye,po,apo,gpo) endif goto 1 99 close(9) write(6,600)ig,nat 600 format(i6,' geometries,',i6,' atoms') call wrc(g,np,rmin,dr,nmol,ig,cc,icc,ci,lper,perx,pery,perz) call aar(ig,nmol,we,ue,xe,ye) call wrrec(nats,nrec,drec,rrec,arec,z) 999 if(lpol)call wrp(we,ue,xe,ye,po,apo,gpo) end subroutine wrrec(nats,nrec,drec,rrec,arec,z) implicit none integer*4 nats,nrec,ir,ia,ix,ii,z(*) real*8 drec(*),rrec(*),arec(*),pi pi=4.0d0*datan(1.0d0) open(9,file='REC.X') do 1 ir=1,nrec write(9,90)ir,drec(ir),arec(ir),acos(arec(ir))*180.0d0/pi,nats*2 90 format(i6,' dist: ',3f12.4,/,i6) do 2 ia=1,nats ii=3*(ia-1)+3*nats*(ir-1) 2 write(9,91)z(ia),(rrec(ix+ii),ix=1,3) do 1 ia=1,nats ii=3*(ia-1)+3*nats*(ir-1)+nrec*3*nats 1 write(9,91)z(ia),(rrec(ix+ii),ix=1,3) 91 format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0') close(9) return end subroutine ldp(nats,po,apo,gpo) implicit none integer*4 nats real*8 apo(3,3,3),po(3,3),gpo(3,3),fr real*8,allocatable::RS(:,:) integer*4,allocatable::q(:) allocate(RS(3,nats),q(nats)) CALL READRXS(nats,RS,q) call readpolal(po,'SMALL.POL',fr) call readpolg(gpo,'SMALL.G.POL') call readpola(apo,'SMALL.A.POL') c transform G, A to small molecule center call trs(nats,RS,po,apo,gpo,q) c transform al, G, A to molecular axes: call tr2(nats,RS,po,apo,gpo,q) return end subroutine wrp(we,ue,xe,ye,po,apo,gpo) implicit none integer*4 a,b,c,d,e,ii real*8 we(*),ue(*),xe(*),ye(*),apo(3,3,3),po(3,3),gpo(3,3),AMU, 1EXCA,BOHR,gpisvejc,aa,Gaa,ab,Gab,s1,s3, 1s1a,S0e,S31e,S32e,T0,T31,T32,T33,T3,ynorm EXCA=5320.0d0 AMU=1822.0d0 BOHR=0.529177d0 gpisvejc=(AMU*BOHR**5)*2.0d0*4.0d0*atan(1.0d0)/EXCA aa= po(1,1)+ po(2,2)+ po(3,3) Gaa=gpo(1,1)+gpo(2,2)+gpo(3,3) ab =po(1,1)* po(1,1)+po(1,2)* po(1,2)+po(1,3)* po(1,3) 1 +po(2,1)* po(2,1)+po(2,2)* po(2,2)+po(2,3)* po(2,3) 1 +po(3,1)* po(3,1)+po(3,2)* po(3,2)+po(3,3)* po(3,3) Gab=po(1,1)*Gpo(1,1)+po(1,2)*Gpo(1,2)+po(1,3)*Gpo(1,3) 1 +po(2,1)*Gpo(2,1)+po(2,2)*Gpo(2,2)+po(2,3)*Gpo(2,3) 1 +po(3,1)*Gpo(3,1)+po(3,2)*Gpo(3,2)+po(3,3)*Gpo(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 s3 =(2.0d0* aa* aa +4.0d0* ab)/30.0d0 c EPS_abc a_cd A_bda / 15: s1a=0.0d0 do 14 d=1,3 14 s1a=s1a+(po(3,d)*apo(2,d,1) -po(2,d)*apo(3,d,1) 1 +po(1,d)*apo(3,d,2) -po(3,d)*apo(1,d,2) 2 +po(2,d)*apo(1,d,3) -po(1,d)*apo(2,d,3))/15.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 write(6,*)s1,s3 S0e=s1+s3 c (7 ab ab + aa bb ) / 30 S31e=2.0d0/15.0d0*(3.0d0*Gab-aa*Gaa) c 2/15 (3 a_ab G_ab-a_aa*G_bb) S32e=2.0d0*s1a/3.0d0 S0e =S0e*180.0d0*AMU*BOHR**4 S31e=S31e*gpisvejc*360.0d0 S32e=S32e*gpisvejc*360.0d0 write(6,6003)S0e,S31e,S32e,S31e+S32e 6003 format(' RamSCP180 :',e11.3,/,/, 1 ' ROASCP180G:',e11.3,/, 1 ' A:',e11.3,/, 1 ' ',11(1h-),11(1h-),/, 1 ' T:',e11.3,/) T0=0.0d0 T31=0.0d0 T32=0.0d0 T33=0.0d0 ii=0 do 1 a=1,3 do 1 b=1,3 do 1 d=1,3 do 1 e=1,3 ii=ii+1 T0=T0+po(a,b)*po(d,e)*we(ii) T31=T31-po(a,b)*po(d,e)*xe(ii) 1 T32=T32+po(a,b)*gpo(d,e)*ue(ii) ynorm=0.0d0 ii=0 do 2 a=1,3 do 2 b=1,3 do 2 c=1,3 do 2 d=1,3 do 2 e=1,3 ii=ii+1 ynorm=ynorm+ye(ii)**2 2 T33=T33+po(a,b)*apo(d,e,c)*ye(ii) c apo-dipole index last write(6,6018)dsqrt(ynorm) 6018 format(' ynorm',e12.3) T0=T0*6.0d0*AMU*BOHR**4 T31=T31/15.0d0*gpisvejc*360.0d0 T32=T32* 0.4d0*gpisvejc*360.0d0 T33=T33/90.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,/) return end subroutine trs(nats,RS,po,apo,gpo,q) implicit none integer*4 nats,ix,ia,MENDELEV,q(*) real*8 RS(3,nats),po(3,3),apo(3,3,3),gpo(3,3),m,c(3),bohr parameter (MENDELEV=89) real*4 mass(MENDELEV) data mass/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ bohr=0.529177d0 do 313 ix=1,3 c(ix)=0.0d0 m=0.0d0 do 3131 ia=1,NATS m=m+dble(mass(q(ia))) 3131 c(ix)=c(ix)+RS(ix,ia)*dble(mass(q(ia))) 313 c(ix)=c(ix)/bohr/m call TTT2(po,apo,gpo,1,c,1) return end subroutine tr2(nats,RS,po,apo,gpo,q) implicit none integer*4 nats,q(*),i,ix,ia,MENDELEV,jx,IERR,kx real*8 RS(3,nats),po(3,3),apo(3,3,3),gpo(3,3),tm,c(3),ma, 1p(3),u(3,3),A(3,3),e(3),av(3),bv(3),cv(3),sp,au(3,3),xi, 1t(3,3),y(3,3,3) parameter (MENDELEV=89) real*4 mass(MENDELEV) data mass/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ c total mass tm, mass center c: tm=0.0d0 c=0.0d0 do 3 ia=1,nats ma=dble(mass(q(ia))) tm=tm+ma do 3 ix=1,3 3 c(ix)=c(ix)+ma*RS(ix,ia) do 2 ix=1,3 2 c(ix)=c(ix)/tm c p = Qi xi : dipole moment: c u = Qi xi yi: moment of inertia: p=0.0d0 u=0.0d0 do 5 ia=1,nats ma=dble(mass(q(ia))) do 5 ix=1,3 xi=ma*(RS(ix,ia)-c(ix)) p(ix)=p(ix)+xi do 5 jx=1,3 5 u(ix,jx)=u(ix,jx)+xi*(RS(jx,ia)-c(jx)) c principal axes CALL TRED12(3,u,A,e,2,IERR) do 6 ix=1,3 av(ix)=A(ix,1) 6 bv(ix)=A(ix,2) c check that a-axis points along p if(sp(av,p).lt.0.0d0)then av(1)=-av(1) av(2)=-av(2) av(3)=-av(3) endif c check that b-axis points along p if(sp(bv,p).lt.0.0d0)then bv(1)=-bv(1) bv(2)=-bv(2) bv(3)=-bv(3) endif c c = a x b: call vp(av,bv,cv) do 10 ix=1,3 au(ix,1)=av(ix) au(ix,2)=bv(ix) 10 au(ix,3)=cv(ix) do 4 ix=1,3 do 4 jx=1,3 t(ix,jx)=0.0d0 do 4 i=1,3 4 t(ix,jx)=t(ix,jx)+au(i,jx)*po(ix,i) do 7 ix=1,3 do 7 jx=1,3 po(ix,jx)=0.0d0 do 7 i=1,3 7 po(ix,jx)=po(ix,jx)+au(i,ix)*t(i,jx) do 8 ix=1,3 do 8 jx=1,3 t(ix,jx)=0.0d0 do 8 i=1,3 8 t(ix,jx)=t(ix,jx)+au(i,jx)*gpo(ix,i) do 9 ix=1,3 do 9 jx=1,3 gpo(ix,jx)=0.0d0 do 9 i=1,3 9 gpo(ix,jx)=gpo(ix,jx)+au(i,ix)*t(i,jx) do 81 ix=1,3 do 81 jx=1,3 do 81 kx=1,3 y(ix,jx,kx)=0.0d0 do 81 i=1,3 81 y(ix,jx,kx)=y(ix,jx,kx)+au(i,kx)*apo(ix,jx,i) do 82 ix=1,3 do 82 jx=1,3 do 82 kx=1,3 apo(ix,jx,kx)=0.0d0 do 82 i=1,3 82 apo(ix,jx,kx)=apo(ix,jx,kx)+au(i,jx)*y(ix,i,kx) do 83 ix=1,3 do 83 jx=1,3 do 83 kx=1,3 y(ix,jx,kx)=0.0d0 do 83 i=1,3 83 y(ix,jx,kx)=y(ix,jx,kx)+au(i,ix)*apo(i,jx,kx) apo=y 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,IC,iwr real*8 ALPHA(3,3),G(3,3),A(3,3,3),c(3),SIGN,SUM 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 G(IA,1)=G(IA,1)+SIGN*0.5d0*(c(2)*ALPHA(IA,3)-c(3)*ALPHA(IA,2)) G(IA,2)=G(IA,2)+SIGN*0.5d0*(c(3)*ALPHA(IA,1)-c(1)*ALPHA(IA,3)) 2 G(IA,3)=G(IA,3)+SIGN*0.5d0*(c(1)*ALPHA(IA,2)-c(2)*ALPHA(IA,1)) IF(ICO.EQ.3)G=0.0d0 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 READRXS(N,R,q) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) real*8 R(3,N) integer*4 N,i,q(*) open(2,file='SMALL.X',status='old') READ(2,*) read(2,*)N do 1 i=1,N 1 READ(2,*)q(i),R(1,i),R(2,i),R(3,i) close(2) write(6,*)'SMALL.X read' 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,status='old') 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,status='old') 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,status='old') 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 subroutine wcos(cg,fg,cs,ic,ir,is,lext,ue,xe,ye, 1we,ur,xr,yr,wr,nmol,ig) implicit none real*8 cg(3,3),fg(3,3),cs(3,3),ue(*),xe(*),ye(*),am,an, 1ur(*),xr(*),yr(*),we(*),wr(*),ic,ir,z0,f integer*4 is,i,j,ii,ie(4),i2,nmol,ig real*8,allocatable::te(:) data ie/81,81,81,243/ character*2 ce(4) data ce/'W:','U:','X:','Y:'/ character*7 de(2) data de/' all:',' close:'/ logical lext f=dble(nmol*ig) z0=0.0d0 if(ic.gt.z0)write(6,600)' close',((cg(i,j)/ic,i=1,3),j=1,3) 600 format(a6,9f6.3) if(ir.gt.z0)write(6,600)' far',((fg(i,j)/ir,i=1,3),j=1,3) if(is.ne.0)write(6,600)' space',((cs(i,j)/dble(is),i=1,3),j=1,3) if(lext.and.ir.gt.z0.and.ic.gt.z0)then allocate(te(243)) do 1 i2=1,2 do 1 ii=1,4 an=0.0d0 am=0.0d0 do 2 i=1,ie(ii) if(i2.eq.1)then if(ii.eq.1)te(i)=we(i)/f if(ii.eq.2)te(i)=ue(i)/f if(ii.eq.3)te(i)=xe(i)/f if(ii.eq.4)te(i)=ye(i)/f else if(ii.eq.1)te(i)=wr(i)/f if(ii.eq.2)te(i)=ur(i)/f if(ii.eq.3)te(i)=xr(i)/f if(ii.eq.4)te(i)=yr(i)/f endif 2 continue do 3 i=1,ie(ii) an=an+te(i)**2 3 if(dabs(te(i)).gt.am)am=dabs(te(i)) an=dsqrt(an) 1 write(6,601)de(i2),ce(ii),an,am 601 format(a7,a2,' norm ',e10.4,', max ',e10.4) endif return end subroutine wrc(g,np,rmin,dr,nmol,ig,cc,icc,ci,lper,perx,pery,perz) implicit none real*8 g(*),rmin,dr,r,pi,cc(np,3,3),f,ci(np,3),perx,pery,perz,ra logical lper integer*4 np,i,nmol,icc(*),ix,jx,ig character*1 AB(3) data AB/'A','B','C'/ c average density: ra=1.0d0 if(lper)ra=dble(nmol)/perx/pery/perz c divide by dumber of pairs, number of geometries if(nmol.gt.0.and.ig.gt.0)ra=ra*dble(nmol-1)*dble(ig)/2.0d0 pi=4.0d0*datan(1.0d0) do 3 i=1,np if(icc(i).gt.0)then ci(i,1)=ci(i,1)/dble(icc(i)) ci(i,2)=ci(i,2)/dble(icc(i)) ci(i,3)=ci(i,3)/dble(icc(i)) endif 3 g(i)=g(i)/ra c write open(9,file='RDFM.PRN') r=rmin-dr do 1 i=1,np r=r+dr 1 write(9,900)r,g(i) 900 format(2f12.6) close(9) write(6,*)' RDFM.PRN' do 5 ix=1,3 open(9,file=AB(ix)//'.PRN') r=rmin-dr do 6 i=1,np r=r+dr 6 write(9,900)r,ci(i,ix) close(9) write(6,*)AB(ix)//'.PRN' do 5 jx=1,3 open(9,file=AB(ix)//AB(jx)//'.PRN') r=rmin-dr do 4 i=1,np r=r+dr f=cc(i,ix,jx) if(icc(i).gt.0)f=f/dble(icc(i)) 4 write(9,900)r,acos(f)*180.0d0/pi close(9) 5 write(6,*)AB(ix)//AB(jx)//'.PRN' return end subroutine incr(g,z,r,nats,nmol,lper,perx,pery,perz,dr,np,rmin, 1rlim,cg,fg,ic,ir,cc,icc,cs,ics,ci,ue,xe,ye,ur,xr,yr,lext, 1we,wr,drec,nrec,rrec,arec) implicit none logical lper,lext real*8 g(*),r(*),perx,pery,perz,dr,d,pi,rmin,ci(np,3),ai(3), 1xi,u(3,3),A(3,3),e(3),av(3),bv(3),sp,p(3),ma,tm,aj(3),rij(3), 1cg(3,3),fg(3,3),rlim,cv(3),cc(np,3,3),ct(3,3),cs(3,3),rj(3), 1bi(3),gi(3),bj(3),cj(3),ue(81),xe(81),ye(243),ri(3),arec(*), 1ur(81),xr(81),yr(243),wr(*),we(*),ic,ir,rrec(*),drec(*),qa integer*4 np,ni,nats,nmol,i,ix,ia,j,jx,IERR,icc(*),ics, 1MENDELEV,z(*),i1,i2,ii,nrec,i3 parameter (MENDELEV=89) real*8,allocatable::c(:,:),au(:,:,:) real*4 mass(MENDELEV) data mass/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ allocate(c(nmol,3),au(nmol,3,3)) c=0.0d0 pi=4.0d0*datan(1.0d0) IERR=0 do 4 i=1,nmol c mass tm, mass center c: tm=0.0d0 do 3 ia=1,nats i1=nats*(i-1)+ia ma=dble(mass(z(i1))) tm=tm+ma do 3 ix=1,3 3 c(i,ix)=c(i,ix)+ma*r(ix+3*(i1-1)) do 2 ix=1,3 2 c(i,ix)=c(i,ix)/tm c p = Qi xi : dipole moment: c u = Qi xi yi: quadrupole moment: p=0.0d0 u=0.0d0 do 5 ia=1,nats i1=nats*(i-1)+ia qa=z(i1) do 5 ix=1,3 xi=qa*(r(ix+3*(i1-1))-c(i,ix)) p(ix)=p(ix)+xi do 5 jx=1,3 5 u(ix,jx)=u(ix,jx)+xi*(r(jx+3*(i1-1))-c(i,jx)) c principal axes CALL TRED12(3,u,A,e,2,IERR) do 6 ix=1,3 av(ix)=A(ix,1) 6 bv(ix)=A(ix,2) c check that a-axis points along p if(sp(av,p).lt.0.0d0)then av(1)=-av(1) av(2)=-av(2) av(3)=-av(3) endif c check that b-axis points along p if(sp(bv,p).lt.0.0d0)then bv(1)=-bv(1) bv(2)=-bv(2) bv(3)=-bv(3) endif c c = a x b: call vp(av,bv,cv) ics=ics+1 do 4 ix=1,3 c save axes' vectors: au(i,ix,1)=av(ix) au(i,ix,2)=bv(ix) au(i,ix,3)=cv(ix) c average cosinus between space axis: cs(ix,1)=cs(ix,1)+av(ix) cs(ix,2)=cs(ix,2)+bv(ix) 4 cs(ix,3)=cs(ix,3)+cv(ix) c positions of two molecules: do 1 i=1,nmol do 14 ix=1,3 ai(ix)=au(i,ix,1) bi(ix)=au(i,ix,2) 14 gi(ix)=au(i,ix,3) do 1 j=i+1,nmol do 15 ix=1,3 aj(ix)=au(j,ix,1) bj(ix)=au(j,ix,2) 15 cj(ix)=au(j,ix,3) rij(1)=c(i,1)-c(j,1) rij(2)=c(i,2)-c(j,2) rij(3)=c(i,3)-c(j,3) c c calculate distance and correct to periodicity: call dp(d,rij,lper,perx,pery,perz) c projections to axes of i: rj(1)=rij(1)*ai(1)+rij(2)*ai(2)+rij(3)*ai(3) rj(2)=rij(1)*bi(1)+rij(2)*bi(2)+rij(3)*bi(3) rj(3)=rij(1)*gi(1)+rij(2)*gi(2)+rij(3)*gi(3) c projections to axes of j: ri(1)=-rij(1)*aj(1)-rij(2)*aj(2)-rij(3)*aj(3) ri(2)=-rij(1)*bj(1)-rij(2)*bj(2)-rij(3)*bj(3) ri(3)=-rij(1)*cj(1)-rij(2)*cj(2)-rij(3)*cj(3) c cosinuses between moment axes do 9 ix=1,3 do 9 jx=1,3 9 ct(ix,jx)=au(i,1,ix)*au(j,1,jx) 1 +au(i,2,ix)*au(j,2,jx) 1 +au(i,3,ix)*au(j,3,jx) if(d.lt.rlim)then c close molecules do 7 ix=1,3 do 7 jx=1,3 7 cg(ix,jx)=cg(ix,jx)+ct(ix,jx) ic=ic+1.0d0 if(lext)call addue(ur,xr,yr,wr,ct,ri,rj) else c far molecules do 8 ix=1,3 do 8 jx=1,3 8 fg(ix,jx)=fg(ix,jx)+ct(ix,jx) ir=ir+1.0d0 endif if(lext)call addue(ue,xe,ye,we,ct,ri,rj) c radial distribution: ni=nint((d-rmin)/dr)+1 if(ni.ge.1.and.ni.le.np)then c function: g(ni)=g(ni)+1.0d0/(4.0d0*pi*d**2*dr) c cosinus: do 10 ix=1,3 do 10 jx=1,3 10 cc(ni,ix,jx)=cc(ni,ix,jx)+ct(ix,jx) c chirality index call vp(rij,aj,cv) ci(ni,1)=ci(ni,1)+sp(ai,cv)/dsqrt(sp(rij,rij)) call vp(rij,bj,cv) ci(ni,2)=ci(ni,2)+sp(bi,cv)/dsqrt(sp(rij,rij)) call vp(rij,cj,cv) ci(ni,3)=ci(ni,3)+sp(gi,cv)/dsqrt(sp(rij,rij)) icc(ni)=icc(ni)+1 endif c record nrec closest pairs do 12 ii=nrec,1,-1 if(d.lt.drec(ii))then drec(ii)=d arec(ii)=ct(1,1) do 16 ia=1,nats i1=nats*(i-1)+ia i2=nats*(j-1)+ia do 16 ix=1,3 i3=ix+3*(ia-1) rrec(i3+3*nats*(ii-1) )=r(ix+3*(i1-1)) 16 rrec(i3+3*nats*(ii-1)+nrec*3*nats)=r(ix+3*(i2-1)) call rorder(nats,nrec,drec,rrec,arec) goto 1 endif 12 continue 1 continue return end subroutine rorder(nats,nrec,drec,rrec,arec) implicit none integer*4 nats,nrec,ir,irp,ii,ia,ix real*8 drec(*),rrec(*),t,arec(*) do 1 ir=1,nrec do 1 irp=ir+1,nrec if(drec(irp).lt.drec(ir))then t=drec(ir) drec(ir)=drec(irp) drec(irp)=t t=arec(ir) arec(ir)=arec(irp) arec(irp)=t do 2 ia=1,nats do 2 ix=1,3 ii=ix+3*(ia-1) t=rrec(ii+3*nats*(ir-1)) rrec(ii+3*nats*(ir-1))=rrec(ii+3*nats*(irp-1)) rrec(ii+3*nats*(irp-1))=t ii=ii+nrec*3*nats t=rrec(ii+3*nats*(ir-1)) rrec(ii+3*nats*(ir-1))=rrec(ii+3*nats*(irp-1)) 2 rrec(ii+3*nats*(irp-1))=t endif 1 continue return end subroutine vp(a,b,c) implicit none real*8 a(3),b(3),c(3) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end function sp(a,b) implicit none real*8 a(3),b(3),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT integer*4(I-N) DIMENSION A(N,N),Z(N,N),D(N) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) JP1 = J + 1 IF (L .LT. JP1) GO TO 220 DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE HH = F / (H + H) DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) 290 D(I) = H 300 CONTINUE 320 D(1) = 0.0D0 E(1) = 0.0D0 IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 DO 360 J = 1, L G = 0.0D0 DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE 500 CONTINUE 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT integer*4(I-N) DIMENSION Z(N,N),D(N),E(N) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS IERR = 0 IF (N .EQ. 1) GO TO 1001 DO 100 I = 2, N 100 E(I-1) = E(I) F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 110 CONTINUE 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) DO 140 I = L1, N 140 D(I) = D(I) - H F = F + H P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE 200 CONTINUE E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE GO TO 1001 1000 IERR = L 1001 RETURN END subroutine cre(x,perx) implicit none real*8 x,perx if(dabs(x+perx).lt.dabs(x))then x=x+perx else if(dabs(x-perx).lt.dabs(x))x=x-perx endif return end subroutine dp(d,x,lper,perx,pery,perz) implicit none logical lper real*8 d,x(3),perx,pery,perz if(lper)then call cre(x(1),perx) call cre(x(2),pery) call cre(x(3),perz) endif d=dsqrt(x(1)**2+x(2)**2+x(3)**2) return end subroutine readopt(nats,np,rmin,rmax,rlim,lper,perx,pery,perz, 1irep,lext,lpol,lrest,nrec) implicit none logical lper,lext,lpol,lrest integer*4 nats,np,irep,nrec real*8 rmin,rmax,rlim,perx,pery,perz character*4 key irep=100 nats=0 np=100 rmin=0.0d0 rmax=20.0d0 rlim=5.0d0 lper=.false. perx=0.0d0 pery=0.0d0 perz=0.0d0 lext=.true. lrest=.false. lpol=.false. c record nrec smallest distances: nrec=10 open(8,file='RDFU.OPT',status='old') 1 read(8,80,end=88,err=88)key 80 format(a4) if(key.eq.'NATS')read(8,*)nats if(key.eq.'RMIN')read(8,*)rmin if(key.eq.'RMAX')read(8,*)rmax if(key.eq.'RLIM')read(8,*)rlim if(key.eq.'PERX')read(8,*)perx if(key.eq.'PERY')read(8,*)pery if(key.eq.'PERZ')read(8,*)perz if(key.eq.'LPER')read(8,*)lper if(key.eq.'LEXT')read(8,*)lext if(key.eq.'LPOL')read(8,*)lpol if(key.eq.'IREP')read(8,*)irep if(key.eq.'REST')read(8,*)lrest if(key.eq.'NREC')read(8,*)nrec if(key.eq.'NPOI')read(8,*)np goto 1 88 close(8) if(nats.eq.0)call report('NATS = 0') write(6,600)nats 600 format(' NATS = ',I6) return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine addue(ue,xe,ye,we,ct,ri,rj) implicit none real*8 ue(*),xe(*),ye(*),ct(3,3),ri(*),rj(*),we(*) integer*4 ii,ix,jx,kx,lx,mx,w,f,iip c U_abde: ii=0 do 1 ix=1,3 do 1 jx=1,3 do 1 kx=1,3 do 1 lx=1,3 ii=ii+1 if(ix.eq.jx.and.kx.eq.lx)ue(ii)=ue(ii)-2.0d0/3.0d0 1 ue(ii)=ue(ii)+ct(ix,kx)*ct(jx,lx)+ct(kx,ix)*ct(lx,jx) c sum(j>i) ( if(ix.eq.jx) 1xe(ii)=xe(ii)+(ct(2,kx)*ct(3,lx)-ct(3,kx)*ct(2,lx))*rj(1) 1 +(ct(3,kx)*ct(1,lx)-ct(1,kx)*ct(3,lx))*rj(2) 1 +(ct(1,kx)*ct(2,lx)-ct(2,kx)*ct(1,lx))*rj(3) 1 +(ct(kx,2)*ct(lx,3)-ct(kx,3)*ct(lx,2))*ri(1) 1 +(ct(kx,3)*ct(lx,1)-ct(kx,1)*ct(lx,3))*ri(2) 1 +(ct(kx,1)*ct(lx,2)-ct(kx,2)*ct(lx,1))*ri(3) c EPS_awf rw w=iip(ix) f=iip( w) xe(ii)=xe(ii)+ct(jx,kx)*(ct(f,lx)*rj(w)-ct(w,lx)*rj(f)) 1 +ct(kx,jx)*(ct(lx,f)*ri(w)-ct(lx,w)*ri(f)) w=iip(jx) f=iip( w) c EPS_bwf rw 2 xe(ii)=xe(ii)+ct(ix,kx)*(ct(w,lx)*rj(f)-ct(f,lx)*rj(w)) 1 +ct(kx,ix)*(ct(lx,w)*ri(f)-ct(lx,f)*ri(w)) c Y_abcde: ii=0 do 3 ix=1,3 do 3 jx=1,3 do 3 kx=1,3 do 3 lx=1,3 do 3 mx=1,3 ii=ii+1 if(kx.eq.lx)then w=iip(jx) f=iip( w) if(jx.eq.w)ye(ii)=ye(ii)-2.0d0*(ct(f,mx)+ct(mx,f)) if(jx.eq.f)ye(ii)=ye(ii)+2.0d0*(ct(w,mx)+ct(mx,w)) endif w=iip(jx) f=iip( w) ye(ii)=ye(ii)+3.0d0*( 1(ct(ix,lx)*ct(w,kx)+ct(w,lx)*ct(ix,kx))*ct(f,mx)- 1(ct(ix,lx)*ct(f,kx)+ct(f,lx)*ct(ix,kx))*ct(w,mx)+ 1(ct(lx,ix)*ct(kx,w)+ct(lx,w)*ct(kx,ix))*ct(mx,f)- 1(ct(lx,ix)*ct(kx,f)+ct(lx,f)*ct(kx,ix))*ct(mx,w)) w=iip(ix) f=iip( w) 3 ye(ii)=ye(ii)+ 1(ct(jx,lx)*ct(w,kx)+ct(w,lx)*ct(jx,kx))*ct(f,mx)- 1(ct(jx,lx)*ct(f,kx)+ct(f,lx)*ct(jx,kx))*ct(w,mx)+ 1(ct(lx,jx)*ct(kx,w)+ct(lx,w)*ct(kx,jx))*ct(mx,f)- 1(ct(lx,jx)*ct(kx,f)+ct(lx,f)*ct(kx,jx))*ct(mx,w) c = sum(j>i) (3 eps_wxb + 3 eps_fxb c +eps_wxa +eps_fxa -2 eps_bxa del_dc] c W_abde: ii=0 do 4 ix=1,3 do 4 jx=1,3 do 4 kx=1,3 do 4 lx=1,3 ii=ii+1 if(ix.eq.jx.and.kx.eq.lx)we(ii)=we(ii)-7.0d0/3.0d0 4 we(ii)=we(ii)+6.0d0*ct(ix,lx)*ct(jx,kx)+ct(jx,lx)*ct(ix,kx) c 6 +-(7/3) d_ab d_ed return end function iip(i) implicit none integer*4 iip,i if(i.lt.3)then iip=i+1 else iip=1 endif return end subroutine aar(ig,nmol,we,ue,xe,ye) implicit none real*8 f,we(81),ue(81),xe(81),ye(243) integer*4 i,ig,nmol f=dble(ig*nmol) do 1 i=1,81 we(i)=we(i)/f ue(i)=ue(i)/f 1 xe(i)=xe(i)/f do 2 i=1,243 2 ye(i)=ye(i)/f open(9,file='AV.SCR',form='unformatted') write(9)we,ue,xe,ye close(9) write(6,*)' Averages written to AV.SCR.' return end