c program to provide phonon frequencies c this program assumes that the input file is c from "config.text", which is a text-file equivalent c of Li Ju's "config" output file (see the read section below) c c c **** due to laziness, I have not yet changed c this code to be general for defects for which n is not c conserved, such as vacancies. In that case, the number 648 c must be changed to 645. (8 places in the code). ******* c implicit real(a-h,o-z) real lloopxx,lloopxy,nloopxx,nloopxy, a lloop2xx,lloop2xy,nloop2xx,nloop2xy, a nloop,lloopyx,nloopyx,lloop2yx, a nloop2yx,lloopyy,nloopyy,lloop2yy, a nloop2yy,lloopzz,lloopxz,nloopzz, a nloopxz,lloop2zz,lloop2xz,nloop2zz, a nloop2xz,lloopyz, lloopzx, lloopzy, a nloopyz, nloopzx, nloopzy, lloop2yz, a lloop2zx, lloop2zy, nloop2yz,nloop2zx, a nloop2zy, a nloop2,lloop3xx,lloop3xy,lloop3xz, a lloop3yx,lloop3yy,lloop3yz,lloop3zx, a lloop3zy,lloop3zz,lloop4xx,lloop4xy, a lloop4xz,lloop4yx,lloop4yy,lloop4yz, a lloop4zx,lloop4zy,lloop4zz, a nloop3xx,nloop3xy,nloop3xz,nloop3yx, a nloop3yy,nloop3yz,nloop3zx,nloop3zy, a nloop3zz,nloop4xx,nloop4xy,nloop4xz, a nloop4yx,nloop4yy,nloop4yz,nloop4zx, a nloop4zy,nloop4zz, a lloop5xx,lloop5yy,lloop5zz,lloop5xy, a lloop5xz,lloop5yx,lloop5yz,lloop5zx, a lloop5zy c PARAMETER ( nm=600, NX = 20 ,pip=3.1415926536d0) c parameter ( usigma=2.352127d0, ueps=2.167285d0 ) parameter ( usigma=1.85762, ueps=2.1672850 ) c parameter ( usigma=1.541972610, ueps=2.1672850 ) c c DIMENSION FX(NM),FY(NM),FZ(NM) DIMENSION PI(3,3) DIMENSION SX(NM), SY(NM), SZ(NM), VX(NM), VY(NM), VZ(NM), 1 AX(NM), AY(NM), AZ(NM), BX(NM), BY(NM), BZ(NM), 1 CX(NM), CY(NM), CZ(NM) DIMENSION H(3,3), HV(3,3), HA(3,3), HB(3,3), HC(3,3) DIMENSION HI(3,3), ntyp(nm), sqm(2) c dimension x(nm,nx),y(nm,nx),z(nm,nx),xt(nm),yt(nm),zt(nm), a r2(nm,nx),r(nm,nx),list2(nx) dimension zet2(nm,nx,nx),zet3(nm,nx,nx),epp(nm),fc(nm), a fcp(nm),fcpp(nm) dimension zet12(nx,nx),zet13(nx,nx), a zet22(nx,nx),zet23(nx,nx),zet33(nx,nx) dimension rjk2(nx,nx),rjk(nx,nx),rjkx(nx,nx),rjky(nx,nx), a rjkz(nx,nx), a prjkx(nx,nx),prjky(nx,nx),prjkz(nx,nx) dimension pa(2), pb(2),pl1(2),pl2(2), pbt(2),pn(2),pc(2),pwr2(2) a ,pd(2), ph(2), pr(2), ps(2), c2d2(2), pwr(2),pwr1(2) dimension chij(nx),rlm1(nx), rlm2(nx),ca(nx),cb(nx), a cr(nx),cs(nx) dimension h00(3,3),cpb(3,4) dimension zet(nx),szet1(nx),szet11(nx),sftm2x(nx),sftm2y(nx), a sftm2z(nx),sftm3x(nx),sftm3y(nx),sftm3z(nx) dimension pf(nx),pf2(nx),pf3(nm,nx),frst(nx),frst2(nx), a ff11(nm,nx), strs21(nx),strs22(nx),strs23(nx),strs24(nx), a strs25(nx),strs26(nx),strs31(nx),strs32(nx),strs33(nx), a strs34(nx),strs35(nx),strs36(nx),ff22(nm,nx,nx), a ff23(nm,nx,nx),ff13(nm,nx,nx),ff1(nm,nx), a ff12(nm,nx,nx),ff3(nm,nx,nx),ff2(nm,nx,nx),ff33(nm,nx,nx) dimension c11tm2(nx),c11tm3(nx),c11tm4(nx),c11tm5(nx),c11tm6(nx) dimension c11tm22(nx),c11tm33(nx),c11tm23(nx),c11tm32(nx) dimension c12tm2(nx),c12tm3(nx),c12tm4(nx),c12tm5(nx),c12tm6(nx) dimension c12tm22(nx),c12tm33(nx),c12tm23(nx),c12tm32(nx) dimension c44tm2(nx),c44tm3(nx),c44tm4(nx),c44tm5(nx),c44tm6(nx) dimension c44tm22(nx),c44tm33(nx),c44tm23(nx),c44tm32(nx) dimension nn(nm),idpart(nm,nx),negno(nm,nm) dimension DXX(nm,nm), DXY(nm,nm),DYX(nm,nm),DYY(nm,nm), a DZZ(nm,nm), DXZ(nm,nm),DYZ(nm,nm),DZX(nm,nm) a , DZY(nm,nm) c C real HESS(648,648) real eval(648) data mm,nmax/648,648/ c open(unit=98,file='config.text') read(98,*) n read(98,'(f15.8)')h do 7 ix = 1, n read(98,211) ix, sx(ix),sy(ix),sz(ix),ntyp(ix) 7 continue 211 format(1x,i4,3(1x,f15.8),2x,i1,3(1x,i3)) c do 41 i = 1,3 do 42 j = 1,3 H(i,j) = H(i,j)/1.85762 42 continue 41 continue c write(*,*) n write(*,*) h c C LI JU DEFINES SI = 1 AND C = 2, WHICH IS THE OPPOSITE OF MEIJIE C SO, SWITCH THEM C DO 43 I = 1,N IF (NTYP(I).EQ.1) THEN NTYP(I) = 2 ELSE NTYP(I) = 1 ENDIF 43 CONTINUE c c test to make sure reading okay c do 47 ij = 1,n c write(71,*) ij, sx(ij),sy(ij),sz(ij),ntyp(ij) c 47 continue c s=1.0 sv = 0.0 sa = 0.0 sb = 0.0 sc = 0.0 c np = n sqm(2)=1.0 sqm(1)=sqrt(12.0/27.97720) ratl_cut = 1.0 c pa(2)=1.8308e+3/ueps pb(2)=4.7118e+2/ueps pl1(2)=2.4799*usigma pl2(2)=1.7322*usigma pn(2)=7.8734e-1 pbt(2)=(1.1e-6)**pn(2) pc(2)=(1.0039e+5)**2. pd(2)=(1.6217e+1)**2. ph(2)=-5.9825e-1 pr(2)=2.70/usigma*ratl_cut ps(2)=3.00/usigma*ratl_cut c2d2(2)=1.0+pc(2)/pd(2) pwr(2)=-1.0/(2.0*pn(2)) pwr1(2)=pwr(2)-1.0 pwr2(2)=pn(2)-1.0 c 1st set of paramters for C. PRB 39,5566(1989) c pa(1)=1.3936d3/ueps c pb(1)=3.4670d2/ueps c pl1(1)=3.4879d0*usigma c pl2(1)=2.2119d0*usigma c pn(1)=7.2751d-1 c pbt(1)=(1.5724d-7)**pn(1) c pc(1)=(3.8049d4)**2. c pd(1)=(4.384d0)**2. c ph(1)=-5.7058d-1 c pr(1)=1.93d0/usigma*ratl_cut c ps(1)=2.13d0/usigma*ratl_cut c c2d2(1)=1.d0+pc(1)/pd(1) c pwr(1)=-1.d0/(2.d0*pn(1)) c pwr1(1)=pwr(1)-1.d0 c pwr2(1)=pn(1)-1.d0 c 2nd set of paramters for C. PRL 64,1757(1990) pa(1)=1.5448e+3/ueps pb(1)=3.8963e+2/ueps pl1(1)=3.4653*usigma pl2(1)=2.3064*usigma pn(1)=9.9054e-1 pbt(1)=(4.1612e-6)**pn(1) pc(1)=(1.9981e+4)**2. pd(1)=(7.034)**2. ph(1)=-3.9953e-1 pr(1)=1.80/usigma*ratl_cut ps(1)=2.10/usigma*ratl_cut c2d2(1)=1.0+pc(1)/pd(1) pwr(1)=-1.0/(2.0*pn(1)) pwr1(1)=pwr(1)-1.0 pwr2(1)=pn(1)-1.0 c pssic=2.56/usigma*ratl_cut prsic=2.36/usigma*ratl_cut c c************************** c Find square of cut-off c Zero force and energy arrays c do 10 i=1,np fx(i)=0.0 fy(i)=0.0 fz(i)=0.0 epp(i)=0.0 10 continue c goto 55 c if(mod(istep,iprint).eq.0) then c11p = 0.0 c12p = 0.0 c44p = 0.0 c11s = 0.0 c12s = 0.0 c44s = 0.0 c11t = 0.0 c12t = 0.0 c44t = 0.0 do 54 j=1,4 cpb(1,j)=0.0 cpb(2,j)=0.0 cpb(3,j)=0.0 54 continue c endif 55 continue c c zero the stress components c pxx =0.00 pxy =0.00 pxz =0.00 pyx =0.00 pyy =0.00 pyz =0.00 pzx =0.00 pzy =0.00 pzz =0.00 c do 57 i = 1,np do 58 j = 1,np negno(i,j) = 0 58 continue 57 continue c do 80 i=1,nx 80 list2(i)=0 c c loop over atoms c _____________________ c do 500 ipt=1,np n1=ntyp(ipt) c i = ipt c SXI = SX(IPT) SYI = SY(IPT) SZI = SZ(IPT) c c c ================ c first inner loop c ================ c c Run through all particles, determining those inside cutoff c neg = 0 do 600 jpt = 1,np if(ipt.eq.jpt) go to 600 c n2=ntyp(jpt) if(n1.eq.n2) then RSMAX = ps(n1)*ps(n2) else rsmax=pssic**2 endif c c distance between atom i and close image of atom j c and r(ij) c SXIJ = SX(JPT) - SXI SYIJ = SY(JPT) - SYI SZIJ = SZ(JPT) - SZI c SXIJ = SXIJ - ANINT ( SXIJ ) SYIJ = SYIJ - ANINT ( SYIJ ) SZIJ = SZIJ - ANINT ( SZIJ ) c c XIJ = (H(1,1)*SXIJ + H(1,2)*SYIJ + H(1,3)*SZIJ)*usigma c YIJ = (H(2,1)*SXIJ + H(2,2)*SYIJ + H(2,3)*SZIJ)*usigma c ZIJ = (H(3,1)*SXIJ + H(3,2)*SYIJ + H(3,3)*SZIJ)*usigma c XIJ = H(1,1)*SXIJ + H(1,2)*SYIJ + H(1,3)*SZIJ YIJ = H(2,1)*SXIJ + H(2,2)*SYIJ + H(2,3)*SZIJ ZIJ = H(3,1)*SXIJ + H(3,2)*SYIJ + H(3,3)*SZIJ c c rij2 = xij**2. + yij**2. + zij**2. c c if(rij2.gt.rsmax) go to 600 c neg = neg + 1 rij = sqrt(rij2) x(i,neg) = xij y(i,neg) = yij z(i,neg) = zij xt(neg) = xij/rij yt(neg) = yij/rij zt(neg) = zij/rij r2(i,neg) = rij2 r(i,neg) = rij list2(neg) = jpt idpart(i,neg) = jpt negno(i,jpt) = neg if(ipt.eq.38) write (6,*) idpart(i,neg) c***************************************************************** c Below is added for Tersoff if(n1.eq.n2) then rlm1(neg)=pl1(n1) rlm2(neg)=pl2(n1) ca(neg)=pa(n1) cb(neg)=pb(n1) cr(neg)=pr(n1) cs(neg)=ps(n1) chij(neg)=1.0 c if(mod(istep,1000).eq.0) then c write(89,*) istep,ipt,jpt,rij,cr(neg),cs(neg) c endif else rlm1(neg)=(pl1(1)+pl1(2))/2.0 rlm2(neg)=(pl2(1)+pl2(2))/2.0 ca(neg)=sqrt(pa(1)*pa(2)) cb(neg)=sqrt(pb(1)*pb(2)) cr(neg)=prsic cs(neg)=pssic c for second set of parameter: chij=0.9972 c use Li Ju's parameter of 1.0086 c chij(neg)=1.00860 endif if(rij.lt.cr(neg)) then fc(neg)=1.0 fcp(neg)=0.0 fcpp(neg)=0.0 else fc(neg)=0.50*(1.0+cos(pip*(rij-cr(neg)) a /(cs(neg)-cr(neg)))) fcp(neg)=-0.50*pip*sin(pip*(rij-cr(neg)) a /(cs(neg)-cr(neg)))/(cs(neg)-cr(neg)) fcpp(neg)=(0.50-fc(neg))*pip*pip/(cs(neg)-cr(neg)) a /(cs(neg)-cr(neg)) c if(mod(istep,1000).eq.0) then c write(88,*) ipt,jpt,rij,fc(neg) c endif endif 600 continue c nn(i) = neg c write(6,*) ipt,neg,i,nn(i) c if (neg.gt. 4) then c write(6,*) 'i = ',ipt c write(6,*) 'neighbors = ',neg c endif c ================= c second inner loop c ================= c do 700 j = 1,neg jpt = list2(j) n22=ntyp(jpt) c c************************* c zero accumulators zet(j)=0.0 szet1(j)=0.0 szet11(j)=0.0 sftm2x(j)=0.0 sftm2y(j)=0.0 sftm2z(j)=0.0 sftm3x(j)=0.0 sftm3y(j)=0.0 sftm3z(j)=0.0 c************************** do 800 k = 1,neg if(k.eq.j) goto 800 kp = list2(k) c234567 rmu = (x(i,j)*x(i,k)+y(i,j)*y(i,k) a +z(i,j)*z(i,k))/(r(i,j)*r(i,k)) rjk2(j,k)=r2(i,j)+r2(i,k)-2.0*rmu*r(i,j)*r(i,k) rjk(j,k)=sqrt(rjk2(j,k)) rjkx(j,k)=x(i,j)-x(i,k) rjky(j,k)=y(i,j)-y(i,k) rjkz(j,k)=z(i,j)-z(i,k) prjkx(j,k)=rjkx(j,k)/rjk(j,k) prjky(j,k)=rjky(j,k)/rjk(j,k) prjkz(j,k)=rjkz(j,k)/rjk(j,k) c rmu1=rmu-ph(n1) rmu2=rmu1*rmu1 rmu3=pd(n1)+rmu2 gmu=c2d2(n1)-pc(n1)/rmu3 gmup=2.0*rmu1*pc(n1)/rmu3**2 c zet(j)=zet(j)+fc(k)*gmu com1=fcp(k)*gmu com2=fc(k)*gmup com3=fcpp(k)*gmu c tip1=rmu/r(i,j)-1.0/r(i,k) tip2=rmu/r(i,k)-1.0/r(i,j) tip3=rjk(j,k)/r(i,j)/r(i,k) zet1 = com2*(-tip1) szet1(j)=szet1(j)+zet1 zet2(i,j,k)=com1+com2*(-tip2) zet3(i,j,k)=com2*(-tip3) c goto 66 c if(mod(istep,iprint).eq.0) then gmupp=2.0*pc(n1)/rmu3**2 a - 8.0*pc(n1)*rmu2/rmu3**3 com4=fcp(k)*gmup com5=fc(k)*gmupp tip11=2.0*rmu/r2(i,j)-1.0/r(i,j)/r(i,k) tip12=rmu/r(i,j)/r(i,k)-1.0/r2(i,k)-1.0/r2(i,j) tip13=rjk(j,k)/r2(i,j)/r(i,k) tip22=2.0*rmu/r2(i,k)-1.0/r(i,j)/r(i,k) tip23=rjk(j,k)/r2(i,k)/r(i,j) tip33=-1.0/r(i,j)/r(i,k) zet11 = com5*tip1*tip1+com2*tip11 szet11(j) = szet11(j) + zet11 zet12(j,k) = com4*(-tip1)+com5*tip2*tip1+com2*tip12 zet13(j,k) = com5*tip3*tip1+com2*tip13 zet22(j,k) = com3+2.0*com4*(-tip2)+com5*tip2*tip2+com2*tip22 zet23(j,k) = com4*(-tip3)+com5*tip2*tip3+com2*tip23 zet33(j,k) = com5*tip3*tip3+com2*tip33 c endif 66 continue c sftm2x(j) = sftm2x(j) + zet2(i,j,k)*xt(k) sftm2y(j) = sftm2y(j) + zet2(i,j,k)*yt(k) sftm2z(j) = sftm2z(j) + zet2(i,j,k)*zt(k) sftm3x(j) = sftm3x(j) + zet3(i,j,k)*prjkx(j,k) sftm3y(j) = sftm3y(j) + zet3(i,j,k)*prjky(j,k) sftm3z(j) = sftm3z(j) + zet3(i,j,k)*prjkz(j,k) c 800 continue 700 continue c************************ c calculate the summation prefactor, constant for each force c pwr=-1/2n, pn=n, pwr2=n-1,pwr1=-1/2n-1 c do 704 j=1,neg jpt=list2(j) n22=ntyp(jpt) elm2=exp(-rlm2(j)*r(i,j)) elm2b=elm2*cb(j) elm1=exp(-rlm1(j)*r(i,j)) elm1a=elm1*ca(j) fcd1=fcp(j)-rlm1(j)*fc(j) fcd2=fcp(j)-rlm2(j)*fc(j) fcd1p=fcpp(j)-2.0*rlm1(j)*fcp(j)+rlm1(j)**2.*fc(j) fcd2p=fcpp(j)-2.0*rlm2(j)*fcp(j)+rlm2(j)**2.*fc(j) c fcd1=-rlm1(j)*fc(j) c fcd2=-rlm2(j)*fc(j) c fcd1p=-rlm1(j)*fcp(j)+rlm1(j)**2.*fc(j) c fcd2p=-rlm2(j)*fcp(j)+rlm2(j)**2.*fc(j) if(zet(j).ne.0.0) then bezet=1.0+ pbt(n1)*(zet(j)**pn(n1)) bij=(bezet**pwr(n1))*chij(j) c if(zet(j).lt.1.e-5) write(*,*) 'zet=',zet(j) bij1=0.50*pbt(n1)*(zet(j)**pwr2(n1)) bijp=bij1*(bezet**pwr1(n1))*chij(j) pf(j) = fc(j)*elm2b*bijp pf2(j) = fcd2*elm2b*bijp pf3(i,j)=pf(j)*(pwr2(n1)/zet(j)- a (1.+2.*pn(n1))*bij1/bezet) else bij=1.0 pf(j)=0.0 pf2(j)=0.0 pf3(i,j)=0.0 endif c frst(j) = fcd1*elm1a - bij*fcd2*elm2b frst2(j)=fcd1p*elm1a -bij*fcd2p*elm2b ff1(i,j)=frst(j)+pf(j)*szet1(j) ff11(i,j)=frst2(j)+2.0*pf2(j)*szet1(j) a +pf3(i,j)*szet1(j)**2+pf(j)*szet11(j) c c calculate potential energy per particle c vijr=fc(j)*elm1a vija=-fc(j)*bij*elm2b c write(*,*) ipt,jpt,vijr*ueps,vija*ueps,bij,fc(j) vij=vijr+vija epp(ipt)=epp(ipt) + vij/2.0 704 continue c do 701 j=1,neg jpt=list2(j) n22=ntyp(jpt) strs21(j) = 0.0 strs22(j) = 0.0 strs23(j) = 0.0 strs24(j) = 0.0 strs25(j) = 0.0 strs26(j) = 0.0 c strs31(j) = 0.0 strs32(j) = 0.0 strs33(j) = 0.0 strs34(j) = 0.0 strs35(j) = 0.0 strs36(j) = 0.0 c goto 77 c if(mod(istep,iprint).eq.0) then c11tm2(j) = 0.0 c11tm3(j) = 0.0 c11tm4(j) = 0.0 c11tm5(j) = 0.0 c11tm6(j) = 0.0 c12tm2(j) = 0.0 c12tm3(j) = 0.0 c12tm4(j) = 0.0 c12tm5(j) = 0.0 c12tm6(j) = 0.0 c44tm2(j) = 0.0 c44tm3(j) = 0.0 c44tm4(j) = 0.0 c44tm5(j) = 0.0 c44tm6(j) = 0.0 c11tm22(j) = 0.0 c12tm22(j) = 0.0 c44tm22(j) = 0.0 c11tm33(j) = 0.0 c12tm33(j) = 0.0 c44tm33(j) = 0.0 c11tm23(j) = 0.0 c12tm23(j) = 0.0 c44tm23(j) = 0.0 c11tm32(j) = 0.0 c12tm32(j) = 0.0 c44tm32(j) = 0.0 c endif 77 continue 701 continue c c three-body stress terms and three-body forces. c do 702 j=1,neg jpt=list2(j) n22=ntyp(jpt) do 840 k=1,neg kp=list2(k) if(k.eq.j) goto 840 ff2(i,j,k)=pf(j)*zet2(i,j,k) ff3(i,j,k)=pf(j)*zet3(i,j,k) ff12(i,j,k)=pf2(j)*zet2(i,j,k)+ a pf3(i,j)*szet1(j)*zet2(i,j,k) a +pf(j)*zet12(j,k) ff13(i,j,k)=pf2(j)*zet3(i,j,k)+ a pf3(i,j)*szet1(j)*zet3(i,j,k) a +pf(j)*zet13(j,k) ff22(i,j,k)=pf3(i,j)*zet2(ipt,j,k)**2+pf(j)*zet22(j,k) ff23(i,j,k)=pf3(i,j)*zet2(ipt,j,k)*zet3(ipt,j,k) a +pf(j)*zet23(j,k) ff33(i,j,k)=pf3(i,j)*zet3(i,j,k)**2+pf(j)*zet33(j,k) c strs21(j)=strs21(j)+xt(k)*x(i,k)*ff2(i,j,k) strs22(j)=strs22(j)+xt(k)*y(i,k)*ff2(i,j,k) strs23(j)=strs23(j)+xt(k)*z(i,k)*ff2(i,j,k) strs24(j)=strs24(j)+yt(k)*y(i,k)*ff2(i,j,k) strs25(j)=strs25(j)+yt(k)*z(i,k)*ff2(i,j,k) strs26(j)=strs26(j)+zt(k)*z(i,k)*ff2(i,j,k) c strs31(j)=strs31(j)+prjkx(j,k)*rjkx(j,k)*ff3(i,j,k) strs32(j)=strs32(j)+prjkx(j,k)*rjky(j,k)*ff3(i,j,k) strs33(j)=strs33(j)+prjkx(j,k)*rjkz(j,k)*ff3(i,j,k) strs34(j)=strs34(j)+prjky(j,k)*rjky(j,k)*ff3(i,j,k) strs35(j)=strs35(j)+prjky(j,k)*rjkz(j,k)*ff3(i,j,k) strs36(j)=strs36(j)+prjkz(j,k)*rjkz(j,k)*ff3(i,j,k) c goto 555 c if(mod(istep,iprint).eq.0) then cc1tm2=0.0 cc2tm2=0.0 cc4tm2=0.0 cc1tm3=0.0 cc2tm3=0.0 cc4tm3=0.0 do 410 kk=1,neg if(kk.eq.j.or.kk.eq.k) goto 410 cc1tm2=cc1tm2+zet2(i,j,kk)*x(i,kk)*x(i,kk)/r(i,kk) cc2tm2=cc2tm2+zet2(i,j,kk)*y(i,kk)*y(i,kk)/r(i,kk) cc4tm2=cc4tm2+zet2(i,j,kk)*y(i,kk)*z(i,kk)/r(i,kk) cc1tm3=cc1tm3+zet3(i,j,kk)*rjkx(j,kk)*rjkx(j,kk)/rjk(j,kk) cc2tm3=cc2tm3+zet3(i,j,kk)*rjky(j,kk)*rjky(j,kk)/rjk(j,kk) cc4tm3=cc4tm3+zet3(i,j,kk)*rjky(j,kk)*rjkz(j,kk)/rjk(j,kk) 410 continue c23456789012345678901234567890123456789012345678901234567890123456789012 c11tm22(j)=c11tm22(j)+zet2(i,j,k)*x(i,k)*x(i,k)/r(i,k)*cc1tm2 c12tm22(j)=c12tm22(j)+zet2(i,j,k)*x(i,k)*x(i,k)/r(i,k)*cc2tm2 c44tm22(j)=c44tm22(j)+zet2(i,j,k)*y(i,k)*z(i,k)/r(i,k)*cc4tm2 c c11tm33(j)=c11tm33(j)+zet3(i,j,k)*rjkx(j,k)*rjkx(j,k) a /rjk(j,k)*cc1tm3 c12tm33(j)=c12tm33(j)+zet3(i,j,k)*rjkx(j,k)*rjkx(j,k) a /rjk(j,k)*cc2tm3 c44tm33(j)=c44tm33(j)+zet3(i,j,k)*rjky(j,k)*rjkz(j,k) a /rjk(j,k)*cc4tm3 c234567 c11tm23(j)=c11tm23(j)+zet2(i,j,k)*x(i,k)*x(i,k)/r(i,k)*cc1tm3 c12tm23(j)=c12tm23(j)+zet2(i,j,k)*x(i,k)*x(i,k)/r(i,k)*cc2tm3 c44tm23(j)=c44tm23(j)+zet2(i,j,k)*y(i,k)*z(i,k)/r(i,k)*cc4tm3 cc c11tm32(j)=c11tm32(j)+zet3(i,j,k)*rjkx(j,k)*rjkx(j,k) a /rjk(j,k)*cc1tm2 c12tm32(j)=c12tm32(j)+zet3(i,j,k)*rjkx(j,k)*rjkx(j,k) a /rjk(j,k)*cc2tm2 c44tm32(j)=c44tm32(j)+zet3(i,j,k)*rjky(j,k)*rjkz(j,k) a /rjk(j,k)*cc4tm2 c t=(ff22(i,j,k)-ff2(ipt,j,k)/r(i,k))/r2(i,k) c11tm2(j)=c11tm2(j)+t*x(i,k)*x(i,k)*x(i,k)*x(i,k) c12tm2(j)=c12tm2(j)+t*x(i,k)*x(i,k)*y(i,k)*y(i,k) c44tm2(j)=c44tm2(j)+t*y(i,k)*z(i,k)*y(i,k)*z(i,k) c t=(ff33(i,j,k)-ff3(ipt,j,k)/rjk(j,k))/rjk2(j,k) c11tm3(j)=c11tm3(j)+t*rjkx(j,k)*rjkx(j,k)*rjkx(j,k)*rjkx(j,k) c12tm3(j)=c12tm3(j)+t*rjkx(j,k)*rjkx(j,k)*rjky(j,k)*rjky(j,k) c44tm3(j)=c44tm3(j)+t*rjky(j,k)*rjkz(j,k)*rjky(j,k)*rjkz(j,k) c t=(x(i,j)*x(i,j)*x(i,k)*x(i,k)+x(i,k)*x(i,k)*x(i,j)*x(i,j)) c11tm4(j)=c11tm4(j)+ff12(ipt,j,k)*t/r(i,j)/r(i,k) t=(x(i,j)*x(i,j)*y(i,k)*y(i,k)+x(i,k)*x(i,k)*y(i,j)*y(i,j)) c12tm4(j)=c12tm4(j)+ff12(ipt,j,k)*t/r(i,j)/r(i,k) t=(y(i,j)*z(i,j)*y(i,k)*z(i,k)*2.0) c44tm4(j)=c44tm4(j)+ff12(ipt,j,k)*t/r(i,j)/r(i,k) c t=x(i,j)*x(i,j)*rjkx(j,k)*rjkx(j,k)*2.0 c11tm5(j)=c11tm5(j)+ff13(ipt,j,k)*t/rjk(j,k)/r(i,j) t=x(i,j)*x(i,j)*rjky(j,k)*rjky(j,k)+ a rjkx(j,k)*rjkx(j,k)*y(i,j)*y(i,j) c12tm5(j)=c12tm5(j)+ff13(ipt,j,k)*t/rjk(j,k)/r(i,j) t=y(i,j)*z(i,j)*rjky(j,k)*rjkz(j,k)*2.0 c44tm5(j)=c44tm5(j)+ff13(ipt,j,k)*t/rjk(j,k)/r(i,j) c c11tm6(j)=c11tm6(j)+ a ff23(i,j,k)*(x(i,k)*x(i,k)*rjkx(j,k)*rjkx(j,k) a +rjkx(j,k)*rjkx(j,k)*x(i,k)*x(i,k))/rjk(j,k)/r(i,k) c12tm6(j)=c12tm6(j)+ a ff23(i,j,k)*(x(i,k)*x(i,k)*rjky(j,k)*rjky(j,k) a +rjkx(j,k)*rjkx(j,k)*y(i,k)*y(i,k))/r(i,k)/rjk(j,k) c44tm6(j)=c44tm6(j)+ a ff23(i,j,k)*(y(i,k)*z(i,k)*rjky(j,k)*rjkz(j,k) a +rjky(j,k)*rjkz(j,k)*y(i,k)*z(i,k))/rjk(j,k)/r(i,k) c c endif 555 continue c c calc force on particle k c234567 fkx=-(ff2(i,j,k)*x(i,k)/r(i,k)+ff3(ipt,j,k)*(-prjkx(j,k)))/2.0 fky=-(ff2(i,j,k)*y(i,k)/r(i,k)+ff3(ipt,j,k)*(-prjky(j,k)))/2.0 fkz=-(ff2(i,j,k)*z(i,k)/r(i,k)+ff3(ipt,j,k)*(-prjkz(j,k)))/2.0 c fx(kp)=fx(kp) + fkx fy(kp)=fy(kp) + fky fz(kp)=fz(kp) + fkz 840 continue c c calc. force on particle i c fix=(frst(j)+pf(j)*szet1(j))*xt(j)+pf(j)*sftm2x(j) fx(ipt)=fx(ipt) + fix*0.50 fiy=(frst(j)+pf(j)*szet1(j))*yt(j)+pf(j)*sftm2y(j) fy(ipt)=fy(ipt) + fiy*0.50 fiz=(frst(j)+pf(j)*szet1(j))*zt(j)+pf(j)*sftm2z(j) fz(ipt)=fz(ipt) + fiz*0.50 c c calc force on particle j c fjx=-((frst(j)+pf(j)*szet1(j))*xt(j)+sftm3x(j)*pf(j))*0.50 fjy=-((frst(j)+pf(j)*szet1(j))*yt(j)+sftm3y(j)*pf(j))*0.50 fjz=-((frst(j)+pf(j)*szet1(j))*zt(j)+sftm3z(j)*pf(j))*0.50 c fx(jpt)=fx(jpt) + fjx fy(jpt)=fy(jpt) + fjy fz(jpt)=fz(jpt) + fjz c c calc. the stress tensor of the system. c pxx=pxx-(ff1(ipt,j)*xt(j)*x(i,j)+strs21(j)+strs31(j))/2.0 pxy=pxy-(ff1(ipt,j)*xt(j)*y(i,j)+strs22(j)+strs32(j))/2.0 pxz=pxz-(ff1(ipt,j)*xt(j)*z(i,j)+strs23(j)+strs33(j))/2.0 pyy=pyy-(ff1(ipt,j)*yt(j)*y(i,j)+strs24(j)+strs34(j))/2.0 pyz=pyz-(ff1(ipt,j)*yt(j)*z(i,j)+strs25(j)+strs35(j))/2.0 pzz=pzz-(ff1(ipt,j)*zt(j)*z(i,j)+strs26(j)+strs36(j))/2.0 702 continue do 703 j=1,neg jpt=list2(j) n22=ntyp(jpt) c goto 703 c if(mod(istep,iprint).eq.0) then tm1=ff11(ipt,j)-ff1(ipt,j)/r(i,j) c ttm=x(i,j)*x(i,j)*x(i,j)*x(i,j)/r2(i,j) c11p=c11p+tm1*ttm+pf3(i,j)*(c11tm22(j)+c11tm33(j) a +c11tm23(j)+c11tm32(j)) a +(c11tm2(j)+c11tm3(j))+c11tm4(j)+c11tm5(j)+c11tm6(j) c if(n1.eq.1.and.n22.eq.1) then cpb(1,1)=cpb(1,1)+tm1*ttm+pf3(i,j)*(c11tm22(j)+c11tm33(j) a +c11tm23(j)+c11tm32(j)) a +(c11tm2(j)+c11tm3(j))+c11tm4(j)+c11tm5(j)+c11tm6(j) endif if(n1.eq.1.and.n22.eq.2) then cpb(1,2)=cpb(1,2)+tm1*ttm+pf3(i,j)*(c11tm22(j)+c11tm33(j) a +c11tm23(j)+c11tm32(j)) a +(c11tm2(j)+c11tm3(j))+c11tm4(j)+c11tm5(j)+c11tm6(j) endif if(n1.eq.2.and.n22.eq.1) then cpb(1,3)=cpb(1,3)+tm1*ttm+pf3(i,j)*(c11tm22(j)+c11tm33(j) a +c11tm23(j)+c11tm32(j)) a +(c11tm2(j)+c11tm3(j))+c11tm4(j)+c11tm5(j)+c11tm6(j) endif if(n1.eq.2.and.n22.eq.2) then cpb(1,4)=cpb(1,4)+tm1*ttm+pf3(i,j)*(c11tm22(j)+c11tm33(j) a + c11tm23(j)+c11tm32(j)) a +(c11tm2(j)+c11tm3(j))+c11tm4(j)+c11tm5(j)+c11tm6(j) endif c ttm=x(i,j)*x(i,j)*y(i,j)*y(i,j)/r2(i,j) c12p=c12p+tm1*ttm+pf3(i,j)*(c12tm22(j)+c12tm33(j) a +c12tm23(j)+c12tm32(j)) a +(c12tm2(j)+c12tm3(j))+c12tm4(j)+c12tm5(j)+c12tm6(j) c if(n1.eq.1.and.n22.eq.1) then cpb(2,1)=cpb(2,1)+tm1*ttm+pf3(i,j)*(c12tm22(j)+c12tm33(j) a +c12tm23(j)+c12tm32(j)) a +(c12tm2(j)+c12tm3(j))+c12tm4(j)+c12tm5(j)+c12tm6(j) endif if(n1.eq.1.and.n22.eq.2) then cpb(2,2)=cpb(2,2)+tm1*ttm+pf3(i,j)*(c12tm22(j)+c12tm33(j) a +c12tm23(j)+c12tm32(j)) a +(c12tm2(j)+c12tm3(j))+c12tm4(j)+c12tm5(j)+c12tm6(j) endif if(n1.eq.2.and.n22.eq.1) then cpb(2,3)=cpb(2,3)+tm1*ttm+pf3(i,j)*(c12tm22(j)+c12tm33(j) a +c12tm23(j)+c12tm32(j)) a +(c12tm2(j)+c12tm3(j))+c12tm4(j)+c12tm5(j)+c12tm6(j) endif if(n1.eq.2.and.n22.eq.2) then cpb(2,4)=cpb(2,4)+tm1*ttm+pf3(i,j)*(c12tm22(j)+c12tm33(j) a +c12tm23(j)+c12tm32(j)) a +(c12tm2(j)+c12tm3(j))+c12tm4(j)+c12tm5(j)+c12tm6(j) endif c ttm=y(i,j)*z(i,j)*y(i,j)*z(i,j)/r2(i,j) c44p=c44p+tm1*ttm+pf3(i,j)*(c44tm22(j)+c44tm33(j) a +c44tm23(j)+c44tm32(j)) a +(c44tm2(j)+c44tm3(j))+c44tm4(j)+c44tm5(j)+c44tm6(j) c if(n1.eq.1.and.n22.eq.1) then cpb(3,1)=cpb(3,1)+tm1*ttm+pf3(i,j)*(c44tm22(j)+c44tm33(j) a +c44tm23(j)+c44tm32(j)) a +(c44tm2(j)+c44tm3(j))+c44tm4(j)+c44tm5(j)+c44tm6(j) endif if(n1.eq.1.and.n22.eq.2) then cpb(3,2)=cpb(3,2)+tm1*ttm+pf3(i,j)*(c44tm22(j)+c44tm33(j) a +c44tm23(j)+c44tm32(j)) a +(c44tm2(j)+c44tm3(j))+c44tm4(j)+c44tm5(j)+c44tm6(j) endif if(n1.eq.2.and.n22.eq.1) then cpb(3,3)=cpb(3,3)+tm1*ttm+pf3(i,j)*(c44tm22(j)+c44tm33(j) a +c44tm23(j)+c44tm32(j)) a +(c44tm2(j)+c44tm3(j))+c44tm4(j)+c44tm5(j)+c44tm6(j) endif if(n1.eq.2.and.n22.eq.2) then cpb(3,4)=cpb(3,4)+tm1*ttm+pf3(i,j)*(c44tm22(j)+c44tm33(j) a +c44tm23(j)+c44tm32(j)) a +(c44tm2(j)+c44tm3(j))+c44tm4(j)+c44tm5(j)+c44tm6(j) endif c endif c 703 continue c 500 continue C C C THIS SUBROUTINE WILL CALCULATE THE DYNAMIC MATRIX C THE NOTATION FOLLOWS THE ANALYTIC DERIVATION DONE BY C L.J.PORTER. IE, LLOOP AND LLOOP2 REFER TO THE FIRST SUM OVER L C AND THE SECOND SUM OVER L FOUND IN THOSE NOTES. C C THE DERIVATIVE NOTATION IS A "HYBRID EINSTEIN NOTATION". C FOR EXAMPLE, DIJDIJVIL = THE SECOND DERIVATIVE OF VIL WITH C RESPECT TO RIJ. DIJDIKVIJ = THE SECOND DERIVATIVE OF VIJ WITH C RESPECT TO RIJ AND RIK. ETC C C put the open statement in the beginning of the code C c open(unit=28,file='dynamorph.out') open(unit=29,file='rigrot.out') open(unit=91,file='omegas.out') C234567 c IF(MOD(ISTEP,IPRINT).EQ.0) THEN C c write(28,*) ' the step # = ',istep c write(28,*) ' ' C C zero the dynamic matrix do 17 i = 1,np do 18 j = 1,np DXX(i,j) = 0.0 DXY(i,j) = 0.0 DXZ(i,j) = 0.0 DYY(i,j) = 0.0 DZZ(i,j) = 0.0 DYX(i,j) = 0.0 DYZ(i,j) = 0.0 DZX(i,j) = 0.0 DZY(i,j) = 0.0 18 continue 17 continue C C C MEIJIE DEFINED XIJ,YIJ,ZIJ AS XJ - XI, ETC C NOTE ALSO THAT BELOW, I HAVE QUANTITIES LIKE C RJNX = -1*(X(I,J) - X(I,N)) BECAUSE AFTER THE C LOOPS BELOW, X(I,J) WILL = XI - XJ. THUS, C X(I,J) - X(I,N) WOULD EQUAL XN - XJ = -RJNX. C c do 21 ii= 1,np do 22 jj = 1,nn(ii) x(ii,jj) = -1.0*x(ii,jj) y(ii,jj) = -1.0*y(ii,jj) z(ii,jj) = -1.0*z(ii,jj) 22 continue 21 continue c c C IN THIS 1002 LOOP, WE EVENTUALLY WANT TO LOOP OVER ALL PARTICLES, C SO DO 1002 I = 1, NP SHOULD EVENTUALLY BE USED C do 1002 i = 1, np c c write(6,*) 'inside D loop: i = ',i c do 1100 j = 1, np c if (j .eq. i) goto 1100 c c zero everything firstxx = 0.0 firstyy = 0.0 firstzz = 0.0 firstxy = 0.0 firstxz = 0.0 firstyx = 0.0 firstyz = 0.0 firstzx = 0.0 firstzy = 0.0 c lloopxx = 0.0 lloopyy = 0.0 lloopzz = 0.0 lloopxy = 0.0 lloopxz = 0.0 lloopyx = 0.0 lloopyz = 0.0 lloopzx = 0.0 lloopzy = 0.0 c lloop2xx = 0.0 lloop2yy = 0.0 lloop2zz = 0.0 lloop2xy = 0.0 lloop2xz = 0.0 lloop2yx = 0.0 lloop2yz = 0.0 lloop2zx = 0.0 lloop2zy = 0.0 c lloop5xx = 0.0 lloop5yy = 0.0 lloop5zz = 0.0 lloop5xy = 0.0 lloop5xz = 0.0 lloop5yx = 0.0 lloop5yz = 0.0 lloop5zx = 0.0 lloop5zy = 0.0 c jid = j negji = negno(i,j) negi = negno(jid,i) if (negji.eq.0) goto 1800 c rij = r(i,negji) rijx = x(i,negji) rijy = y(i,negji) rijz = z(i,negji) rij2 = r2(i,negji) rij3 = rij2*rij C C SIMPLE CHECK C rji = r(jid,negi) if(rji .ne. rij) then write(6,*) 'something is wrong ' stop c else c write(6,*) 'rij = rji checks' endif C C dijvij = ff1(i,negji) dijdijvij = ff11(i,negji) dijvji = ff1(jid,negi) dijdijvji = ff11(jid,negi) c23456 firstxx = (rijx*rijx/rij3 - 1.0/rij)*(dijvij + dijvji) + a (rijx*(-rijx)/rij2)*dijdijvij + (rijx*(-rijx)/rij2)*dijdijvji c firstyy = (rijy*rijy/rij3 - 1.0/rij)*(dijvij + dijvji) + a (rijy*(-rijy)/rij2)*dijdijvij + (rijy*(-rijy)/rij2)*dijdijvji c firstxy = (rijx*rijy/rij3)*(dijvij + dijvji) + a (rijx*(-rijy)/rij2)*dijdijvij + (rijx*(-rijy)/rij2)*dijdijvji c firstyx = (rijy*rijx/rij3)*(dijvij + dijvji) + a (rijy*(-rijx)/rij2)*dijdijvij + (rijy*(-rijx)/rij2)*dijdijvji c firstzz = (rijz*rijz/rij3 - 1.0/rij)*(dijvij + dijvji) + a (rijz*(-rijz)/rij2)*dijdijvij + (rijz*(-rijz)/rij2)*dijdijvji c firstxz = (rijx*rijz/rij3)*(dijvij + dijvji) + a (rijx*(-rijz)/rij2)*dijdijvij + (rijx*(-rijz)/rij2)*dijdijvji c firstyz = (rijy*rijz/rij3)*(dijvij + dijvji) + a (rijy*(-rijz)/rij2)*dijdijvij + (rijy*(-rijz)/rij2)*dijdijvji c firstzx = (rijz*rijx/rij3)*(dijvij + dijvji) + a (rijz*(-rijx)/rij2)*dijdijvij + (rijz*(-rijx)/rij2)*dijdijvji c firstzy = (rijz*rijy/rij3)*(dijvij + dijvji) + a (rijz*(-rijy)/rij2)*dijdijvij + (rijz*(-rijy)/rij2)*dijdijvji c c c do 1150 l = 1, nn(i) if (l .eq. negji) goto 1150 ril = r(i,l) ril2 = ril*ril rilx = x(i,l) rily = y(i,l) rilz = z(i,l) rmujl = (x(i,negji)*x(i,l) + y(i,negji)*y(i,l) + a z(i,negji)*z(i,l))/(rij*ril) rjl2 = rij2 + ril2 -2.0*rmujl*rij*ril rjlx = -1.0*(x(i,negji) - x(i,l)) rjly = -1.0*(y(i,negji) - y(i,l)) rjlz = -1.0*(z(i,negji) - z(i,l)) rjl = sqrt(rjl2) c nloopxx = 0.0 nloopxy = 0.0 nloopyx = 0.0 nloopyy = 0.0 nloopzz = 0.0 nloopxz = 0.0 nloopyz = 0.0 nloopzx = 0.0 nloopzy = 0.0 c do 1200 m = 1,nn(i) if ((m .eq. negji).or.(m.eq.l)) goto 1200 rinx = x(i,m) riny = y(i,m) rinz = z(i,m) rin = r(i,m) c dindjlvij = pf3(i,negji)*zet2(i,negji,m)*zet3(i,negji,l) dildijvin = pf3(i,m)*zet2(i,m,negji)*zet2(i,m,l) dindjlvil = pf3(i,l)*zet2(i,l,m)*zet3(i,l,negji) c nloopxx = nloopxx + (rilx*(-rijx)/(rij*ril))*dildijvin + a (rinx*rjlx/(rin*rjl))*(dindjlvij + dindjlvil) c nloopyy = nloopyy + (rily*(-rijy)/(rij*ril))*dildijvin + a (riny*rjly/(rin*rjl))*(dindjlvij + dindjlvil) c nloopxy = nloopxy + (rilx*(-rijy)/(rij*ril))*dildijvin + a (rinx*rjly/(rin*rjl))*(dindjlvij + dindjlvil) c nloopyx = nloopyx + (rily*(-rijx)/(rij*ril))*dildijvin + a (riny*rjlx/(rin*rjl))*(dindjlvij + dindjlvil) c nloopzz = nloopzz + (rilz*(-rijz)/(rij*ril))*dildijvin + a (rinz*rjlz/(rin*rjl))*(dindjlvij + dindjlvil) c nloopxz = nloopxz + (rilx*(-rijz)/(rij*ril))*dildijvin + a (rinx*rjlz/(rin*rjl))*(dindjlvij + dindjlvil) c nloopyz = nloopyz + (rily*(-rijz)/(rij*ril))*dildijvin + a (riny*rjlz/(rin*rjl))*(dindjlvij + dindjlvil) c nloopzx = nloopzx + (rilz*(-rijx)/(rij*ril))*dildijvin + a (rinz*rjlx/(rin*rjl))*(dindjlvij + dindjlvil) c nloopzy = nloopzy + (rilz*(-rijy)/(rij*ril))*dildijvin + a (rinz*rjly/(rin*rjl))*(dindjlvij + dindjlvil) c 1200 continue c c dijvil = ff2(i,l,negji) dijdijvil = ff22(i,l,negji) dijdjlvij = ff13(i,negji,l) dijdjlvil = ff23(i,l,negji) dijdilvil = ff12(i,l,negji) dijdilvij = ff12(i,negji,l) djldilvij = ff23(i,negji,l) djldilvil = ff13(i,l,negji) c23456 lloopxx = lloopxx + (rijx*rijx/rij3 - 1.0/rij)*dijvil + a (rijx*(-rijx)/rij2)*dijdijvil + a (rijx*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rilx*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) + a (rilx*rjlx/(ril*rjl))*(djldilvil + djldilvij) + a nloopxx c lloopyy = lloopyy + (rijy*rijy/rij3 - 1.0/rij)*dijvil + a (rijy*(-rijy)/rij2)*dijdijvil + a (rijy*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rily*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) + a (rily*rjly/(ril*rjl))*(djldilvil + djldilvij) + a nloopyy c lloopxy = lloopxy + (rijx*rijy/rij3)*dijvil + a (rijx*(-rijy)/rij2)*dijdijvil + a (rijx*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rilx*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) + a (rilx*rjly/(ril*rjl))*(djldilvil + djldilvij) + a nloopxy c lloopyx = lloopyx + (rijy*rijx/rij3)*dijvil + a (rijy*(-rijx)/rij2)*dijdijvil + a (rijy*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rily*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) + a (rily*rjlx/(ril*rjl))*(djldilvil + djldilvij) + a nloopyx c lloopzz = lloopzz + (rijz*rijz/rij3 - 1.0/rij)*dijvil + a (rijz*(-rijz)/rij2)*dijdijvil + a (rijz*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rilz*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) + a (rilz*rjlz/(ril*rjl))*(djldilvil + djldilvij) + a nloopzz c lloopxz = lloopxz + (rijx*rijz/rij3)*dijvil + a (rijx*(-rijz)/rij2)*dijdijvil + a (rijx*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rilx*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) + a (rilx*rjlz/(ril*rjl))*(djldilvil + djldilvij) + a nloopxz c lloopyz = lloopyz + (rijy*rijz/rij3)*dijvil + a (rijy*(-rijz)/rij2)*dijdijvil + a (rijy*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rily*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) + a (rily*rjlz/(ril*rjl))*(djldilvil + djldilvij) + a nloopyz c lloopzx = lloopzx + (rijz*rijx/rij3)*dijvil + a (rijz*(-rijx)/rij2)*dijdijvil + a (rijz*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rilz*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) + a (rilz*rjlx/(ril*rjl))*(djldilvil + djldilvij) + a nloopzx c lloopzy = lloopzy + (rijz*rijy/rij3)*dijvil + a (rijz*(-rijy)/rij2)*dijdijvil + a (rijz*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) + a (rilz*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) + a (rilz*rjly/(ril*rjl))*(djldilvil + djldilvij) + a nloopzy c 1150 continue c c c do 1250 l = 1,nn(jid) lid = idpart(jid,l) if(lid .eq. i) goto 1250 rjl = r(jid,l) rjl2 = rjl*rjl rjlx = x(jid,l) rjly = y(jid,l) rjlz = z(jid,l) rmuil = (x(jid,negi)*x(jid,l) + y(jid,negi)*y(jid,l) + a z(jid,negi)*z(jid,l))/(rij*rjl) ril2 = rij2 + rjl2 - 2.0*rmuil*rij*rjl ril = sqrt(ril2) rilx = -1.0*(x(jid,negi) - x(jid,l)) rily = -1.0*(y(jid,negi) - y(jid,l)) rilz = -1.0*(z(jid,negi) - z(jid,l)) c nloop2xx = 0.0 nloop2xy = 0.0 nloop2yx = 0.0 nloop2yy = 0.0 nloop2zz = 0.0 nloop2xz = 0.0 nloop2yz = 0.0 nloop2zx = 0.0 nloop2zy = 0.0 c do 1300 m = 1, nn(jid) nid = idpart(jid,m) if ((nid.eq.i) .or.(m.eq.l)) goto 1300 rjn = r(jid,m) rjn2 = rjn*rjn rmuin = (x(jid,negi)*x(jid,m) + y(jid,negi)*y(jid,m) + a z(jid,negi)*z(jid,m))/(rij*rjn) rin2 = rij2 + rjn2 - 2.0*rmuin*rij*rjn rin = sqrt(rin2) rinx = -1.0*(x(jid,negi) - x(jid,m)) riny = -1.0*(y(jid,negi) - y(jid,m)) rinz = -1.0*(z(jid,negi) - z(jid,m)) c djldijvjn = pf3(jid,m)*zet2(jid,m,negi)*zet2(jid,m,l) djldinvji = pf3(jid,negi)*zet2(jid,negi,l)*zet3(jid,negi,m) djldinvjn = pf3(jid,m)*zet2(jid,m,l)*zet3(jid,m,negi) c nloop2xx = nloop2xx + (rijx*rjlx/(rij*rjl))*djldijvjn a +(rinx*rjlx/(rin*rjl))*(djldinvji + djldinvjn) c nloop2yy = nloop2yy + (rijy*rjly/(rij*rjl))*djldijvjn a +(riny*rjly/(rin*rjl))*(djldinvji + djldinvjn) c nloop2xy = nloop2xy + (rijx*rjly/(rij*rjl))*djldijvjn a +(rinx*rjly/(rin*rjl))*(djldinvji + djldinvjn) c nloop2yx = nloop2yx + (rijy*rjlx/(rij*rjl))*djldijvjn a +(riny*rjlx/(rin*rjl))*(djldinvji + djldinvjn) c nloop2zz = nloop2zz + (rijz*rjlz/(rij*rjl))*djldijvjn a +(rinz*rjlz/(rin*rjl))*(djldinvji + djldinvjn) c nloop2xz = nloop2xz + (rijx*rjlz/(rij*rjl))*djldijvjn a +(rinx*rjlz/(rin*rjl))*(djldinvji + djldinvjn) c nloop2yz = nloop2yz + (rijy*rjlz/(rij*rjl))*djldijvjn a +(riny*rjlz/(rin*rjl))*(djldinvji + djldinvjn) c nloop2zx = nloop2zx + (rijz*rjlx/(rij*rjl))*djldijvjn a +(rinz*rjlx/(rin*rjl))*(djldinvji + djldinvjn) c nloop2zy = nloop2zy + (rijz*rjly/(rij*rjl))*djldijvjn a +(rinz*rjly/(rin*rjl))*(djldinvji + djldinvjn) c 1300 continue c dijvjl = ff2(jid,l,negi) djldijvji = ff12(jid,negi,l) dijdijvjl = ff22(jid,l,negi) dildjlvjl = ff13(jid,l,negi) dijdjlvjl = ff12(jid,l,negi) dildjlvji = ff23(jid,negi,l) dijdilvji = ff13(jid,negi,l) dijdilvjl = ff23(jid,l,negi) c lloop2xx = lloop2xx +(rijx*rijx/rij3 - 1.0/rij)*dijvjl + a (rijx*rjlx/(rij*rjl))*djldijvji + a (rijx*(-rijx)/rij2)*dijdijvjl + a (rilx*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijx*rjlx/(rij*rjl))*dijdjlvjl + a (rilx*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2xx c lloop2yy = lloop2yy +(rijy*rijy/rij3 - 1.0/rij)*dijvjl + a (rijy*rjly/(rij*rjl))*djldijvji + a (rijy*(-rijy)/rij2)*dijdijvjl + a (rily*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijy*rjly/(rij*rjl))*dijdjlvjl + a (rily*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2yy c lloop2xy = lloop2xy +(rijx*rijy/rij3)*dijvjl + a (rijx*rjly/(rij*rjl))*djldijvji + a (rijx*(-rijy)/rij2)*dijdijvjl + a (rilx*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijx*rjly/(rij*rjl))*dijdjlvjl + a (rilx*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2xy c lloop2yx = lloop2yx +(rijy*rijx/rij3)*dijvjl + a (rijy*rjlx/(rij*rjl))*djldijvji + a (rijy*(-rijx)/rij2)*dijdijvjl + a (rily*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijy*rjlx/(rij*rjl))*dijdjlvjl + a (rily*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2yx c lloop2zz = lloop2zz +(rijz*rijz/rij3 - 1.0/rij)*dijvjl + a (rijz*rjlz/(rij*rjl))*djldijvji + a (rijz*(-rijz)/rij2)*dijdijvjl + a (rilz*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijz*rjlz/(rij*rjl))*dijdjlvjl + a (rilz*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2zz c lloop2xz = lloop2xz +(rijx*rijz/rij3)*dijvjl + a (rijx*rjlz/(rij*rjl))*djldijvji + a (rijx*(-rijz)/rij2)*dijdijvjl + a (rilx*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijx*rjlz/(rij*rjl))*dijdjlvjl + a (rilx*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2xz c lloop2yz = lloop2yz +(rijy*rijz/rij3)*dijvjl + a (rijy*rjlz/(rij*rjl))*djldijvji + a (rijy*(-rijz)/rij2)*dijdijvjl + a (rily*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijy*rjlz/(rij*rjl))*dijdjlvjl + a (rily*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2yz c lloop2zx = lloop2zx +(rijz*rijx/rij3)*dijvjl + a (rijz*rjlx/(rij*rjl))*djldijvji + a (rijz*(-rijx)/rij2)*dijdijvjl + a (rilz*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijz*rjlx/(rij*rjl))*dijdjlvjl + a (rilz*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2zx c lloop2zy = lloop2zy +(rijz*rijy/rij3)*dijvjl + a (rijz*rjly/(rij*rjl))*djldijvji + a (rijz*(-rijy)/rij2)*dijdijvjl + a (rilz*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+ a (rijz*rjly/(rij*rjl))*dijdjlvjl + a (rilz*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+ a nloop2zy c 1250 continue c c C THIS PART WILL GIVE NONZERO CONTRIBUTIONS WHEN J = NN(I) IF C THE STRUCTURE IS AMORPHOUS. IT WILL GIVE NONZERO CONTRIBUTIONS C IF J = NNN(I) FOR BOTH TETRAHEDRAL AND AMORPHOUS STRUCTURES. C IN THE ORIGINAL CALCULATION, IT REPRESENTS THE CALCULATION OF C D(I,K), WHERE K = NNN(I) C 1800 continue c23456789012345678901234567890123456789012345678901234567890123456789012 c do 1700 l = 1,nn(i) lid = idpart(i,l) if (lid .eq. jid) goto 1700 if (negno(jid,lid) .eq. 0) goto 1700 negil = negno(lid,i) negjl = negno(lid,jid) c ril = r(i,l) rilx = x(i,l) rily = y(i,l) rilz = z(i,l) ril2 = ril*ril rlj = r(lid,negjl) rlj2 = rlj*rlj rljx = x(lid,negjl) rljy = y(lid,negjl) rljz = z(lid,negjl) c c23456789012345678901234567890123456789012345678901234567890123456789012 c The following is a calculation of the distance rij if j is not c a nearest neighbor of i. From above, l must be a neighbor of both c i and j, so we can get rij through l. c if (negji .eq. 0) then rmuij = (x(lid,negil)*x(lid,negjl) + y(lid,negil)*y(lid,negjl) a +z(lid,negil)*z(lid,negjl))/(ril*rlj) rij2 = ril2 + rlj2 -2.0*rmuij*ril*rlj rij = sqrt(rij2) rij3 = rij2*rij rijx = -1.0*(x(lid,negil) - x(lid,negjl)) rijy = -1.0*(y(lid,negil) - y(lid,negjl)) rijz = -1.0*(z(lid,negil) - z(lid,negjl)) endif c nloop3xx = 0.0 nloop3xy = 0.0 nloop3xz = 0.0 nloop3yx = 0.0 nloop3yy = 0.0 nloop3yz = 0.0 nloop3zx = 0.0 nloop3zy = 0.0 nloop3zz = 0.0 c do 1750 m = 1,nn(lid) nid = idpart(lid,m) if((nid.eq.i).or.(nid.eq.jid)) goto 1750 rln = r(lid,m) rln2 = rln*rln rmuin = (x(lid,negil)*x(lid,m) + y(lid,negil)*y(lid,m) + a z(lid,negil)*z(lid,m))/(ril*rln) rin2 = ril2 + rln2 - 2.0*rmuin*ril*rln rin = sqrt(rin2) rinx = -1.0*(x(lid,negil) - x(lid,m)) riny = -1.0*(y(lid,negil) - y(lid,m)) rinz = -1.0*(z(lid,negil) - z(lid,m)) rmujn = (x(lid,negjl)*x(lid,m) + y(lid,negjl)*y(lid,m) + a z(lid,negjl)*z(lid,m))/(rlj*rln) rjn2 = rln2 + rlj2 - 2.0*rmujn*rln*rlj rjn = sqrt(rjn2) rjnx = -1.0*(x(lid,negjl) - x(lid,m)) rjny = -1.0*(y(lid,negjl) - y(lid,m)) rjnz = -1.0*(z(lid,negjl) - z(lid,m)) c c23456789012345678901234567890123456789012345678901234567890123456789012 c dindjlvli = pf3(lid,negil)*zet2(lid,negil,negjl) a *zet3(lid,negil,m) dindjlvln = pf3(lid,m)*zet2(lid,m,negjl)*zet3(lid,m,negil) dildjnvlj= pf3(lid,negjl)*zet2(lid,negjl,negil) a *zet3(lid,negjl,m) dildjnvln = pf3(lid,m)*zet2(lid,m,negil)*zet3(lid,m,negjl) dindijvli = pf3(lid,negil)*zet3(lid,negil,m) a *zet3(lid,negil,negjl) dildjlvln = pf3(lid,m)*zet2(lid,m,negil)*zet2(lid,m,negjl) dijdjnvlj = pf3(lid,negjl)*zet3(lid,negjl,negil) a *zet3(lid,negjl,m) dindjnvln = pf3(lid,m)*zet3(lid,m,negil)*zet3(lid,m,negjl) c nloop3xx = (rinx*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rilx*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) + a (rinx*(-rijx)/(rin*rij))*dindijvli + a (rilx*(-rljx)/(ril*rlj))*dildjlvln + a (rijx*rjnx/(rij*rjn))*dijdjnvlj + a (rinx*rjnx/(rin*rjn))*dindjnvln + nloop3xx c nloop3yy = (riny*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rily*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) + a (riny*(-rijy)/(rin*rij))*dindijvli + a (rily*(-rljy)/(ril*rlj))*dildjlvln + a (rijy*rjny/(rij*rjn))*dijdjnvlj + a (riny*rjny/(rin*rjn))*dindjnvln + nloop3yy c nloop3zz = (rinz*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rilz*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) + a (rinz*(-rijz)/(rin*rij))*dindijvli + a (rilz*(-rljz)/(ril*rlj))*dildjlvln + a (rijz*rjnz/(rij*rjn))*dijdjnvlj + a (rinz*rjnz/(rin*rjn))*dindjnvln + nloop3zz c nloop3xy = (rinx*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rilx*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) + a (rinx*(-rijy)/(rin*rij))*dindijvli + a (rilx*(-rljy)/(ril*rlj))*dildjlvln + a (rijx*rjny/(rij*rjn))*dijdjnvlj + a (rinx*rjny/(rin*rjn))*dindjnvln + nloop3xy c nloop3xz = (rinx*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rilx*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) + a (rinx*(-rijz)/(rin*rij))*dindijvli + a (rilx*(-rljz)/(ril*rlj))*dildjlvln + a (rijx*rjnz/(rij*rjn))*dijdjnvlj + a (rinx*rjnz/(rin*rjn))*dindjnvln + nloop3xz c nloop3yx = (riny*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rily*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) + a (riny*(-rijx)/(rin*rij))*dindijvli + a (rily*(-rljx)/(ril*rlj))*dildjlvln + a (rijy*rjnx/(rij*rjn))*dijdjnvlj + a (riny*rjnx/(rin*rjn))*dindjnvln + nloop3yx c nloop3yz = (riny*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rily*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) + a (riny*(-rijz)/(rin*rij))*dindijvli + a (rily*(-rljz)/(ril*rlj))*dildjlvln + a (rijy*rjnz/(rij*rjn))*dijdjnvlj + a (riny*rjnz/(rin*rjn))*dindjnvln + nloop3yz c nloop3zx = (rinz*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rilz*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) + a (rinz*(-rijx)/(rin*rij))*dindijvli + a (rilz*(-rljx)/(ril*rlj))*dildjlvln + a (rijz*rjnx/(rij*rjn))*dijdjnvlj + a (rinz*rjnx/(rin*rjn))*dindjnvln + nloop3zx c nloop3zy = (rinz*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) + a (rilz*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) + a (rinz*(-rijy)/(rin*rij))*dindijvli + a (rilz*(-rljy)/(ril*rlj))*dildjlvln + a (rijz*rjny/(rij*rjn))*dijdjnvlj + a (rinz*rjny/(rin*rjn))*dindjnvln + nloop3zy c 1750 continue c c dijvli = ff3(lid,negil,negjl) dijvlj = ff3(lid,negjl,negil) dildjlvli = ff12(lid,negil,negjl) dildjlvlj = ff12(lid,negjl,negil) dildijvli = ff13(lid,negil,negjl) dildijvlj = ff23(lid,negjl,negil) dijdjlvli = ff23(lid,negil,negjl) dijdjlvlj = ff13(lid,negjl,negil) dijdijvli = ff33(lid,negil,negjl) dijdijvlj = ff33(lid,negjl,negil) c c23456789012345678901234567890123456789012345678901234567890123456789012 c lloop5xx = (rijx*rijx/rij3 - 1.0/rij)*(dijvli + dijvlj) + a (rilx*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rilx*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) + a (rijx*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijx*(-rijx)/rij2)*(dijdijvli + dijdijvlj) + a nloop3xx + lloop5xx c lloop5yy = (rijy*rijy/rij3 - 1.0/rij)*(dijvli + dijvlj) + a (rily*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rily*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) + a (rijy*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijy*(-rijy)/rij2)*(dijdijvli + dijdijvlj) + a nloop3yy + lloop5yy c lloop5zz = (rijz*rijz/rij3 - 1.0/rij)*(dijvli + dijvlj) + a (rilz*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rilz*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) + a (rijz*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijz*(-rijz)/rij2)*(dijdijvli + dijdijvlj) + a nloop3zz + lloop5zz c lloop5xy = (rijx*rijy/rij3)*(dijvli + dijvlj) + a (rilx*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rilx*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) + a (rijx*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijx*(-rijy)/rij2)*(dijdijvli + dijdijvlj) + a nloop3xy + lloop5xy c lloop5xz = (rijx*rijz/rij3)*(dijvli + dijvlj) + a (rilx*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rilx*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) + a (rijx*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijx*(-rijz)/rij2)*(dijdijvli + dijdijvlj) + a nloop3xz + lloop5xz c lloop5yx = (rijy*rijx/rij3)*(dijvli + dijvlj) + a (rily*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rily*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) + a (rijy*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijy*(-rijx)/rij2)*(dijdijvli + dijdijvlj) + a nloop3yx + lloop5yx c lloop5yz = (rijy*rijz/rij3)*(dijvli + dijvlj) + a (rily*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rily*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) + a (rijy*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijy*(-rijz)/rij2)*(dijdijvli + dijdijvlj) + a nloop3yz + lloop5yz c lloop5zx = (rijz*rijx/rij3)*(dijvli + dijvlj) + a (rilz*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rilz*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) + a (rijz*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijz*(-rijx)/rij2)*(dijdijvli + dijdijvlj) + a nloop3zx + lloop5zx c lloop5zy = (rijz*rijy/rij3)*(dijvli + dijvlj) + a (rilz*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) + a (rilz*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) + a (rijz*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) + a (rijz*(-rijy)/rij2)*(dijdijvli + dijdijvlj) + a nloop3zy + lloop5zy c 1700 continue c c c23456789012345678901234567890123456789012345678901234567890123456789012 DXX(i,jid) = 0.50*(firstxx + lloopxx + lloop2xx + lloop5xx) DXY(i,jid) = 0.50*(firstxy + lloopxy + lloop2xy + lloop5xy) DYX(i,jid) = 0.50*(firstyx + lloopyx + lloop2yx + lloop5yx) DYY(i,jid) = 0.50*(firstyy + lloopyy + lloop2yy + lloop5yy) DZZ(i,jid) = 0.50*(firstzz + lloopzz + lloop2zz + lloop5zz) DXZ(i,jid) = 0.50*(firstxz + lloopxz + lloop2xz + lloop5xz) DYZ(i,jid) = 0.50*(firstyz + lloopyz + lloop2yz + lloop5yz) DZX(i,jid) = 0.50*(firstzx + lloopzx + lloop2zx + lloop5zx) DZY(i,jid) = 0.50*(firstzy + lloopzy + lloop2zy + lloop5zy) c23456 c for purposes of checking, only write if D is nonzero c if (DXX(i,jid) .ne. 0) then c write(28,*) i,jid c write(28,*) ' DXX = ',DXX(i,jid),'DYY = ',DYY(i,jid), c a 'DZZ = ',DZZ(i,jid) c write(28,*)'DXY = ',DXY(i,jid),'DXZ = ',DXZ(i,jid) c write(28,*)'DYX = ',DYX(i,jid),'DYZ = ',DYZ(i,jid) c write(28,*)'DZX = ',DZX(i,jid),'DZY = ',DZY(i,jid) c write(28,*) ' ' c endif c 1100 continue c C THIS PART CALCULATES DXX(I,I) DXY(I,I) ETC, IE, J = I C sumkxx = 0.0 sumkyy = 0.0 sumkzz = 0.0 sumkxy = 0.0 sumkxz = 0.0 sumkyx = 0.0 sumkyz = 0.0 sumkzx = 0.0 sumkzy = 0.0 c do 1500 k = 1,nn(i) kid = idpart(i,k) negi = negno(kid,i) rik = r(i,k) rik2 = rik*rik rik3 = rik2*rik rikx = x(i,k) riky = y(i,k) rikz = z(i,k) dikvik = ff1(i,k) dikvki = ff1(kid,negi) dikdikvik = ff11(i,k) dikdikvki = ff11(kid,negi) c firstxx = (1./rik - rikx*rikx/rik3)*(dikvik + dikvki) + a (rikx*rikx/rik2)*(dikdikvik + dikdikvki) c firstyy = (1./rik - riky*riky/rik3)*(dikvik + dikvki) + a (riky*riky/rik2)*(dikdikvik + dikdikvki) c firstzz = (1./rik - rikz*rikz/rik3)*(dikvik + dikvki) + a (rikz*rikz/rik2)*(dikdikvik + dikdikvki) c firstxy = -(rikx*riky/rik3)*(dikvik + dikvki) + a (rikx*riky/rik2)*(dikdikvik + dikdikvki) c firstxz = -(rikx*rikz/rik3)*(dikvik + dikvki) + a (rikx*rikz/rik2)*(dikdikvik + dikdikvki) c firstyx = -(riky*rikx/rik3)*(dikvik + dikvki) + a (riky*rikx/rik2)*(dikdikvik + dikdikvki) c firstyz = -(riky*rikz/rik3)*(dikvik + dikvki) + a (riky*rikz/rik2)*(dikdikvik + dikdikvki) c firstzx = -(rikz*rikx/rik3)*(dikvik + dikvki) + a (rikz*rikx/rik2)*(dikdikvik + dikdikvki) c firstzy = -(rikz*riky/rik3)*(dikvik + dikvki) + a (rikz*riky/rik2)*(dikdikvik + dikdikvki) c lloop3xx = 0.0 lloop3yy = 0.0 lloop3zz = 0.0 lloop3xy = 0.0 lloop3xz = 0.0 lloop3yx = 0.0 lloop3yz = 0.0 lloop3zx = 0.0 lloop3zy = 0.0 c do 1550 l = 1, nn(i) if (l.eq.k) goto 1550 ril = r(i,l) rilx = x(i,l) rily = y(i,l) rilz = z(i,l) nloop2 = 0.0 do 1580 m = 1, nn(i) if ((m.eq.k).or.(m.eq.l)) goto 1580 dildikvin = pf3(i,m)*zet2(i,m,l)*zet2(i,m,k) nloop2 = nloop2 + dildikvin 1580 continue c dikvil = ff2(i,l,k) dikdikvil = ff22(i,l,k) dikdilvik = ff12(i,k,l) dikdilvil = ff12(i,l,k) c lloop3xx = lloop3xx + (1./rik - rikx*rikx/rik3)*dikvil + a (rikx*rikx/rik2)*dikdikvil + a (rikx*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3yy = lloop3yy + (1./rik - riky*riky/rik3)*dikvil + a (riky*riky/rik2)*dikdikvil + a (riky*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3zz = lloop3zz + (1./rik - rikz*rikz/rik3)*dikvil + a (rikz*rikz/rik2)*dikdikvil + a (rikz*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3xy = lloop3xy - (rikx*riky/rik3)*dikvil + a (rikx*riky/rik2)*dikdikvil + a (rikx*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3xz = lloop3xz - (rikx*rikz/rik3)*dikvil + a (rikx*rikz/rik2)*dikdikvil + a (rikx*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3yx = lloop3yx - (riky*rikx/rik3)*dikvil + a (riky*rikx/rik2)*dikdikvil + a (riky*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3yz = lloop3yz - (riky*rikz/rik3)*dikvil + a (riky*rikz/rik2)*dikdikvil + a (riky*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3zx = lloop3zx - (rikz*rikx/rik3)*dikvil + a (rikz*rikx/rik2)*dikdikvil + a (rikz*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c lloop3zy = lloop3zy - (rikz*riky/rik3)*dikvil + a (rikz*riky/rik2)*dikdikvil + a (rikz*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2) c 1550 continue lloop4xx = 0.0 lloop4yy = 0.0 lloop4zz = 0.0 lloop4xy = 0.0 lloop4xz = 0.0 lloop4yx = 0.0 lloop4yz = 0.0 lloop4zx = 0.0 lloop4zy = 0.0 c do 1600 l = 1, nn(kid) lid = idpart(kid,l) if(lid.eq.i) goto 1600 rkl = r(kid,l) rkl2 = rkl*rkl rklx = x(kid,l) rkly = y(kid,l) rklz = z(kid,l) rmuil = (x(kid,negi)*x(kid,l) + y(kid,negi)*y(kid,l) + a z(kid,negi)*z(kid,l))/(rik*rkl) ril2 = rik2 + rkl2 -2.0*rmuil*rik*rkl ril = sqrt(ril2) ril3 = ril2*ril rilx = -1.0*(x(kid,negi) - x(kid,l)) rily = -1.0*(y(kid,negi) - y(kid,l)) rilz = -1.0*(z(kid,negi) - z(kid,l)) c nloop4xx = 0.0 nloop4yy = 0.0 nloop4zz = 0.0 nloop4xy = 0.0 nloop4xz = 0.0 nloop4yx = 0.0 nloop4yz = 0.0 nloop4zx = 0.0 nloop4zy = 0.0 c do 1650 m = 1, nn(kid) nid = idpart(kid,m) if ((nid.eq.i) .or.(m.eq.l)) goto 1650 rkn = r(kid,m) rkn2 = rkn*rkn rmuin = (x(kid,negi)*x(kid,m) + y(kid,negi)*y(kid,m) + a z(kid,negi)*z(kid,m))/(rik*rkn) rin2 = rik2 + rkn2 - 2.0*rmuin*rik*rkn rin = sqrt(rin2) rinx = -1.0*(x(kid,negi) - x(kid,m)) riny = -1.0*(y(kid,negi) - y(kid,m)) rinz = -1.0*(z(kid,negi) - z(kid,m)) c234567 dildinvki = pf3(kid,negi)*zet3(kid,negi,l)*zet3(kid,negi,m) c nloop4xx = nloop4xx + (rilx*rinx/(ril*rin))*dildinvki nloop4yy = nloop4yy + (rily*riny/(ril*rin))*dildinvki nloop4zz = nloop4zz + (rilz*rinz/(ril*rin))*dildinvki nloop4xy = nloop4xy + (rilx*riny/(ril*rin))*dildinvki nloop4xz = nloop4xz + (rilx*rinz/(ril*rin))*dildinvki nloop4yx = nloop4yx + (rily*rinx/(ril*rin))*dildinvki nloop4yz = nloop4yz + (rily*rinz/(ril*rin))*dildinvki nloop4zx = nloop4zx + (rilz*rinx/(ril*rin))*dildinvki nloop4zy = nloop4zy + (rilz*riny/(ril*rin))*dildinvki c 1650 continue c dikvkl = ff2(kid,l,negi) dilvki = ff3(kid,negi,l) dilvkl = ff3(kid,l,negi) dildilvkl = ff33(kid,l,negi) dildikvki = ff13(kid,negi,l) dildikvkl = ff23(kid,l,negi) dildilvki = ff33(kid,negi,l) dikdikvkl = ff22(kid,l,negi) c C23456789012345678901234567890123456789012345678901234567890123456789012 c lloop4xx = (1./rik - rikx*rikx/rik3)*dikvkl + a (1./ril -rilx*rilx/ril3)*(dilvki + dilvkl) + a ((rilx*rikx + rikx*rilx)/(rik*ril))*(dildikvki + dildikvkl) + a (rikx*rikx/rik2)*dikdikvkl + a (rilx*rilx/ril2)*(dildilvki + dildilvkl) + a nloop4xx + lloop4xx c lloop4yy = (1./rik - riky*riky/rik3)*dikvkl + a (1./ril -rily*rily/ril3)*(dilvki + dilvkl) + a ((rily*riky + riky*rily)/(rik*ril))*(dildikvki + dildikvkl) + a (riky*riky/rik2)*dikdikvkl + a (rily*rily/ril2)*(dildilvki + dildilvkl) + a nloop4yy + lloop4yy c lloop4zz = (1./rik - rikz*rikz/rik3)*dikvkl + a (1./ril -rilz*rilz/ril3)*(dilvki + dilvkl) + a ((rilz*rikz + rikz*rilz)/(rik*ril))*(dildikvki + dildikvkl) + a (rikz*rikz/rik2)*dikdikvkl + a (rilz*rilz/ril2)*(dildilvki + dildilvkl) + a nloop4zz + lloop4zz c lloop4xy = -(rikx*riky/rik3)*dikvkl a -(rilx*rily/ril3)*(dilvki + dilvkl) + a ((rilx*riky + rikx*rily)/(rik*ril))*(dildikvki + dildikvkl) + a (rikx*riky/rik2)*dikdikvkl + a (rilx*rily/ril2)*(dildilvki + dildilvkl) + a nloop4xy + lloop4xy c lloop4xz = -(rikx*rikz/rik3)*dikvkl a -(rilx*rilz/ril3)*(dilvki + dilvkl) + a ((rilx*rikz + rikx*rilz)/(rik*ril))*(dildikvki + dildikvkl) + a (rikx*rikz/rik2)*dikdikvkl + a (rilx*rilz/ril2)*(dildilvki + dildilvkl) + a nloop4xz + lloop4xz c lloop4yx = -(riky*rikx/rik3)*dikvkl a -(rily*rilx/ril3)*(dilvki + dilvkl) + a ((rily*rikx + riky*rilx)/(rik*ril))*(dildikvki + dildikvkl) + a (riky*rikx/rik2)*dikdikvkl + a (rily*rilx/ril2)*(dildilvki + dildilvkl) + a nloop4yx + lloop4yx c lloop4yz = -(riky*rikz/rik3)*dikvkl a -(rily*rilz/ril3)*(dilvki + dilvkl) + a ((rily*rikz + riky*rilz)/(rik*ril))*(dildikvki + dildikvkl) + a (riky*rikz/rik2)*dikdikvkl + a (rily*rilz/ril2)*(dildilvki + dildilvkl) + a nloop4yz + lloop4yz c lloop4zx = -(rikz*rikx/rik3)*dikvkl a -(rilz*rilx/ril3)*(dilvki + dilvkl) + a ((rilz*rikx + rikz*rilx)/(rik*ril))*(dildikvki + dildikvkl) + a (rikz*rikx/rik2)*dikdikvkl + a (rilz*rilx/ril2)*(dildilvki + dildilvkl) + a nloop4zx + lloop4zx c lloop4zy = -(rikz*riky/rik3)*dikvkl a -(rilz*rily/ril3)*(dilvki + dilvkl) + a ((rilz*riky + rikz*rily)/(rik*ril))*(dildikvki + dildikvkl) + a (rikz*riky/rik2)*dikdikvkl + a (rilz*rily/ril2)*(dildilvki + dildilvkl) + a nloop4zy + lloop4zy c 1600 continue c sumkxx = sumkxx + firstxx + lloop3xx + lloop4xx sumkyy = sumkyy + firstyy + lloop3yy + lloop4yy sumkzz = sumkzz + firstzz + lloop3zz + lloop4zz sumkxy = sumkxy + firstxy + lloop3xy + lloop4xy sumkxz = sumkxz + firstxz + lloop3xz + lloop4xz sumkyx = sumkyx + firstyx + lloop3yx + lloop4yx sumkyz = sumkyz + firstyz + lloop3yz + lloop4yz sumkzx = sumkzx + firstzx + lloop3zx + lloop4zx sumkzy = sumkzy + firstzy + lloop3zy + lloop4zy c 1500 continue DXX(i,i) = 0.50*sumkxx DYY(i,i) = 0.50*sumkyy DZZ(i,i) = 0.50*sumkzz DXY(i,i) = 0.50*sumkxy DXZ(i,i) = 0.50*sumkxz DYX(i,i) = 0.50*sumkyx DYZ(i,i) = 0.50*sumkyz DZX(i,i) = 0.50*sumkzx DZY(i,i) = 0.50*sumkzy c c write(28,*) i,i c write(28,*)'DXX = ',DXX(i,i),'DYY = ',DYY(i,i), c a 'DZZ = ',DZZ(i,i) c write(28,*)'DXY = ',DXY(i,i),'DXZ = ',DXZ(i,i) c write(28,*)'DYX = ',DYX(i,i),'DYZ = ',DYZ(i,i) c write(28,*)'DZX = ',DZX(i,i),'DZY = ',DZY(i,i) c write(28,*) ' ' c C END OF BIG "I" LOOP THROUGH PARTICLES c 1002 continue C C CHECK THE RIGID ROTATION CONDITION C write(29,*)'Rigid rotation gives: Code value:' C write(29,*)'---------------------------------------' C do 2002 i = 37,38 C sumDxx = 0.0 sumDyy = 0.0 sumDzz = 0.0 sumDxy = 0.0 sumDxz = 0.0 sumDyx = 0.0 sumDyz = 0.0 sumDzx = 0.0 sumDzy = 0.0 C do 2003 j = 1, np if(j.eq.i) goto 2003 sumDxx = sumDxx + Dxx(i,j) sumDyy = sumDyy + Dyy(i,j) sumDzz = sumDzz + Dzz(i,j) sumDxy = sumDxy + Dxy(i,j) sumDxz = sumDxz + Dxz(i,j) sumDyx = sumDyx + Dyx(i,j) sumDyz = sumDyz + Dyz(i,j) sumDzx = sumDzx + Dzx(i,j) sumDzy = sumDzy + Dzy(i,j) 2003 continue C write(29,*) 'i = ',i write(29,*)'Dxx: ', -1.0*sumDxx, Dxx(i,i) write(29,*)'Dyy: ', -1.0*sumDyy, Dyy(i,i) write(29,*)'Dzz: ', -1.0*sumDzz, Dzz(i,i) write(29,*)'Dxy: ', -1.0*sumDxy, Dxy(i,i) write(29,*)'Dxz: ', -1.0*sumDxz, Dxz(i,i) write(29,*)'Dyx: ', -1.0*sumDyx, Dyx(i,i) write(29,*)'Dyz: ', -1.0*sumDyz, Dyz(i,i) write(29,*)'Dzx: ', -1.0*sumDzx, Dzx(i,i) write(29,*)'Dzy: ', -1.0*sumDzy, Dzy(i,i) write(29,*) ' ' 2002 continue C C This creates an output file in the same format as Meijie's C output filename = dynamorph.format C C NOTE: THE FOLLOWING IS COMMENTED OUT ONLY TO SAVE RUN TIME AND C DISK SPACE. IT SHOULD BE USED IF THE FILE IS DESIRED. C C c do 3000 i = 1,np c do 3100 j = 1,np c write(19,3111) i,1,j,1,Dxx(i,j) c write(19,3111) i,1,j,2,Dxy(i,j) c write(19,3111) i,1,j,3,Dxz(i,j) c3100 continue c do 3200 j = 1,np c write(19,3111) i,2,j,1,Dyx(i,j) c write(19,3111) i,2,j,2,Dyy(i,j) c write(19,3111) i,2,j,3,Dyz(i,j) c3200 continue c do 3300 j = 1,np c write(19,3111) i,3,j,1,Dzx(i,j) c write(19,3111) i,3,j,2,Dzy(i,j) c write(19,3111) i,3,j,3,Dzz(i,j) c3300 continue c3000 continue c3111 format(4(i4,1x),e15.8) C C C C THE FOLLOWING SECTION CALCULATES THE HESSIAN FROM C THE DXX(I,J),DXY(I,J),ETC, OBTAINED ABOVE. ALL IT C DOES IS FORM HESS(A,B) FROM THE FOLLOWING: C [X1,Y1,Z1,X2,Y2,Z2...XN,YN,ZN]^T*[X1,Y1,Z1,X2,Y2,Z2...XN,YN,ZN] C WHICH GIVES A 3NP X 3NP MATRIX. FOR EXAMPLE, C HESS(1,1) = DXX(1,1), HESS(1,2) = DXY(1,1), HESS(2,9) = DYZ(1,3) C C234567 C do 28 ii = 1,np ntp1=ntyp(ii) do 27 jj = 1,np ntp2=ntyp(jj) ia = 3*(ii-1) ib = 3*(jj-1) HESS(1+ia,1+ib) = Dxx(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(1+ia,2+ib) = Dxy(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(1+ia,3+ib) = Dxz(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(2+ia,1+ib) = Dyx(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(2+ia,2+ib) = Dyy(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(2+ia,3+ib) = Dyz(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(3+ia,1+ib) = Dzx(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(3+ia,2+ib) = Dzy(ii,jj)/sqm(ntp1)/sqm(ntp2) HESS(3+ia,3+ib) = Dzz(ii,jj)/sqm(ntp1)/sqm(ntp2) c if ((sqm(ntp1).eq.0).or.(sqm(ntp2).eq.0)) then write(*,*) 'something is wrong with sqms' write(*,*) sqm(ntp1),sqm(ntp2),ntp1,ntp2 endif c 27 continue 28 continue c do 78 ii = 1,648 c do 79 jj = 1,648 c write(91,*) hess(ii,jj) c 911 format(8(1x,f9.5)) c 912 format(f9.5) c 79 continue c 78 continue c234567 C ******************************** call jacobi(hess,mm,nmax,eval) do 78 i = 1,648 if(eval(i).lt.0.0) eval(i) = 0.0 eval(i)=sqrt(eval(i))*2.32525 write(91,*) eval(i) 78 continue C************************************ C end c SUBROUTINE jacobi(a,n,npp,d) INTEGER n,npp,nrot,NMAX real a(npp,npp),d(npp) PARAMETER (NMAX=1000) INTEGER i,ip,iq,j real c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) c do 12 ip=1,n c do 11 iq=1,n c v(ip,iq)=0.d0 c11 continue c v(ip,ip)=1.d0 c12 continue do 13 ip=1,n b(ip)=a(ip,ip) d(ip)=b(ip) z(ip)=0.d0 13 continue nrot=0 do 24 i=1,50 sm=0.d0 do 15 ip=1,n-1 do 14 iq=ip+1,n sm=sm+abs(a(ip,iq)) 14 continue 15 continue c if(sm.eq.0.)return if(sm.lt.1.e-5) return if(i.lt.4)then tresh=0.2*sm/n**2 else tresh=0. endif do 22 ip=1,n-1 do 21 iq=ip+1,n g=100.*abs(a(ip,iq)) if((i.gt.4).and.(abs(d(ip))+ *g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then a(ip,iq)=0. else if(abs(a(ip,iq)).gt.tresh)then h=d(iq)-d(ip) if(abs(h)+g.eq.abs(h))then t=a(ip,iq)/h else theta=0.5*h/a(ip,iq) t=1./(abs(theta)+sqrt(1.+theta**2)) if(theta.lt.0.)t=-t endif c=1./sqrt(1+t**2) s=t*c tau=s/(1.+c) h=t*a(ip,iq) z(ip)=z(ip)-h z(iq)=z(iq)+h d(ip)=d(ip)-h d(iq)=d(iq)+h a(ip,iq)=0. do 16 j=1,ip-1 g=a(j,ip) h=a(j,iq) a(j,ip)=g-s*(h+g*tau) a(j,iq)=h+s*(g-h*tau) 16 continue do 17 j=ip+1,iq-1 g=a(ip,j) h=a(j,iq) a(ip,j)=g-s*(h+g*tau) a(j,iq)=h+s*(g-h*tau) 17 continue do 18 j=iq+1,n g=a(ip,j) h=a(iq,j) a(ip,j)=g-s*(h+g*tau) a(iq,j)=h+s*(g-h*tau) 18 continue c do 19 j=1,n c g=v(j,ip) c h=v(j,iq) c v(j,ip)=g-s*(h+g*tau) c v(j,iq)=h+s*(g-h*tau) c19 continue nrot=nrot+1 endif 21 continue 22 continue do 23 ip=1,n b(ip)=b(ip)+z(ip) d(ip)=b(ip) z(ip)=0. 23 continue 24 continue pause 'too many iterations in jacobi' return END