C ------------------------------------------------------------------------ C PROGRAM MultiChannel (double complex version) C 1. USE TERSOFF POTENTIAL C 2. CALCULATE DYNAMICAL MATRIX FOR ARBITRARY CONFIGURATION C 3. USE MULTICHANNEL PERTURBATION AND SUPERCELL K METHOD TO C GET THE VIBRATIONAL LDOS FOR ANY ATOM IN ANY DIRECTION C 4. HIGH-PRECISION (ORDER-12) INTEGRATION OF EQUATION OF MOTION C 5. "SPIN-WAVE" ENCODED AMPLITUDE SERIES TO MINIMIZE THE ERROR C 6. SHIFT THE SPECTRUM FOR BETTER CONVERGENCE IN THE LOW C FREQUENCY END C C SEPT.21.1996 C DEVELOPED BY JU LI AT MIT C ------------------------------------------------------------------------ C USE THE FOLLOWING UNITS: C ENERGY UNIT = eV = 1.602E-19 JOULE C MASS UNIT = AMU = 1.66054E-27 KG C LENGTH UNIT = ANGSTROM C TIME UNIT = 1.018043E-14 SECOND = 10.18043 FS C STRESS UNIT = 1.602192 MBAR C ------------------------------------------------------------------------- PROGRAM MultiChannel IMPLICIT REAL(A-H,O-Z) PARAMETER (NM=4200,NXX=10,NX=10,ND=80,TOL=1E-5) COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM), A NTYPE(NM),LIST(NM,NXX),INDEX(NM),PMASS(NM),X(NM), A Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM), A DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE, A ZOLD(NM),LIST2(NM,NX),INDEX2(NM),N,MODE,NEW, A HSUM(3,3),RRX(2,NM),RRY(2,NM),RRZ(2,NM),WRITE_LP, A POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0, A VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2), A PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), A C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX), A RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST, A CUBE2,CHI_SI_C,FCUBE2 COMMON /BLOCK3/ IDX(NM*3),ID(ND,NM*3),HESS(ND,NM*3),PI,A(ND,NM*3) DOUBLE COMPLEX A DOUBLE PRECISION PI LOGICAL WRITE_LP,UPDATE,NEW PARAMETER (MM=500,MIP=30000) DOUBLE PRECISION SX(NM),SY(NM),SZ(NM),S(NM),H(3,3),HH(3,3),DELTA, A HA(NM*3,3),HB(NM*3,3),Q(3,3),HBIJ,vec,vx,vy,vz,W,MIN_W, A MAX_W,DEL2,DEL3,DEL4,DEL5,DEL6,AA,DEL7,DEL8,DEL9,DEL10, A DEL11,PHASE,D_PHASE,ALPHA,SUM_WEIGHT,DDSIN(MIP),DDCOS(MIP), A VA(MM),VB(MM),VC(MM),VL(MM),VM(MM),VN(MM),VALUE(MM),WEIGHT, A W1(MM),W2(MM),W4(MM),W6(MM),W8(MM),W10(MM),DOS(3,2000),FPM, A KX,KY,KZ,ar(MM),ai(MM),damax,phi(mm),shift_w,ratio,uu,vv DOUBLE COMPLEX U(3*NM),U2(3*NM),U4(3*NM),U6(3*NM),U8(3*NM),u8i,cc, A U10(3*NM),UOLD(3*NM),FUP(MM),ui,u1i,u2i,u3i,u4i,u5i,u6i,u7i, A u9i,u10i,u11i,uni,un1i,un2i,un3i,un4i,un5i,un6i,un7i,un8i,ii, A un9i,un10i,un11i,uoi,uo2i,uo4i,uo6i,uo8i,uo10i,b(MM),zero,one INTEGER IBASE(200) LOGICAL CONVERGED(MM),ALL_CONVERGED,ALREADY DATA LP,LP_CONFIG,LP_DOS /25,26,27/ CHARACTER *70 BUF # define _ORDER 12 double precision seconds_spent external subroutine startwatch, checkwatch PI = dacos(-1.d0) II = DCMPLX(0.d0,1.d0) ONE = DCMPLX(1.d0,0.d0) ZERO = DCMPLX(0.d0,0.d0) UENERGY_IN_JOULE = 1.602192E-19 ULENGTH_IN_METER = 1.E-10 UMASS_IN_KG = 1.66053E-27 UTIME_IN_SECOND = A ULENGTH_IN_METER*SQRT(UMASS_IN_KG/UENERGY_IN_JOULE) UTIME_IN_PS = UTIME_IN_SECOND/1.E-12 UTIME_IN_FS = UTIME_IN_SECOND/1.E-15 USTRESS_IN_PASCAL = UENERGY_IN_JOULE/ULENGTH_IN_METER**3 USTRESS_IN_BAR = USTRESS_IN_PASCAL/1.E5 USTRESS_IN_MBAR = USTRESS_IN_PASCAL/1.E11 UFREQ_IN_HZ = 1./UTIME_IN_SECOND/2/pi UFREQ_IN_THZ = UFREQ_IN_HZ / 1E12 AVO = 6.0221367E+23 BOLZ = 1.380658E-23 HBAR = 1.054589E-34 PM(1) = 28.086 PM(2) = 12.011 C SET OF PARAMTERS FOR Si. PA(1) = 1.8308E3 PB(1) = 4.7118E2 PL1(1) = 2.4799 PL2(1) = 1.7322 PN(1) = 7.8734E-1 PBT(1) = (1.1E-6)**PN(1) PC(1) = (1.0039E5)**2.E0 PD(1) = (1.6217E1)**2.E0 PH(1) = -5.9825E-1 PR(1) = 2.70 PS(1) = 3.00 C2D2(1) = 1.E0+PC(1)/PD(1) PWR(1) = -1.E0/(2.E0*PN(1)) PWR1(1) = PWR(1)-1.E0 PWR2(1) = PN(1)-1.E0 C SET OF PARAMTERS FOR C. PA(2) = 1.5448E3 PB(2) = 3.8963E2 PL1(2) = 3.4653 PL2(2) = 2.3064 PN(2) = 9.9054E-1 PBT(2) = (4.1612E-6)**PN(2) PC(2) = (1.9981E4)**2.E0 PD(2) = (7.0340)**2.E0 PH(2) = -3.9953E-1 PR(2) = 1.80 PS(2) = 2.10 C2D2(2) = 1.E0+PC(2)/PD(2) PWR(2) = -1.E0/(2.E0*PN(2)) PWR1(2) = PWR(2)-1.E0 PWR2(2) = PN(2)-1.E0 CHI_SI_C = 1.0086 RCUT = 3.00 C for static problems RLIST = RCUT UPDATE = .TRUE. c ====================================================================== c Program input parameters READ (*,'(A70)') buf C "Start from configuration file config(1) or start anew(0):" READ *,LAST READ (*,'(/,A70)') buf C "Lattice constant in Angstrom:" READ *,CUBE READ (*,'(/,A70)') buf C "N: Number of atoms in the supercell:" READ *,N READ (*,'(/,A70,/,A70,/,A70,/,A70)') buf,buf,buf,buf C "Mode of operation: C (0 Stop at maximum iteration MAXITER) C (1 Stop when the base frequency channel converges) C (2 Stop when all frequency channels converge) READ *,MODE READ (*,'(/,A70,/,A70,/,A70)') buf,buf,buf C "Number of supercell k-points that you want to sample C using Monte Carlo (1 and 4 would invoke Cohen-Chadi C k-points for simple cubic as default): " READ *,NSAMPLE READ (*,'(/,A70)') buf C "Random number generator seed:" READ *,ISEED c initiate random number generator ran1 do i=1,mod(iseed,63346) xx = ran1(iseed) enddo READ (*,'(/,A70)') buf C "LDOS atom index (1-N):" READ *,IATOM READ (*,'(/,A70)') buf C "Cartesian direction on that atom:" READ *,Q(1,1),Q(1,2),Q(1,3) READ (*,'(/,A70)') buf C "Number of frequency groups" READ *,NPA READ (*,'(/,A70)') buf C "Mimimum(base) frequency in THz:" READ *,MIN_W READ (*,'(/,A70)') buf C "Maximum frequency(THz) of the spectrum:" READ *,MAX_W READ (*,'(/,A70)') buf C "Shift the spectrum by (THz):" READ *,SHIFT_W C make such that the 2nd freq. is ratio times the 1st freq. w1(1) = (int(shift_w/min_w)+1)*min_w w1(2) = w1(1)+min_w ratio = 2.5 shift_w = dsqrt(((ratio*w1(1))**2-w1(2)**2)/(ratio**2-1)) MAX_W = dsqrt(MAX_W**2+SHIFT_W**2) READ (*,'(/,A70)') buf C "KTWIST: twist period of the "spin-wave" amplitude series:" READ *, KTWIST READ (*,'(/,A70)') buf C "Number of steps of a period for the maximum frequency:" READ *, MINCYCLE READ (*,'(/,A70,/,A70)') buf,buf C "MINITER: number of iterations after which to start C monitoring convergence:" READ *,MINITER READ (*,'(/,A70)') buf C "MAXITER: maximum number of iterations:" READ *,MAXITER READ (*,'(/,A70)') buf C "Fluctuation ratio in three steps to decide convergence:" READ *,TOLE READ (*,'(/,A70)') buf C "Lower bound below which (after MINITER) the LDOS is considered 0:" READ *,TOLM c ====================================================================== NC = NINT((N/8.)**(1.E0/3.E0)) NOLD = N IF (LAST.EQ.1) THEN OPEN (UNIT=LP_CONFIG, STATUS='UNKNOWN', FORM='UNFORMATTED', a FILE="config") READ (LP_CONFIG) N READ (LP_CONFIG) (NTYPE(I),I=1,N) READ (LP_CONFIG) ((H(I,J),J=1,3),I=1,3) READ (LP_CONFIG) (SX(I),SY(I),SZ(I),I=1,N) CLOSE (LP_CONFIG) ELSE C -------------------------------------------------- C DIAMOND CUBIC STRUCTURE AA = 1.E0/FLOAT(NC) X(1) = 0.E0 Y(1) = 0.E0 Z(1) = 0.E0 X(3) = 1.E0/2.E0 Y(3) = 1.E0/2.E0 Z(3) = 0.E0 X(5) = 0.E0 Y(5) = 1.E0/2.E0 Z(5) = 1.E0/2.E0 X(7) = 1.E0/2.E0 Y(7) = 0.E0 Z(7) = 1.E0/2.E0 IP = 0 DO 217 I=1,NC DO 217 J=1,NC DO 217 K=1,NC DO 218 L=1,7,2 SX(IP+L) = (X(L)+I-7.E0/8.E0)*AA SY(IP+L) = (Y(L)+J-7.E0/8.E0)*AA SZ(IP+L) = (Z(L)+K-7.E0/8.E0)*AA NTYPE (IP+L) = 1 SX(IP+L+1) = (X(L)+I-5.E0/8.E0)*AA SY(IP+L+1) = (Y(L)+J-5.E0/8.E0)*AA SZ(IP+L+1) = (Z(L)+K-5.E0/8.E0)*AA NTYPE (IP+L+1) = 2 218 CONTINUE IP = IP+8 217 CONTINUE C -------------------------------------------------- ENDIF CUBE = CUBE*NC CUBE2 = CUBE*0.5 FCUBE2 = -CUBE2 DO 52 I=1,N X(I) = CUBE*SX(I) Y(I) = CUBE*SY(I) Z(I) = CUBE*SZ(I) PMASS(I) = PM(NTYPE(I)) 52 CONTINUE IF (UPDATE) THEN C --------------------------------------------------------- C UPDATE THE VERLET LIST DO 57 IPT = 1,N INDEX(IPT) = 0 IF (IPT.LE.2) THEN RRX(IPT,IPT) = 0. RRY(IPT,IPT) = 0. RRZ(IPT,IPT) = 0. ENDIF 57 CONTINUE DO 170 IPT = 1,N N1 = NTYPE(IPT) XI = X(IPT) YI = Y(IPT) ZI = Z(IPT) DO 270 JPT = IPT+1,N N2 = NTYPE(JPT) RXX = X(JPT) - XI RYY = Y(JPT) - YI RZZ = Z(JPT) - ZI IF(RXX.GT.CUBE2) RXX=RXX-CUBE IF(RXX.LT.FCUBE2) RXX=RXX+CUBE IF(RYY.GT.CUBE2) RYY=RYY-CUBE IF(RYY.LT.FCUBE2) RYY=RYY+CUBE IF(RZZ.GT.CUBE2) RZZ=RZZ-CUBE IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE R1 = SQRT(RXX*RXX+RYY*RYY+RZZ*RZZ) C STORE THE RELATIVE POSITION FOR PHONON DISPERSION CURVE IF (IPT.LE.2) THEN RRX(IPT,JPT) = RXX RRY(IPT,JPT) = RYY RRZ(IPT,JPT) = RZZ ENDIF IF (R1.LT.RLIST) THEN INDEX(IPT)=INDEX(IPT)+1 LIST(IPT,INDEX(IPT))=JPT INDEX(JPT)=INDEX(JPT)+1 LIST(JPT,INDEX(JPT))=IPT ENDIF 270 CONTINUE 170 CONTINUE RRX(2,1) = -RRX(1,2) RRY(2,1) = -RRY(1,2) RRZ(2,1) = -RRZ(1,2) C --------------------------------------------------------- ENDIF call tersoff_dynamical print *,' dynamical matrix construction complete.' print *,' ' C --------------------------------------------------------- C Multi-Channel Perturbation Method N3 = 3*N AA = 2*PI/CUBE IBB = 3*(IATOM-1) c ================================================================= c do a unitary transformation on the dynamical matrix such that c our defined degree of freedom is a normal coordinate. vec = dsqrt(Q(1,1)**2+Q(1,2)**2+Q(1,3)**2) Q(1,1) = Q(1,1)/vec Q(1,2) = Q(1,2)/vec Q(1,3) = Q(1,3)/vec vx = ran1(iseed) vy = ran1(iseed) vz = ran1(iseed) c take a cross-product V*C1 Q(2,1) = vy*Q(1,3)-Q(1,2)*vz Q(2,2) = Q(1,1)*vz-vx*Q(1,3) Q(2,3) = vx*Q(1,2)-Q(1,1)*vy vec = dsqrt(Q(2,1)**2+Q(2,2)**2+Q(2,3)**2) Q(2,1) = Q(2,1)/vec Q(2,2) = Q(2,2)/vec Q(2,3) = Q(2,3)/vec c take the cross-product C1*C2 Q(3,1) = Q(1,2)*Q(2,3)-Q(2,2)*Q(1,3) Q(3,2) = Q(2,1)*Q(1,3)-Q(1,1)*Q(2,3) Q(3,3) = Q(1,1)*Q(2,2)-Q(2,1)*Q(1,2) do j=1,3 c scatter the part of sparse matrix HESS concerning c atom IATOM into an N3 by 3 matrix HA do 730 i=1,N3 730 HA(i,j) = 0. do 76 k=1,idx(j+ibb) 76 HA(id(k,j+ibb),j) = HESS(k,j+ibb) enddo c HB = HA*QT do 8776 i=1,N3 do 8776 j=1,3 HB(i,j) = 0. do 8776 k=1,3 8776 HB(i,j) = HB(i,j)+HA(i,k)*Q(j,k) c get the middle part of HA do 987 i=1,3 do 987 j=1,3 987 HA(i,j) = HB(i+ibb,j) c multiply by Q, put back into HB do 988 i=1,3 do 988 j=1,3 HB(ibb+i,j) = 0. do 988 k=1,3 988 HB(ibb+i,j) = HB(ibb+i,j)+Q(i,k)*HA(k,j) c re-insert the transformed matrix into HESS do 876 j=ibb+1,ibb+3 idx(j) = 0 do 876 i=1,N3 HBIJ = HB(i,j-ibb) c first insert the columns if (dabs(HBIJ).lt.tol) goto 259 idx(j) = idx(j)+1 id(idx(j),j) = i HESS(idx(j),j) = HBIJ c then symmetrically insert the rows. 259 if (i.ge.ibb+1.and.i.le.ibb+3) goto 876 do k=1,idx(i) if (j.eq.id(k,i)) then hess(k,i) = HBIJ goto 876 endif enddo if (dabs(HBIJ).lt.tol) goto 876 idx(i) = idx(i)+1 id(idx(i),i)=j hess(idx(i),i) = HBIJ 876 continue c ================================================================= IBB = IBB+1 C Count the number of non-zero elements in the sparse matrix: nonzero = 0 do I=1,N3 nonzero = nonzero+IDX(I) enddo call startwatch() IP = 0 DO 765 I=1,NPA c construct base frequency within (min_w, 2*min_w) IP = IP+1 IF (I.EQ.1) THEN DOS(1,IP) = MIN_W ELSE DOS(1,IP) = MIN_W+RAN1(ISEED)*MIN_W ENDIF DOS(2,IP) = 0.D0 IBASE(I) = IP ALPHA = DOS(1,IP) c construct multiples of this base frequency within max_w W = 2.D0*ALPHA 155 IF (W.GT.MAX_W) GOTO 765 IP = IP+1 DOS(1,IP) = W DOS(2,IP) = 0.D0 W = W+ALPHA GOTO 155 765 CONTINUE IBASE(NPA+1) = IP+1 DO 764 IW=1,NPA PRINT *,'frequency group',iw,' in THz' DO 763 J=IBASE(IW),IBASE(IW+1)-1 IF (DOS(1,J).LT.SHIFT_W) THEN DOS(3,J) = DOS(1,J)-SHIFT_W ELSE c store the true frequency in dos(3,j) DOS(3,J) = DSQRT(DOS(1,J)**2-SHIFT_W**2) ENDIF PRINT *,'true w = ',DOS(3,J) 763 CONTINUE ALPHA = DOS(1,IBASE(IW))/UFREQ_IN_THZ MMAX = IBASE(IW+1)-IBASE(IW) MINT = MMAX*MINCYCLE print *,'Base Cycle =',mint print *,' ' 764 continue SUM_WEIGHT = 0.D0 DO 209 IK=1,NSAMPLE C --------------------------------------------------------- C Monte Carlo or special-point sampling of k space for c the supercell IF (NSAMPLE.EQ.1) THEN KX = AA*1.D0/4.D0 KY = AA*1.D0/4.D0 KZ = AA*1.D0/4.D0 WEIGHT = 1.D0 ELSEIF (NSAMPLE.EQ.4) THEN IF (IK.EQ.1) THEN KX = AA*3.D0/8.D0 KY = AA*1.D0/8.D0 KZ = AA*1.D0/8.D0 WEIGHT = 3.D0/8.D0 ELSEIF (IK.EQ.2) THEN KX = AA*3.D0/8.D0 KY = AA*3.D0/8.D0 KZ = AA*1.D0/8.D0 WEIGHT = 3.D0/8.D0 ELSEIF (IK.EQ.3) THEN KX = AA*1.D0/8.D0 KY = AA*1.D0/8.D0 KZ = AA*1.D0/8.D0 WEIGHT = 1.D0/8.D0 ELSEIF (IK.EQ.4) THEN KX = AA*3.D0/8.D0 KY = AA*3.D0/8.D0 KZ = AA*3.D0/8.D0 WEIGHT = 1.D0/8.D0 ENDIF ELSE KX = AA*(RAN1(ISEED)-0.5) KY = AA*(RAN1(ISEED)-0.5) KZ = AA*(RAN1(ISEED)-0.5) WEIGHT = 1.D0/NSAMPLE ENDIF print *,' ' print *,'wavevector',IK print *,'(',kx/AA,ky/AA,kz/AA,'): weight=', weight print *,' ' SUM_WEIGHT = SUM_WEIGHT+WEIGHT DO 83 I=1,N3 IPT = (I-1)/3+1 XI = X(IPT) YI = Y(IPT) ZI = Z(IPT) ALREADY = .FALSE. DO 82 K=1,IDX(I) J = ID(K,I) JPT = (J-1)/3+1 RXX = X(JPT) - XI RYY = Y(JPT) - YI RZZ = Z(JPT) - ZI IF(RXX.GT.CUBE2) RXX=RXX-CUBE IF(RXX.LT.FCUBE2) RXX=RXX+CUBE IF(RYY.GT.CUBE2) RYY=RYY-CUBE IF(RYY.LT.FCUBE2) RYY=RYY+CUBE IF(RZZ.GT.CUBE2) RZZ=RZZ-CUBE IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE PHASE = KX*RXX+KY*RYY+KZ*RZZ A(K,I) = HESS(K,I)*DCMPLX(DCOS(PHASE),DSIN(PHASE)) c shift the spectrum by (shift_w)**2 IF (I.EQ.J.AND.(.NOT.ALREADY)) THEN A(K,I) = A(K,I)+(SHIFT_W/UFREQ_IN_THZ)**2 ALREADY = .TRUE. ENDIF 82 CONTINUE IF (.NOT.ALREADY) THEN IDX(I) = IDX(I)+1 ID(IDX(I),I) = I A(IDX(I),I) = (SHIFT_W/UFREQ_IN_THZ)**2 ENDIF 83 CONTINUE DO 462 IW = 1,NPA c different frequency groups ALPHA = DOS(1,IBASE(IW))/UFREQ_IN_THZ MMAX = IBASE(IW+1)-IBASE(IW) MINT = MMAX*MINCYCLE DELTA = 2.D0*dacos(-1.D0)/ALPHA/MINT DEL2 = DELTA**2/2.D0 DEL3 = DELTA**3/6.D0 DEL4 = DELTA**4/24.D0 DEL5 = DELTA**5/120.D0 DEL6 = DELTA**6/720.D0 DEL7 = DELTA**7/5040.D0 DEL8 = DELTA**8/40320.D0 DEL9 = DELTA**9/362880.D0 DEL10 = DELTA**10/3628800.D0 DEL11 = DELTA**11/39916800.D0 D_PHASE = 2*dacos(-1.D0)/MINT DO IP=1,MINT PHASE = D_PHASE*IP DDSIN(IP) = DSIN(PHASE) DDCOS(IP) = DCOS(PHASE) ENDDO DO M = 1,MMAX VA(M) = 1.D0 VB(M) = 2.D0 VC(M) = 5.D0 VL(M) = 1.D0 VM(M) = 1.D0 VN(M) = 1.D0 FUP(M) = ZERO W1(M) = M*ALPHA W2(M) = (M*ALPHA)**2 W4(M) = (M*ALPHA)**4 W6(M) = (M*ALPHA)**6 W8(M) = (M*ALPHA)**8 W10(M)= (M*ALPHA)**10 CONVERGED(M) = .FALSE. ENDDO c ---------------------------------------------------------- c To generate an error-cancelling amplitude series. c ++--++--++--++-- ar(1) = 1.d0 ar(2) = 1.d0 if (ran1(iseed).gt.0.5) ar(2) = -ar(2) do m=3,mmax ar(m)= -ar(m-2) enddo c +-+-+-+-+-+-+-+- parity = 1 if (ran1(iseed).gt.0.5) parity = -parity do m=1,mmax ai(m) = parity*ar(m) parity = -parity enddo c clear the empty channels do m=1,mmax if (dos(1,ibase(iw)+m-1).gt.SHIFT_W) goto 974 ar(m) = 0. ai(m) = 0. enddo 974 if (m.gt.2) then ar(m-2) = ar(m+2) ai(m-2) = ai(m+2) endif c randomly generate a "spin-wave" to avoid numerical c problems in real time phi(1) = 0.d0 damax = 2*dacos(-1.d0)/ktwist if (ran1(iseed).gt.0.5) damax = -damax do m=1,mmax b(m) = DCMPLX(ar(m),ai(m))*DCMPLX(dcos(phi(m)),dsin(phi(m))) phi(m+1)=phi(m)+damax*ran1(iseed) enddo c ---------------------------------------------------------- c pre-condition the integration: c only odd-power derivatives exist at t=0 c to save space, u2 stands for u3, u4 stands for u5 etc. do 612 i=1,n3 u(i) = ZERO u2(i) = ZERO u4(i) = ZERO u6(i) = ZERO u8(i) = ZERO u10(i) = ZERO 612 continue do m=1,mmax u2(ibb) = u2(ibb)+m*alpha*b(m) u4(ibb) = u4(ibb)-(m*alpha)**3*b(m) u6(ibb) = u6(ibb)+(m*alpha)**5*b(m) u8(ibb) = u8(ibb)-(m*alpha)**7*b(m) u10(ibb) = u10(ibb)+(m*alpha)**9*b(m) enddo do i=1,n3 do k=1,idx(i) jd = id(k,i) u4(jd) = u4(jd)-A(k,i)*u2(i) enddo enddo do i=1,n3 do k=1,idx(i) jd = id(k,i) u6(jd) = u6(jd)-A(k,i)*u4(i) enddo enddo do i=1,n3 do k=1,idx(i) jd = id(k,i) u8(jd) = u8(jd)-A(k,i)*u6(i) enddo enddo do i=1,n3 do k=1,idx(i) jd = id(k,i) u10(jd) = u10(jd)-A(k,i)*u8(i) enddo enddo do i=1,n3 u(i) = u2(i)*del3+u4(i)*del5+u6(i)*del7 c +u8(i)*del9+u10(i)*del11 uold(i) = ZERO enddo u1i = ZERO u3i = u2(ibb) u5i = u4(ibb) u7i = u6(ibb) u9i = u8(ibb) u11i = u10(ibb) ui = ZERO u2i = ZERO u4i = ZERO u6i = ZERO u8i = ZERO u10i = ZERO ip = 1 576 ipp = mod(ip-1,mint)+1 #ifndef _NO_PRINT_IP print *, 'ip = ',ip #endif do i=1,n3 u2(i) = ZERO u4(i) = ZERO u6(i) = ZERO u8(i) = ZERO u10(i) = ZERO enddo c add perturbation do m=1,mmax ipm = mod(m*ipp-1,mint)+1 u2(ibb) = u2(ibb)+ddsin(ipm)*b(m) u4(ibb) = u4(ibb)-w2(m)*ddsin(ipm)*b(m) u6(ibb) = u6(ibb)+w4(m)*ddsin(ipm)*b(m) u8(ibb) = u8(ibb)-w6(m)*ddsin(ipm)*b(m) u10(ibb)= u10(ibb)+w8(m)*ddsin(ipm)*b(m) enddo do i=1,n3 cc = u(i) idxi=idx(i) do k=1,idxi jd = id(k,i) u2(jd) = u2(jd)-A(k,i)*cc enddo enddo do i=1,n3 cc = u2(i) idxi=idx(i) do k=1,idxi jd = id(k,i) u4(jd) = u4(jd)-A(k,i)*cc enddo enddo do i=1,n3 cc = u4(i) idxi=idx(i) do k=1,idxi jd = id(k,i) u6(jd) = u6(jd)-A(k,i)*cc enddo enddo do i=1,n3 cc = u6(i) idxi=idx(i) do k=1,idxi jd = id(k,i) u8(jd) = u8(jd)-A(k,i)*cc enddo enddo do i=1,n3 cc = u8(i) idxi=idx(i) do k=1,idxi jd = id(k,i) u10(jd) = u10(jd)-A(k,i)*cc enddo enddo if (ip.ne.1) then u11i = (u10(ibb)-uo10i)/2.d0/delta u9i = ((u8(ibb)-uo8i)/2.d0-del3*u11i)/delta u7i = ((u6(ibb)-uo6i)/2.d0-del3*u9i-del5*u11i)/delta u5i = ((u4(ibb)-uo4i)/2.d0-del3*u7i-del5*u9i-del7*u11i)/delta u3i = ((u2(ibb)-uo2i)/2.d0-del3*u5i-del5*u7i-del7*u9i c -del9*u11i)/delta u1i = ((u(ibb)-uoi)/2.d0-del3*u3i-del5*u5i-del7*u7i c -del9*u9i-del11*u11i)/delta endif c -------------------------------------------------------------------- c calculate the integral uni = ui+delta*u1i+del2*u2i+del3*u3i+del4*u4i+del5*u5i c +del6*u6i+del7*u7i+del8*u8i+del9*u9i+del10*u10i+del11*u11i un1i = u1i+delta*u2i+del2*u3i+del3*u4i+del4*u5i+del5*u6i c +del6*u7i+del7*u8i+del8*u9i+del9*u10i+del10*u11i un2i = u2i+delta*u3i+del2*u4i+del3*u5i+del4*u6i+del5*u7i c +del6*u8i+del7*u9i+del8*u10i+del9*u11i un3i = u3i+delta*u4i+del2*u5i+del3*u6i+del4*u7i+del5*u8i c +del6*u9i+del7*u10i+del8*u11i un4i = u4i+delta*u5i+del2*u6i+del3*u7i+del4*u8i+del5*u9i c +del6*u10i+del7*u11i un5i = u5i+delta*u6i+del2*u7i+del3*u8i+del4*u9i+del5*u10i c +del6*u11i un6i = u6i+delta*u7i+del2*u8i+del3*u9i+del4*u10i+del5*u11i un7i = u7i+delta*u8i+del2*u9i+del3*u10i+del4*u11i un8i = u8i+delta*u9i+del2*u10i+del3*u11i un9i = u9i+delta*u10i+del2*u11i un10i = u10i+delta*u11i un11i = u11i do m=1,mmax ipm = mod(m*ipp-1,mint)+1 ibm = mod(m*(ipp+mint-1)-1,mint)+1 fup(m) = fup(m)+(uni-un2i/w2(m)+un4i/w4(m)-un6i/w6(m)+un8i/w8(m) c -un10i/w10(m))/w1(m)*ddsin(ipm) c -(ui-u2i/w2(m)+u4i/w4(m)-u6i/w6(m)+u8i/w8(m) c -u10i/w10(m))/w1(m)*ddsin(ibm) c +(un1i-un3i/w2(m)+un5i/w4(m)-un7i/w6(m)+un9i/w8(m) c -un11i/w10(m))/w2(m)*ddcos(ipm) c -(u1i-u3i/w2(m)+u5i/w4(m)-u7i/w6(m)+u9i/w8(m) c -u11i/w10(m))/w2(m)*ddcos(ibm) enddo c -------------------------------------------------------------------- c total flops = (4+2) * (# of non-zero matrix elements) c * (# of time steps) * (_ORDER-2)/2 #define _A_CHUNK_OF_WORK 500000000 #define _STEP_TO_INFORM _A_CHUNK_OF_WORK/nonzero/(_ORDER-2) IF (mod(ip,_STEP_TO_INFORM).eq.0) THEN call checkwatch(seconds_spent) totalflops = 3. * nonzero * ip * (_ORDER-2) print *,'Step = ',ip,100.*ip/mint,'% of a base cycle done:' floprate = (1.e-6)*totalflops/seconds_spent print *,' Average Mflops/s up to now = ',floprate ENDIF IF (IPP.EQ.MINT) THEN print *,' ' DO 347 M=1,MMAX im = ibase(iw)+m-1 c for empty channels if (ar(m).eq.0.and.ai(m).eq.0) then converged(m) = .true. value(m) = 0. endif if (.not.converged(m)) then c decode the spin-wave fup(m) = fup(m)*DCMPLX(dcos(phi(m)),-dsin(phi(m))) uu = real(fup(m)) vv = dimag(fup(m)) c cancel out all odd-numbered noise functions fpm = (uu+dsign(ar(m)/ai(m),ai(m)*ar(m))*vv)/ar(m)/2.d0 c plug into the formulae VALUE(M) = -2.*M*ALPHA*ALPHA*MINT/IP/pi/pi/UFREQ_IN_THZ*FPM c scale it because of the spectra shift transformation VALUE(M) = VALUE(M)*DOS(3,IM)/DOS(1,IM) print *,'w =',dos(3,im), value(m) else print *,'w =',dos(3,im), dos(2,im)/sum_weight, a ' * converged' endif 347 CONTINUE print *,' ' IF ((IP/MINT).GT.MAXITER) THEN C ------------------------------------------------------------ C stop if after MAXITER base cycles we still did't converge do 705 m=1,mmax IM = IBASE(IW)+M-1 if (.NOT.CONVERGED(M)) a DOS (2,IM)=DOS(2,IM)+ a ((VALUE(M)+VC(M)+VB(M)+VA(M)+VN(M)+VM(M))/6.D0)*WEIGHT 705 continue GOTO 462 ELSEIF (MODE.EQ.0.AND.(IP/MINT).EQ.MAXITER) THEN C ------------------------------------------------------------ C stop after certain iterations of the base cycle do 703 m=1,mmax IM = IBASE(IW)+M-1 if (value(M).gt.tolm) then DOS (2,IM)=DOS(2,IM)+VALUE(M)*WEIGHT endif 703 continue GOTO 462 ELSEIF (MODE.EQ.1) THEN C ------------------------------------------------------------ C converge when the base frequency component converge IF ((ABS(VALUE(1)/VA(1)-1).LT.TOLE.AND. a ABS(VALUE(1)/VB(1)-1).LT.TOLE.AND. a ABS(VALUE(1)/VC(1)-1).LT.TOLE).or. a (value(1).lt.tolm.and.(IP/MINT).gt.MINITER)) THEN do 738 m=1,mmax IM = IBASE(IW)+M-1 if (value(M).gt.tolm) a DOS (2,IM)=DOS(2,IM)+ a (VALUE(M)*0.5+(VC(M)+VB(M))*0.25)*WEIGHT 738 continue GOTO 462 ENDIF ELSEIF (MODE.EQ.2) THEN C ------------------------------------------------------------ C converge when all frequency components converge ALL_CONVERGED = .TRUE. DO 932 M=1,MMAX IF (.NOT.CONVERGED(M)) THEN CONVERGED(M) = (ABS(VALUE(M)/VA(M)-1).LT.TOLE.AND. a ABS(VALUE(M)/VB(M)-1).LT.TOLE.and. a ABS(VALUE(M)/VC(M)-1).LT.TOLE).or. a (value(M).lt.tolm.and.(IP/MINT).gt.MINITER) IF (.NOT.CONVERGED(M)) THEN ALL_CONVERGED = .FALSE. ELSEIF (value(m).gt.tolm) THEN IM = IBASE(IW)+M-1 DOS (2,IM)=DOS(2,IM)+ a (VALUE(M)*0.5+(VC(M)+VB(M))*0.25)*WEIGHT ENDIF ENDIF 932 CONTINUE IF (ALL_CONVERGED) GOTO 462 C ------------------------------------------------------------ ENDIF DO 679 M=1,MMAX VL(M) = VM(M) VM(M) = VN(M) VN(M) = VA(M) VA(M) = VB(M) VB(M) = VC(M) VC(M) = VALUE(M) 679 CONTINUE ENDIF c for a base frequency cycle uo10i = u10i uo8i = u8i uo6i = u6i uo4i = u4i uo2i = u2i uoi = ui u10i = u10(ibb) u8i = u8(ibb) u6i = u6(ibb) u4i = u4(ibb) u2i = u2(ibb) ui = u(ibb) do i=1,n3 cc = uold(i) uold(i) = u(i) u(i) = 2.d0*(u(i)+del2*u2(i)+del4*u4(i)+del6*u6(i) c +del8*u8(i)+del10*u10(i))-cc enddo IP = IP+1 GOTO 576 462 CONTINUE c the loop through all base frequency sets ends here OPEN (UNIT=LP_DOS, STATUS='UNKNOWN', FORM='FORMATTED', a FILE="dos.out") DO I = 1,IBASE(NPA+1)-1 IF (DOS(3,I).GT.0.E0) a WRITE (LP_DOS,501) DOS(3,I), DOS(2,I)/SUM_WEIGHT 501 FORMAT (2(1X,E13.7)) ENDDO CLOSE (LP_DOS) if (ntype(iatom).eq.1) goto 209 c if it's carbon atom then save an extra copy. open (unit=lp_dos, status='unknown', form='formatted', a file="doc.out") do i = 1,ibase(npa+1)-1 if (dos(3,i).gt.0.e0) a write (lp_dos,501) dos(3,i), dos(2,i)/sum_weight enddo close (lp_dos) c sample more supercell k's 209 CONTINUE call checkwatch(seconds_spent) totalflops = 3. * nonzero * ip * (_ORDER-2) print *,'Step = ',ip,100.*ip/mint,'% of a base cycle done:' floprate = (1.e-6)*totalflops/seconds_spent print *,' Average Mflops/s up to now = ',floprate STOP END SUBROUTINE TERSOFF_DYNAMICAL IMPLICIT REAL(A-H,O-Z) PARAMETER (NM=4200,NXX=10,NX=10,ND=80,TOL=1E-5) COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM), A NTYPE(NM),LIST(NM,NXX),INDEX(NM),PMASS(NM),X(NM), A Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM), A DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE, A ZOLD(NM),LIST2(NM,NX),INDEX2(NM),N,MODE,NEW, A HSUM(3,3),RRX(2,NM),RRY(2,NM),RRZ(2,NM),WRITE_LP, A POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0, A VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2), A PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), A C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX), A RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST, A CUBE2,CHI_SI_C,FCUBE2 COMMON /BLOCK3/ IDX(NM*3),ID(ND,NM*3),HESS(ND,NM*3),PI, a A(ND,NM*3),B(ND,NM*3) DOUBLE PRECISION A,B,PI LOGICAL WRITE_LP DIMENSION ZET2(NM,NX,NX),ZET3(NM,NX,NX),FC(NM),FCP(NM), A ZET12(NX,NX),ZET13(NX,NX),ZET22(NX,NX),ZET23(NX,NX), A ZET33(NX,NX),RX(NM,NX),RY(NM,NX),RZ(NM,NX),R(NM,NX), A R2(NM,NX),RXT(NX),RYT(NX),RZT(NX),RJK2(NX,NX), A RJK(NX,NX),RJKX(NX,NX),RJKY(NX,NX),FF11(NM,NX), A RJKZ(NX,NX),PRJKX(NX,NX),PRJKY(NX,NX),PRJKZ(NX,NX), A ZET(NX),SZET1(NX),SZET11(NX),SFTM2X(NX),SFTM2Y(NX), A SFTM2Z(NX),SFTM3X(NX),SFTM3Y(NX),SFTM3Z(NX),PF(NX), A PF2(NX),PF3(NM,NX),FRST(NX),FRST2(NX),FF1(NM,NX), A FF22(NM,NX,NX),FF23(NM,NX,NX),FF13(NM,NX,NX),FF12(NM,NX,NX), A FF3(NM,NX,NX),FF2(NM,NX,NX),FF33(NM,NX,NX),FCPP(NM), A LIST3(NX+NX*NX) real lloopxx,lloopxy,nloopxx,nloopxy, a lloop2xx,lloop2xy,nloop2xx,nloop2xy, a 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 DO 10 IPT=1,N FX(IPT) = 0.E0 FY(IPT) = 0.E0 FZ(IPT) = 0.E0 EP(IPT) = 0.E0 10 CONTINUE DO 500 IPT = 1,N N1 = NTYPE(IPT) XI = X(IPT) YI = Y(IPT) ZI = Z(IPT) NEG = 0 DO 600 J = 1,INDEX(IPT) JPT = LIST(IPT,J) N2 = NTYPE(JPT) RSMAX = PS(N1)*PS(N2) RXX = X(JPT) - XI RYY = Y(JPT) - YI RZZ = Z(JPT) - ZI IF(RXX.GT.CUBE2) RXX=RXX-CUBE IF(RXX.LT.FCUBE2) RXX=RXX+CUBE IF(RYY.GT.CUBE2) RYY=RYY-CUBE IF(RYY.LT.FCUBE2) RYY=RYY+CUBE IF(RZZ.GT.CUBE2) RZZ=RZZ-CUBE IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE RR = RXX*RXX+RYY*RYY+RZZ*RZZ IF(RR.GT.RSMAX) GOTO 600 NEG = NEG + 1 LIST2(IPT,NEG) = JPT R1 = SQRT(RR) RX(IPT,NEG) = RXX RY(IPT,NEG) = RYY RZ(IPT,NEG) = RZZ R2(IPT,NEG) = RR R(IPT,NEG) = R1 RXT(NEG) = RXX/R1 RYT(NEG) = RYY/R1 RZT(NEG) = RZZ/R1 CA(NEG) = SQRT(PA(N1)*PA(N2)) CB(NEG) = SQRT(PB(N1)*PB(N2)) CR(NEG) = SQRT(PR(N1)*PR(N2)) CS(NEG) = SQRT(PS(N1)*PS(N2)) if (n1.eq.1.and.n2.eq.1) then cr(neg) = 1.6 cs(neg) = 2.9 elseif (n1.eq.2.and.n2.eq.2) then cr(neg) = 1.8 cs(neg) = 2.1 else cr(neg) = 2.2 cs(neg) = 2.5 endif RLM1(NEG) = (PL1(N1)+PL1(N2))/2.E0 RLM2(NEG) = (PL2(N1)+PL2(N2))/2.E0 IF(N1.EQ.N2) THEN CHIJ(NEG) = 1.E0 ELSE CHIJ(NEG) = 1.0109 ENDIF IF(R1.LT.CR(NEG)) THEN FC(NEG) = 1.E0 FCP(NEG) = 0.E0 FCPP(NEG) = 0.E0 ELSE Fc(NEG)=0.5E0*(1.E0+COS(PI*(R1-CR(NEG)) A /(CS(NEG)-CR(NEG)))) FCP(NEG)=-0.5E0*PI*SIN(PI*(R1-CR(NEG)) A /(CS(NEG)-CR(NEG)))/(CS(NEG)-CR(NEG)) FCPP(NEG)=(0.5E0-FC(NEG))*PI*PI/(CS(NEG)-CR(NEG)) A /(CS(NEG)-CR(NEG)) ENDIF 600 CONTINUE INDEX2(IPT) = NEG DO 700 J = 1,NEG JPT = LIST2(IPT,J) N22 = NTYPE(JPT) ZET(J) = 0.E0 SZET1(J) = 0.E0 SZET11(J) = 0.E0 SFTM2X(J) = 0.E0 SFTM2Y(J) = 0.E0 SFTM2Z(J) = 0.E0 SFTM3X(J) = 0.E0 SFTM3Y(J) = 0.E0 SFTM3Z(J) = 0.E0 DO 800 K = 1,NEG IF(K.EQ.J) GOTO 800 RMU = (RX(IPT,J)*RX(IPT,K)+RY(IPT,J)*RY(IPT,K)+ a RZ(IPT,J)*RZ(IPT,K))/R(IPT,J)/R(IPT,K) RJK2(J,K) = R2(IPT,J)+R2(IPT,K) a -2.E0*RMU*R(IPT,J)*R(IPT,K) RJK(J,K) = SQRT(RJK2(J,K)) RJKX(J,K) = RX(IPT,J)-RX(IPT,K) RJKY(J,K) = RY(IPT,J)-RY(IPT,K) RJKZ(J,K) = RZ(IPT,J)-RZ(IPT,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) RMU1 = RMU-PH(N1) RMU2 = RMU1*RMU1 RMU3 = PD(N1)+RMU2 GMU = C2D2(N1)-PC(N1)/RMU3 GMUP = 2.E0*RMU1*PC(N1)/RMU3**2 COM1 = FCP(K)*GMU COM2 = FC(K)*GMUP COM3 = FCPP(K)*GMU ZET(J) = ZET(J)+FC(K)*GMU TIP1 = RMU/R(IPT,J)-1.E0/R(IPT,K) TIP2 = RMU/R(IPT,K)-1.E0/R(IPT,J) TIP3 = RJK(J,K)/R(IPT,J)/R(IPT,K) ZET1 = COM2*(-TIP1) SZET1(J) = SZET1(J)+ZET1 ZET2(IPT,J,K) = COM1+COM2*(-TIP2) ZET3(IPT,J,K) = COM2*(-TIP3) GMUPP = 2.E0*PC(N1)/RMU3**2-8.E0*PC(N1)*RMU2/RMU3**3 COM4 = FCP(K)*GMUP COM5 = FC(K)*GMUPP TIP11 = 2.E0*RMU/R2(IPT,J)-1.E0/R(IPT,J)/R(IPT,K) TIP12 = RMU/R(IPT,J)/R(IPT,K)-1.E0/R2(IPT,K)-1.E0/R2(IPT,J) TIP13 = RJK(J,K)/R2(IPT,J)/R(IPT,K) TIP22 = 2.E0*RMU/R2(IPT,K)-1.E0/R(IPT,J)/R(IPT,K) TIP23 = RJK(J,K)/R2(IPT,K)/R(IPT,J) TIP33 = -1.E0/R(IPT,J)/R(IPT,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.E0*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 SFTM2X(J) = SFTM2X(J) + ZET2(IPT,J,K)*RXT(K) SFTM2Y(J) = SFTM2Y(J) + ZET2(IPT,J,K)*RYT(K) SFTM2Z(J) = SFTM2Z(J) + ZET2(IPT,J,K)*RZT(K) SFTM3X(J) = SFTM3X(J) + ZET3(IPT,J,K)*PRJKX(J,K) SFTM3Y(J) = SFTM3Y(J) + ZET3(IPT,J,K)*PRJKY(J,K) SFTM3Z(J) = SFTM3Z(J) + ZET3(IPT,J,K)*PRJKZ(J,K) 800 CONTINUE 700 CONTINUE DO 704 J=1,NEG JPT = LIST2(IPT,J) N22 = NTYPE(JPT) ELM2 = EXP(-RLM2(J)*R(IPT,J)) ELM2B = ELM2*CB(J) ELM1 = EXP(-RLM1(J)*R(IPT,J)) ELM1A = ELM1*CA(J) FCD1 = FCP(J)-RLM1(J)*FC(J) FCD2 = FCP(J)-RLM2(J)*FC(J) FCD1P = FCPP(J)-2.E0*RLM1(J)*FCP(J)+RLM1(J)**2.*FC(J) FCD2P = FCPP(J)-2.E0*RLM2(J)*FCP(J)+RLM2(J)**2.*FC(J) IF(ZET(J).NE.0.E0) THEN BEZET = 1.E0+ PBT(N1)*(ZET(J)**PN(N1)) BIJ = (BEZET**PWR(N1))*CHIJ(J) BIJ1 = 0.5E0*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(IPT,J) = PF(J)*(PWR2(N1)/ZET(J)- A (1.+2.*PN(N1))*BIJ1/BEZET) ELSE BIJ = 1.E0 PF(J) = 0.E0 PF2(J) = 0.E0 PF3(IPT,J) = 0.E0 ENDIF FRST(J) = FCD1*ELM1A - BIJ*FCD2*ELM2B FRST2(J) = FCD1P*ELM1A -BIJ*FCD2P*ELM2B FF1(IPT,J) = FRST(J)+PF(J)*SZET1(J) FF11(IPT,J) = FRST2(J)+2.E0*PF2(J)*SZET1(J) A +PF3(IPT,J)*SZET1(J)**2+PF(J)*SZET11(J) C CALCULATE POTENTIAL ENERGY PER PARTICLE VIJR = FC(J)*ELM1A VIJA = -FC(J)*BIJ*ELM2B VIJ = VIJR+VIJA EP(IPT) = EP(IPT) + VIJ/4.E0 EP(JPT) = EP(JPT) + VIJ/4.E0 C THREE BODY FORCES AND THREE BODY STRESSES DO 840 K=1,NEG KPT = LIST2(IPT,K) IF (K.EQ.J) GOTO 840 FF2(IPT,J,K) = PF(J)*ZET2(IPT,J,K) FF3(IPT,J,K) = PF(J)*ZET3(IPT,J,K) FF12(IPT,J,K) = PF2(J)*ZET2(IPT,J,K)+PF3(IPT,J)*SZET1(J) A *ZET2(IPT,J,K)+PF(J)*ZET12(J,K) FF13(IPT,J,K) = PF2(J)*ZET3(IPT,J,K)+PF(J)*ZET13(J,K) A +PF3(IPT,J)*SZET1(J)*ZET3(IPT,J,K) FF22(IPT,J,K) = PF3(IPT,J)*ZET2(IPT,J,K)**2+PF(J)*ZET22(J,K) FF23(IPT,J,K) = PF3(IPT,J)*ZET2(IPT,J,K)*ZET3(IPT,J,K) A +PF(J)*ZET23(J,K) FF33(IPT,J,K) = PF3(IPT,J)*ZET3(IPT,J,K)**2+PF(J)*ZET33(J,K) C CALC FORCE ON PARTICLE K FKX = -(FF2(IPT,J,K)*RX(IPT,K)/R(IPT,K)+FF3(IPT,J,K) A *(-PRJKX(J,K)))/2.E0 FKY = -(FF2(IPT,J,K)*RY(IPT,K)/R(IPT,K)+FF3(IPT,J,K) A *(-PRJKY(J,K)))/2.E0 FKZ = -(FF2(IPT,J,K)*RZ(IPT,K)/R(IPT,K)+FF3(IPT,J,K) A *(-PRJKZ(J,K)))/2.E0 FX(KPT) = FX(KPT) + FKX FY(KPT) = FY(KPT) + FKY FZ(KPT) = FZ(KPT) + FKZ 840 CONTINUE C CALC FORCE ON PARTICLE I FIX = ((FRST(J)+PF(J)*SZET1(J))*RXT(J)+PF(J)*SFTM2X(J))/2.E0 FIY = ((FRST(J)+PF(J)*SZET1(J))*RYT(J)+PF(J)*SFTM2Y(J))/2.E0 FIZ = ((FRST(J)+PF(J)*SZET1(J))*RZT(J)+PF(J)*SFTM2Z(J))/2.E0 FX(IPT) = FX(IPT)+FIX FY(IPT) = FY(IPT)+FIY FZ(IPT) = FZ(IPT)+FIZ C CALC FORCE ON PARTICLE J FJX=-((FRST(J)+PF(J)*SZET1(J))*RXT(J)+SFTM3X(J)*PF(J))*0.5E0 FJY=-((FRST(J)+PF(J)*SZET1(J))*RYT(J)+SFTM3Y(J)*PF(J))*0.5E0 FJZ=-((FRST(J)+PF(J)*SZET1(J))*RZT(J)+SFTM3Z(J)*PF(J))*0.5E0 FX(JPT) = FX(JPT)+FJX FY(JPT) = FY(JPT)+FJY FZ(JPT) = FZ(JPT)+FJZ 704 CONTINUE 500 CONTINUE POTE = 0.E0 DO 400 I =1,N POTE = POTE + EP(I) 400 CONTINUE c print *, pote/n,n c print *,' ' C CALCULATE THE DYNAMICAL MATRIX C THE NOTATION FOLLOWS THE ANALYTIC DERIVATION BY L.J.PORTER C I.E, 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 zero the dynamical matrix do 17 i = 1,3*n idx(i) = 0 17 continue C MEIJIE DEFINED XIJ,YIJ,ZIJ AS XJ - XI, ETC C NOTE ALSO THAT BELOW, I HAVE QUANTITIES LIKE C RJNX = RX(I,N) - RX(I,J) BECAUSE AFTER THE C LOOPS BELOW, RX(I,J) WILL = XI - XJ. do 21 ipt= 1,n do 22 j = 1,INDEX2(ipt) RX(ipt,j) = -RX(ipt,j) RY(ipt,j) = -RY(ipt,j) Rz(ipt,j) = -Rz(ipt,j) 22 continue 21 continue do 1002 ipt = 1,n c find the possible neighbours that could have non-zero c dynamical matrix element. ipw = 0 do 257 k=1,index2(ipt) ipw = ipw+1 list3(ipw)=list2(ipt,k) jpt = list2(ipt,k) do 257 kk=1,index2(jpt) ipw = ipw+1 list3(ipw)=list2(jpt,kk) 257 continue do 45 i=1,ipw if (list3(i).eq.0) goto 45 do 48 j=i+1,ipw if (list3(j).eq.list3(i)) list3(j) = 0 48 continue 45 continue do 1100 jkk = 1,ipw jpt = list3(jkk) if (jpt.eq.0.or.jpt.eq.ipt) goto 1100 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 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 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 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 negij = 0 do 654 iw = 1,index2(ipt) 654 if (list2(ipt,iw).eq.jpt) negij = iw negji = 0 do 659 iw = 1,index2(jpt) 659 if (list2(jpt,iw).eq.ipt) negji = iw if (negij.eq.0) goto 1800 rij = r(ipt,negij) rijx = RX(ipt,negij) rijy = RY(ipt,negij) rijz = RZ(ipt,negij) rij2 = r2(ipt,negij) rij3 = rij2*rij dijvij = ff1(ipt,negij) dijdijvij = ff11(ipt,negij) dijvji = ff1(jpt,negji) dijdijvji = ff11(jpt,negji) 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 do 1150 l = 1, INDEX2(ipt) if (l .eq. negij) goto 1150 ril = r(ipt,l) ril2 = ril*ril rilx = RX(ipt,l) rily = RY(ipt,l) rilz = RZ(ipt,l) rmujl = (RX(ipt,negij)*RX(ipt,l) + RY(ipt,negij)*RY(ipt,l) + a RZ(ipt,negij)*RZ(ipt,l))/(rij*ril) rjl2 = rij2 + ril2 -2.0*rmujl*rij*ril rjlx = RX(ipt,l)-RX(ipt,negij) rjly = RY(ipt,l)-RY(ipt,negij) rjlz = RZ(ipt,l)-RZ(ipt,negij) rjl = sqrt(rjl2) 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 do 1200 m = 1,INDEX2(ipt) if ((m .eq. negij).or.(m.eq.l)) goto 1200 rinx = Rx(ipt,m) riny = RY(ipt,m) rinz = RZ(ipt,m) rin = r(ipt,m) dindjlvij = pf3(ipt,negij)*zet2(ipt,negij,m)*zet3(ipt,negij,l) dildijvin = pf3(ipt,m)*zet2(ipt,m,negij)*zet2(ipt,m,l) dindjlvil = pf3(ipt,l)*zet2(ipt,l,m)*zet3(ipt,l,negij) 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) 1200 continue dijvil = ff2(ipt,l,negij) dijdijvil = ff22(ipt,l,negij) dijdjlvij = ff13(ipt,negij,l) dijdjlvil = ff23(ipt,l,negij) dijdilvil = ff12(ipt,l,negij) dijdilvij = ff12(ipt,negij,l) djldilvij = ff23(ipt,negij,l) djldilvil = ff13(ipt,l,negij) 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 1150 continue do 1250 l = 1,INDEX2(jpt) lpt = LIST2(jpt,l) if(lpt .eq. ipt) goto 1250 rjl = r(jpt,l) rjl2 = rjl*rjl rjlx = RX(jpt,l) rjly = RY(jpt,l) rjlz = RZ(jpt,l) rmuil = (RX(jpt,negji)*RX(jpt,l) + RY(jpt,negji)*RY(jpt,l) + a RZ(jpt,negji)*RZ(jpt,l))/(rij*rjl) ril2 = rij2 + rjl2 - 2.0*rmuil*rij*rjl ril = sqrt(ril2) rilx = RX(jpt,l)-RX(jpt,negji) rily = RY(jpt,l)-RY(jpt,negji) rilz = RZ(jpt,l)-RZ(jpt,negji) 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 do 1300 m = 1, INDEX2(jpt) npt = LIST2(jpt,m) if ((npt.eq.ipt) .or.(m.eq.l)) goto 1300 rjn = r(jpt,m) rjn2 = rjn*rjn rmuin = (RX(jpt,negji)*RX(jpt,m) + RY(jpt,negji)*RY(jpt,m) + a RZ(jpt,negji)*RZ(jpt,m))/(rij*rjn) rin2 = rij2 + rjn2 - 2.0*rmuin*rij*rjn rin = sqrt(rin2) rinx = RX(jpt,m)-RX(jpt,negji) riny = RY(jpt,m)-RY(jpt,negji) rinz = RZ(jpt,m)-RZ(jpt,negji) djldijvjn = pf3(jpt,m)*zet2(jpt,m,negji)*zet2(jpt,m,l) djldinvji = pf3(jpt,negji)*zet2(jpt,negji,l)*zet3(jpt,negji,m) djldinvjn = pf3(jpt,m)*zet2(jpt,m,l)*zet3(jpt,m,negji) 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) 1300 continue dijvjl = ff2(jpt,l,negji) djldijvji = ff12(jpt,negji,l) dijdijvjl = ff22(jpt,l,negji) dildjlvjl = ff13(jpt,l,negji) dijdjlvjl = ff12(jpt,l,negji) dildjlvji = ff23(jpt,negji,l) dijdilvji = ff13(jpt,negji,l) dijdilvjl = ff23(jpt,l,negji) 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 1250 continue 1800 continue do 1700 l = 1,INDEX2(ipt) lpt = LIST2(ipt,l) if (lpt .eq. jpt) goto 1700 negjl = 0 do 644 iw = 1,index2(jpt) 644 if (list2(jpt,iw).eq.lpt) negjl = iw if (negjl .eq. 0) goto 1700 negli = 0 do 634 iw = 1,index2(lpt) 634 if (list2(lpt,iw).eq.ipt) negli = iw neglj = 0 do 624 iw = 1,index2(lpt) 624 if (list2(lpt,iw).eq.jpt) neglj = iw ril = r(ipt,l) rilx = RX(ipt,l) rily = RY(ipt,l) rilz = RZ(ipt,l) ril2 = ril*ril rlj = r(lpt,neglj) rlj2 = rlj*rlj rljx = RX(lpt,neglj) rljy = RY(lpt,neglj) rljz = RZ(lpt,neglj) 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. if (negij .eq. 0) then rmuij=(RX(lpt,negli)*RX(lpt,neglj)+RY(lpt,negli)*RY(lpt,neglj) a +RZ(lpt,negli)*RZ(lpt,neglj))/(ril*rlj) rij2 = ril2+rlj2-2.0*rmuij*ril*rlj rij = sqrt(rij2) rij3 = rij2*rij rijx = RX(lpt,neglj)-RX(lpt,negli) rijy = RY(lpt,neglj)-RY(lpt,negli) rijz = RZ(lpt,neglj)-RZ(lpt,negli) endif 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 do 1750 m = 1,INDEX2(lpt) npt = LIST2(lpt,m) if((npt.eq.ipt).or.(npt.eq.jpt)) goto 1750 rln = r(lpt,m) rln2 = rln*rln rmuin = (RX(lpt,negli)*RX(lpt,m) + RY(lpt,negli)*RY(lpt,m) + a RZ(lpt,negli)*RZ(lpt,m))/(ril*rln) rin2 = ril2 + rln2 - 2.0*rmuin*ril*rln rin = sqrt(rin2) rinx = RX(lpt,m)-RX(lpt,negli) riny = RY(lpt,m)-RY(lpt,negli) rinz = RZ(lpt,m)-RZ(lpt,negli) rmujn = (RX(lpt,neglj)*RX(lpt,m) + RY(lpt,neglj)*RY(lpt,m) + a RZ(lpt,neglj)*RZ(lpt,m))/(rlj*rln) rjn2 = rln2 + rlj2 - 2.0*rmujn*rln*rlj rjn = sqrt(rjn2) rjnx = RX(lpt,m)-RX(lpt,neglj) rjny = RY(lpt,m)-RY(lpt,neglj) rjnz = RZ(lpt,m)-RZ(lpt,neglj) dindjlvli = pf3(lpt,negli)*zet2(lpt,negli,neglj) a *zet3(lpt,negli,m) dindjlvln = pf3(lpt,m)*zet2(lpt,m,neglj)*zet3(lpt,m,negli) dildjnvlj= pf3(lpt,neglj)*zet2(lpt,neglj,negli) a *zet3(lpt,neglj,m) dildjnvln = pf3(lpt,m)*zet2(lpt,m,negli)*zet3(lpt,m,neglj) dindijvli = pf3(lpt,negli)*zet3(lpt,negli,m) a *zet3(lpt,negli,neglj) dildjlvln = pf3(lpt,m)*zet2(lpt,m,negli)*zet2(lpt,m,neglj) dijdjnvlj = pf3(lpt,neglj)*zet3(lpt,neglj,negli) a *zet3(lpt,neglj,m) dindjnvln = pf3(lpt,m)*zet3(lpt,m,negli)*zet3(lpt,m,neglj) 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 1750 continue dijvli = ff3(lpt,negli,neglj) dijvlj = ff3(lpt,neglj,negli) dildjlvli = ff12(lpt,negli,neglj) dildjlvlj = ff12(lpt,neglj,negli) dildijvli = ff13(lpt,negli,neglj) dildijvlj = ff23(lpt,neglj,negli) dijdjlvli = ff23(lpt,negli,neglj) dijdjlvlj = ff13(lpt,neglj,negli) dijdijvli = ff33(lpt,negli,neglj) dijdijvlj = ff33(lpt,neglj,negli) 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 smij = sqrt(pmass(ipt)*pmass(jpt)) ibb = 3*(ipt-1)+1 jbb = 3*(jpt-1)+1 DXXij = 0.50*(firstxx + lloopxx + lloop2xx + lloop5xx)/smij if (abs(DXXij).gt.TOL) then idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb hess(idp,ibb) = DXXij endif DXYij = 0.50*(firstxy + lloopxy + lloop2xy + lloop5xy)/smij if (abs(DXYij).gt.TOL) then idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+1 hess(idp,ibb) = DXYij endif DXZij = 0.50*(firstxz + lloopxz + lloop2xz + lloop5xz)/smij if (abs(DXZij).gt.TOL) then idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+2 hess(idp,ibb) = DXZij endif DYXij = 0.50*(firstyx + lloopyx + lloop2yx + lloop5yx)/smij if (abs(DYXij).gt.TOL) then ibb = ibb+1 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb hess(idp,ibb) = DYXij ibb = ibb-1 endif DYYij = 0.50*(firstyy + lloopyy + lloop2yy + lloop5yy)/smij if (abs(DYYij).gt.TOL) then ibb = ibb+1 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+1 hess(idp,ibb) = DYYij ibb = ibb-1 endif DYZij = 0.50*(firstyz + lloopyz + lloop2yz + lloop5yz)/smij if (abs(DYZij).gt.TOL) then ibb = ibb+1 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+2 hess(idp,ibb) = DYZij ibb = ibb-1 endif DZXij = 0.50*(firstzx + lloopzx + lloop2zx + lloop5zx)/smij if (abs(DZXij).gt.TOL) then ibb = ibb+2 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb hess(idp,ibb) = DZXij ibb = ibb-2 endif DZYij = 0.50*(firstzy + lloopzy + lloop2zy + lloop5zy)/smij if (abs(DZYij).gt.TOL) then ibb = ibb+2 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+1 hess(idp,ibb) = DZYij ibb = ibb-2 endif DZZij = 0.50*(firstzz + lloopzz + lloop2zz + lloop5zz)/smij if (abs(DZZij).gt.TOL) then ibb = ibb+2 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+2 hess(idp,ibb) = DZZij ibb = ibb-2 endif 1100 continue C THIS PART CALCULATES DXX(Ipt,I) DXY(Ipt,I) ETC, IE, J = I 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 do 1500 k = 1,INDEX2(ipt) kpt = LIST2(ipt,k) negki = 0 do 604 iw = 1,index2(kpt) 604 if (list2(kpt,iw).eq.ipt) negki = iw rik = r(ipt,k) rik2 = rik*rik rik3 = rik2*rik rikx = RX(ipt,k) riky = RY(ipt,k) rikz = Rz(ipt,k) dikvik = ff1(ipt,k) dikvki = ff1(kpt,negki) dikdikvik = ff11(ipt,k) dikdikvki = ff11(kpt,negki) 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) 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 do 1550 l = 1, INDEX2(ipt) if (l.eq.k) goto 1550 ril = r(ipt,l) rilx = RX(ipt,l) rily = RY(ipt,l) rilz = RZ(ipt,l) nloop2 = 0.0 do 1580 m = 1, INDEX2(ipt) if ((m.eq.k).or.(m.eq.l)) goto 1580 dildikvin = pf3(ipt,m)*zet2(ipt,m,l)*zet2(ipt,m,k) nloop2 = nloop2 + dildikvin 1580 continue dikvil = ff2(ipt,l,k) dikdikvil = ff22(ipt,l,k) dikdilvik = ff12(ipt,k,l) dikdilvil = ff12(ipt,l,k) 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) 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 do 1600 l = 1, INDEX2(kpt) lpt = LIST2(kpt,l) if(lpt.eq.ipt) goto 1600 rkl = r(kpt,l) rkl2 = rkl*rkl rklx = RX(kpt,l) rkly = RY(kpt,l) rklz = RZ(kpt,l) rmuil = (RX(kpt,negki)*RX(kpt,l) + RY(kpt,negki)*RY(kpt,l) + a RZ(kpt,negki)*RZ(kpt,l))/(rik*rkl) ril2 = rik2 + rkl2 -2.0*rmuil*rik*rkl ril = sqrt(ril2) ril3 = ril2*ril rilx = -1.0*(RX(kpt,negki) - RX(kpt,l)) rily = -1.0*(RY(kpt,negki) - RY(kpt,l)) rilz = -1.0*(RZ(kpt,negki) - RZ(kpt,l)) 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 do 1650 m = 1, INDEX2(kpt) npt = LIST2(kpt,m) if ((npt.eq.ipt) .or.(m.eq.l)) goto 1650 rkn = r(kpt,m) rkn2 = rkn*rkn rmuin = (RX(kpt,negki)*RX(kpt,m) + RY(kpt,negki)*RY(kpt,m) + a RZ(kpt,negki)*RZ(kpt,m))/(rik*rkn) rin2 = rik2 + rkn2 - 2.0*rmuin*rik*rkn rin = sqrt(rin2) rinx = RX(kpt,m)-RX(kpt,negki) riny = RY(kpt,m)-RY(kpt,negki) rinz = RZ(kpt,m)-RZ(kpt,negki) dildinvki = pf3(kpt,negki)*zet3(kpt,negki,l)*zet3(kpt,negki,m) 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 1650 continue dikvkl = ff2(kpt,l,negki) dilvki = ff3(kpt,negki,l) dilvkl = ff3(kpt,l,negki) dildilvkl = ff33(kpt,l,negki) dildikvki = ff13(kpt,negki,l) dildikvkl = ff23(kpt,l,negki) dildilvki = ff33(kpt,negki,l) dikdikvkl = ff22(kpt,l,negki) 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 1600 continue 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 1500 continue smij = pmass(ipt) ibb = 3*(ipt-1)+1 jbb = ibb DXXij = 0.50*sumkxx/smij DYYij = 0.50*sumkyy/smij DZZij = 0.50*sumkzz/smij DXYij = 0.50*sumkxy/smij DXZij = 0.50*sumkxz/smij DYXij = 0.50*sumkyx/smij DYZij = 0.50*sumkyz/smij DZXij = 0.50*sumkzx/smij DZYij = 0.50*sumkzy/smij if (abs(DXXij).gt.TOL) then idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb hess(idp,ibb) = DXXij endif if (abs(DXYij).gt.TOL) then idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+1 hess(idp,ibb) = DXYij endif if (abs(DXZij).gt.TOL) then idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+2 hess(idp,ibb) = DXZij endif if (abs(DYXij).gt.TOL) then ibb = ibb+1 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb hess(idp,ibb) = DYXij ibb = ibb-1 endif if (abs(DYYij).gt.TOL) then ibb = ibb+1 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+1 hess(idp,ibb) = DYYij ibb = ibb-1 endif if (abs(DYZij).gt.TOL) then ibb = ibb+1 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+2 hess(idp,ibb) = DYZij ibb = ibb-1 endif if (abs(DZXij).gt.TOL) then ibb = ibb+2 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb hess(idp,ibb) = DZXij ibb = ibb-2 endif if (abs(DZYij).gt.TOL) then ibb = ibb+2 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+1 hess(idp,ibb) = DZYij ibb = ibb-2 endif if (abs(DZZij).gt.TOL) then ibb = ibb+2 idx(ibb) = idx(ibb)+1 idp = idx(ibb) id(idp,ibb) = jbb+2 hess(idp,ibb) = DZZij ibb = ibb-2 endif 1002 continue c read *, i,j c i = 1 c j = 2 c print *,idx(i),idx(j) c do 456 k=1,idx(i) c print *,hess(k,i) c456 continue return end FUNCTION ran1(idum) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL ran1,AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum ran1=min(AM*iy,RNMX) return END