C --------------------------------------------------------- C PROGRAM Tersoff v1.0: C COHESIVE ENERGY CURVES OF BINARY SI-C SYSTEMS C ENERGY UNIT = eV = 1.602D-19 JOULE C LENGTH UNIT = ANGSTROM C Aug. 14, 1998 Ju Li (MIT) C --------------------------------------------------------- PROGRAM Tersoff IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MMA=4,MML=2,MMH=(2*MML+1)**3,NM=MMH*MMA,NX=60) PARAMETER (PI=3.14159265358979323846D0,EPS=1E-4) COMMON /DOUBLE/ TAU(MMA,3),S(NM,3),H(3,3),HI(3,3),EP(NM), A AA,DET,TSF,CHI_SI_C,PRSIC,PSSIC,RCUT,CUBE0, A 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) COMMON /INT/ NPA,NTYPE(NM),LIST(NM,NX),IDX(NM),LIST2(NX) COMMON /CHAR/ BUF,NAME CHARACTER*70 BUF,NAME DATA LP,LP_OUT /25,26/ PA(1) = 1.8308D3 PB(1) = 4.7118D2 PL1(1) = 2.4799 PL2(1) = 1.7322 PN(1) = 7.8734D-1 PBT(1) = (1.1D-6)**PN(1) PC(1) = (1.0039D5)**2. PD(1) = (1.6217D1)**2. PH(1) = -5.9825D-1 PR(1) = 2.7 PS(1) = 3.0 C2D2(1) = 1.+PC(1)/PD(1) PWR(1) = -1./(2.*PN(1)) PWR1(1) = PWR(1)-1. PWR2(1) = PN(1)-1. PA(2) = 1.3936D3 PB(2) = 3.4674D2 PL1(2) = 3.4879 PL2(2) = 2.2119 PN(2) = 7.2751D-1 PBT(2) = (1.5724D-7)**PN(2) PC(2) = (3.8049D4)**2. PD(2) = (4.3484)**2. PH(2) = -5.7058D-1 PR(2) = 1.8 PS(2) = 2.1 C2D2(2) = 1.+PC(2)/PD(2) PWR(2) = -1./(2.*PN(2)) PWR1(2) = PWR(2)-1. PWR2(2) = PN(2)-1. CHI_SI_C = 0.9776 C Meijie's parameter: PRSIC = 2.36D0 PSSIC = 2.56D0 RCUT = MAX(PS(1),PS(2)) READ (*,'("Number of Object Groups: ",I)') NOBJ DO 10 IOBJ = 1,NOBJ READ (*,'(/,"Object Group Name: ",A70)') NAME READ (*,'("Structure Description: ",A70)') BUF PRINT *,'Loading structure from ', BUF(1:INDEX(BUF,' ')-1) OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED',FILE=BUF) C Read file header: "Name:", NAME, "Lattice:" READ (LP,'(A70,/,A70,/,A70)') BUF,BUF,BUF READ (LP,*) ((H(I,J),J=1,3),I=1,3) C Invert the matrix: D11 = H(2,2)*H(3,3)-H(2,3)*H(3,2) D22 = H(3,3)*H(1,1)-H(3,1)*H(1,3) D33 = H(1,1)*H(2,2)-H(1,2)*H(2,1) D12 = H(2,3)*H(3,1)-H(2,1)*H(3,3) D23 = H(3,1)*H(1,2)-H(3,2)*H(1,1) D31 = H(1,2)*H(2,3)-H(1,3)*H(2,2) D13 = H(2,1)*H(3,2)-H(3,1)*H(2,2) D21 = H(3,2)*H(1,3)-H(1,2)*H(3,3) D32 = H(1,3)*H(2,1)-H(2,3)*H(1,1) DET = H(1,1)*D11+H(1,2)*D12+H(1,3)*D13 HI(1,1) = D11/DET HI(2,2) = D22/DET HI(3,3) = D33/DET HI(1,2) = D21/DET HI(2,3) = D32/DET HI(3,1) = D13/DET HI(2,1) = D12/DET HI(3,2) = D23/DET HI(1,3) = D31/DET C "Number of atoms in the unit cell:" READ (LP,'(A70,/,I,/,A70)') BUF,NPA,BUF C "Type X Y Z" DO N = 1,NPA READ (LP,*) NTYPE(N),(TAU(N,I),I=1,3) S(N,1) = HI(1,1)*TAU(N,1)+HI(2,1)*TAU(N,2)+HI(3,1)*TAU(N,3) S(N,2) = HI(1,2)*TAU(N,1)+HI(2,2)*TAU(N,2)+HI(3,2)*TAU(N,3) S(N,3) = HI(1,3)*TAU(N,1)+HI(2,3)*TAU(N,2)+HI(3,3)*TAU(N,3) ENDDO CLOSE (LP) IP = 0 DO II = 0,2*MML DO JJ = 0,2*MML DO KK = 0,2*MML DO N = 1,NPA IP = IP+1 S(IP,1) = II+S(N,1) S(IP,2) = JJ+S(N,2) S(IP,3) = KK+S(N,3) NTYPE(IP) = NTYPE(N) ENDDO ENDDO ENDDO ENDDO READ (*,'("CPS database: ",A70)') BUF PRINT *,'Loading CPS data from ', BUF(1:INDEX(BUF,' ')-1) PRINT *,'Calculating cohesive energy for ', A NAME(1:INDEX(NAME,' ')-1) PRINT *,' ' PRINT *,'Length(A) Volume(A^3) CPS(eV) Tersoff(eV)' OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED',FILE=BUF) READ (LP, *) NDATA, NDATA, NDATA OPEN (UNIT=LP_OUT, STATUS='UNKNOWN', FORM='FORMATTED', A FILE=NAME(1:INDEX(NAME,' ')-1)//'.out') CURV_ERROR = 0.D0 ABS_CURV_ERROR = 0.D0 DO 20 I = 1,NDATA READ (LP, *) AA, VOLUME, CPS CALL TERSOFF_VALUE WRITE (*,'(F9.6,3(F11.6))') AA,ABS(DET)*AA**3/NPA,CPS,TSF WRITE (LP_OUT,'(F9.6,3(F11.6))') AA,ABS(DET)*AA**3/NPA,CPS,TSF IF (I.EQ.1) THEN TSF_GROUND = TSF CPS_GROUND = CPS ELSE CC = (TSF-TSF_GROUND-CPS+CPS_GROUND)/(CPS-CPS_GROUND)/(NDATA-1) CURV_ERROR = CURV_ERROR + CC ABS_CURV_ERROR = ABS_CURV_ERROR + ABS(CC) ENDIF 20 CONTINUE CLOSE (LP_OUT) CLOSE (LP) WRITE (*,'(" Relative ground error = ",f5.1, A "%, curv error = ",f5.1,"% (",f4.1,"%)",/)') A -(TSF_GROUND-CPS_GROUND)/CPS_GROUND*100, CURV_ERROR*100, A ABS_CURV_ERROR*100 10 CONTINUE STOP END SUBROUTINE TERSOFF_VALUE IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MMA=4,MML=2,MMH=(2*MML+1)**3,NM=MMH*MMA,NX=60) PARAMETER (PI=3.14159265358979323846D0) COMMON /DOUBLE/ TAU(MMA,3),S(NM,3),H(3,3),HI(3,3),EP(NM), A AA,DET,TSF,CHI_SI_C,PRSIC,PSSIC,RCUT,CUBE0, A 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) COMMON /INT/ NPA,NTYPE(NM),LIST(NM,NX),IDX(NM),LIST2(NX) DIMENSION FC(NM),RX(NX),RY(NX),RZ(NX),RXT(NX), A RYT(NX),RZT(NX),R(NX),ZET(NX),R2(NX) DO I=1,NPA*MMH EP(I) = 0. ENDDO C --------------------------------------------------------- C ALWAYS UPDATE THE VERLET LIST DO IPT = 1,NPA*MMH IDX(IPT) = 0 ENDDO DO 170 IPT = 1,NPA*MMH N1 = NTYPE(IPT) SXI = S(IPT,1) SYI = S(IPT,2) SZI = S(IPT,3) DO 270 JPT = IPT+1,NPA*MMH N2 = NTYPE(JPT) SXJI = S(JPT,1) - SXI SYJI = S(JPT,2) - SYI SZJI = S(JPT,3) - SZI IF(SXJI.GT.MML+0.5D0) SXJI=SXJI-2*MML-1.D0 IF(SXJI.LT.-MML-0.5D0) SXJI=SXJI+2*MML+1.D0 IF(SYJI.GT.MML+0.5D0) SYJI=SYJI-2*MML-1.D0 IF(SYJI.LT.-MML-0.50) SYJI=SYJI+2*MML+1.D0 IF(SZJI.GT.MML+0.5D0) SZJI=SZJI-2*MML-1.D0 IF(SZJI.LT.-MML-0.5D0) SZJI=SZJI+2*MML+1.D0 RXX = AA * (H(1,1)*SXJI + H(2,1)*SYJI + H(3,1)*SZJI) RYY = AA * (H(1,2)*SXJI + H(2,2)*SYJI + H(3,2)*SZJI) RZZ = AA * (H(1,3)*SXJI + H(2,3)*SYJI + H(3,3)*SZJI) RR = RXX*RXX+RYY*RYY+RZZ*RZZ IF (RR.LT.RCUT**2) THEN IF (IDX(IPT).EQ.NX.OR.IDX(JPT).EQ.NX) THEN PRINT *, 'tersoff_value(): NX out of bound.' STOP ENDIF IF (IDX(IPT).LT.NX.AND.IDX(JPT).LT.NX) THEN IDX(IPT)=IDX(IPT)+1 LIST(IPT,IDX(IPT))=JPT IDX(JPT)=IDX(JPT)+1 LIST(JPT,IDX(JPT))=IPT ENDIF ENDIF 270 CONTINUE 170 CONTINUE C --------------------------------------------------------- DO 500 IPT = 1,NPA*MMH N1 = NTYPE(IPT) SXI = S(IPT,1) SYI = S(IPT,2) SZI = S(IPT,3) NEG = 0 DO 600 J = 1,IDX(IPT) JPT = LIST(IPT,J) N2 = NTYPE(JPT) RSMAX = PS(N1)*PS(N2) SXJI = S(JPT,1) - SXI SYJI = S(JPT,2) - SYI SZJI = S(JPT,3) - SZI IF(SXJI.GT.MML+0.5D0) SXJI=SXJI-2*MML-1.D0 IF(SXJI.LT.-MML-0.5D0) SXJI=SXJI+2*MML+1.D0 IF(SYJI.GT.MML+0.5D0) SYJI=SYJI-2*MML-1.D0 IF(SYJI.LT.-MML-0.50) SYJI=SYJI+2*MML+1.D0 IF(SZJI.GT.MML+0.5D0) SZJI=SZJI-2*MML-1.D0 IF(SZJI.LT.-MML-0.5D0) SZJI=SZJI+2*MML+1.D0 RXX = AA * (H(1,1)*SXJI + H(2,1)*SYJI + H(3,1)*SZJI) RYY = AA * (H(1,2)*SXJI + H(2,2)*SYJI + H(3,2)*SZJI) RZZ = AA * (H(1,3)*SXJI + H(2,3)*SYJI + H(3,3)*SZJI) RR = RXX*RXX+RYY*RYY+RZZ*RZZ IF(RR.GT.RSMAX) GOTO 600 NEG = NEG + 1 R1 = SQRT(RR) RX(NEG) = RXX RY(NEG) = RYY RZ(NEG) = RZZ RXT(NEG) = RXX/R1 RYT(NEG) = RYY/R1 RZT(NEG) = RZZ/R1 R2(NEG) = RR R(NEG) = R1 LIST2(NEG) = JPT 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)) RLM1(NEG) = (PL1(N1)+PL1(N2))*0.5 RLM2(NEG) = (PL2(N1)+PL2(N2))*0.5 IF(N1.EQ.N2) THEN CHIJ(NEG) = 1. ELSE CHIJ(NEG) = CHI_SI_C ENDIF IF(R1.LT.CR(NEG)) THEN FC(NEG) = 1. ELSE FC(NEG)=0.5*(1.+COS(PI*(R1-CR(NEG)) A /(CS(NEG)-CR(NEG)))) ENDIF 600 CONTINUE DO 700 J = 1,NEG JPT = LIST2(J) N22 = NTYPE(JPT) ZET(J) = 0. DO 800 K = 1,NEG IF(K.EQ.J) GOTO 800 KP = LIST2(K) RMU = (RX(J)*RX(K)+RY(J)*RY(K)+RZ(J)*RZ(K))/R(J)/R(K) RMU1 = RMU-PH(N1) RMU2 = RMU1*RMU1 RMU3 = PD(N1)+RMU2 GMU = C2D2(N1)-PC(N1)/RMU3 ZET(J) = ZET(J)+FC(K)*GMU 800 CONTINUE 700 CONTINUE DO J=1,NEG JPT = LIST2(J) N22 = NTYPE(JPT) ELM2 = EXP(-RLM2(J)*R(J)) ELM2B = ELM2*CB(J) ELM1 = EXP(-RLM1(J)*R(J)) ELM1A = ELM1*CA(J) IF(ZET(J).NE.0.) THEN BEZET = 1.+ PBT(N1)*(ZET(J)**PN(N1)) BIJ = (BEZET**PWR(N1))*CHIJ(J) ELSE BIJ = 1. ENDIF C CALCULATE POTENTIAL ENERGY PER PARTICLE VIJR = FC(J)*ELM1A VIJA = -FC(J)*BIJ*ELM2B VIJ = VIJR+VIJA EP(IPT) = EP(IPT) + VIJ/4. EP(JPT) = EP(JPT) + VIJ/4. ENDDO 500 CONTINUE POTE = 0. DO I = 1,NPA*MMH POTE = POTE+EP(I) ENDDO TSF = POTE/(NPA*MMH) RETURN END