program rdfutest implicit none integer*4 nat,ig,ia,n,nats,nmol,np,irep,ics,nrec,pkind parameter (pkind=213) 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,ltest real*8,allocatable::r(:),g(:),cc(:,:,:),ci(:,:),ue(:),xe(:),ye(:), 1ur(:),xr(:),yr(:),we(:),wr(:),rrec(:),drec(:),arec(:),alave(:) 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,ltest) 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) if(ltest)then allocate(alave(pkind)) alave=0.d0 endif 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) if(ltest) call avepol(nmol,nat,nats,po,gpo,apo,r,z, 1alave,perx,pery,perz,lper) 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) if(ltest)call wrave(pkind,alave,nmol,ig,po,gpo,apo,2) 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) if(ltest)call wrave(pkind,alave,nmol,ig,po,gpo,apo,1) end subroutine avepol(NO,NAT,nats,po,gpo,apo,r,q, 1alave,perx,pery,perz,lper) implicit none integer*4 IX,JX,KX,NO,NAT,nats,IO,JO,ii,i1,ME,q(*),ia,IERR real*8 m,pi(3,3),gi(3,3),ai(3,3,3),alave(*),r(*),e(3),A(3,3), 1dk,perx,pery,av(3),bv(3),cv(3),au(3,3),gpo(3,3),sp, 1perz,pj(3,3),gj(3,3),aj(3,3,3),ddp,p(3),u(3,3),po(3,3), 1apo(3,3,3),xi logical lper real*8,allocatable::coo(:,:),opol(:,:,:),gopol(:,:,:), 1aopol(:,:,:,:) parameter (ME=89) real*4 mass(ME) 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/ IERR=0 allocate(coo(NO,3),opol(NO,3,3),gopol(NO,3,3),aopol(NO,3,3,3)) c precalculate variables for each molecule: DO 100 IO=1,NO c mass: m=0.0d0 do 103 ia=1,nats 103 m=m+dble(mass(q(nats*(IO-1)+ia))) c mass center: do 101 ix=1,3 coo(IO,ix)=0.0d0 do 102 ia=1,nats i1=nats*(IO-1)+ia 102 coo(IO,ix)=coo(IO,ix)+r(ix+3*(i1-1))*dble(mass(q(i1))) 101 coo(IO,ix)=coo(IO,ix)/m c p = Qi xi : dipole moment: c u = Qi xi yi: quadrupole: p=0.0d0 u=0.0d0 do 104 ia=1,nats i1=nats*(IO-1)+ia do 104 ix=1,3 xi=dble(q(i1)) * (r(ix+3*(i1-1))-coo(IO,ix)) p(ix)=p(ix)+xi do 104 jx=1,3 104 u(ix,jx)=u(ix,jx)+xi*(r(jx+3*(i1-1))-coo(IO,jx)) c principal axes: CALL TRED12(3,u,A,e,2,IERR) do 105 ix=1,3 av(ix)=A(ix,1) 105 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 106 ix=1,3 au(ix,1)=av(ix) au(ix,2)=bv(ix) 106 au(ix,3)=cv(ix) c Polarizabilities in laboratory axes: call tt2( po, opol,au,NO,IO) call tt2(gpo,gopol,au,NO,IO) 100 call tt3(apo,aopol,au,NO,IO) 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 107 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 107 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 107 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 107 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 subroutine tt2(po,opol,au,NO,IO) implicit none integer*4 ix,jx,ii,jj,NO,IO real*8 t2(3,3),po(3,3),opol(NO,3,3),au(3,3) do 107 ix=1,3 do 107 jx=1,3 t2(ix,jx)=0.0d0 do 107 ii=1,3 107 t2(ix,jx)=t2(ix,jx)+au(ix,ii)*po(ii,jx) do 108 ix=1,3 do 108 jx=1,3 opol(IO,ix,jx)=0.0d0 do 108 jj=1,3 108 opol(IO,ix,jx)=opol(IO,ix,jx)+au(jx,jj)*t2(ix,jj) return end subroutine tt3(po,opol,au,NO,IO) implicit none integer*4 ix,jx,kx,ii,jj,kk,NO,IO real*8 t2(3,3,3),po(3,3,3),opol(NO,3,3,3),au(3,3),t3(3,3,3) do 107 ix=1,3 do 107 jx=1,3 do 107 kx=1,3 t2(ix,jx,kx)=0.0d0 do 107 ii=1,3 107 t2(ix,jx,kx)=t2(ix,jx,kx)+au(ix,ii)*po(ii,jx,kx) do 109 ix=1,3 do 109 jx=1,3 do 109 kx=1,3 t3(ix,jx,kx)=0.0d0 do 109 jj=1,3 109 t3(ix,jx,kx)=t3(ix,jx,kx)+au(jx,jj)*t2(ix,jj,kx) do 108 ix=1,3 do 108 jx=1,3 do 108 kx=1,3 opol(IO,ix,jx,kx)=0.0d0 do 108 kk=1,3 108 opol(IO,ix,jx,kx)=opol(IO,ix,jx,kx)+au(kx,kk)*t2(ix,jx,kk) return 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: p=0.0d0 do 10 ia=1,nats do 10 ix=1,3 10 p(ix)=p(ix)+dble(q(ia))*(RS(ix,ia)-c(ix)) c u = Qi xi yi: quadrupole: u=0.0d0 do 5 ia=1,nats do 5 ix=1,3 xi=dble(q(ia))*(RS(ix,ia)-c(ix)) 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 11 ix=1,3 au(ix,1)=av(ix) au(ix,2)=bv(ix) 11 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,ltest) implicit none logical lper,lext,lpol,lrest,ltest 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. ltest=.true. 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.'TEST')read(8,*)ltest 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) if(ltest.and..not.lpol)call report('Switch on LPOL for TEST!') 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 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 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 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 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 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 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 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 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 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