C -------------------------------------------------------------------- C PROGRAM DEFECT: C VERSION 3.1 (JULY.11,1997): C 1. TERSOFF POTENTIAL FOR BINARY SI/C SYSTEMS. C 2. USE CONJUGATE GRADIENT MINIMIZER. C 3. RELAX CONFIGURATION AT CONSTANT H OR CONSTANT STRESS. C C JU LI C -------------------------------------------------------------------- C USE THE FOLLOWING UNITS: C ENERGY UNIT = eV = 1.602D-19 JOULE C MASS UNIT = AMU = 1.66054D-27 KG C LENGTH UNIT = ANGSTROM C STRESS UNIT = 1.602192 MBAR C -------------------------------------------------------------------- PROGRAM DEFECT INCLUDE 'defect.fh' DATA H0,EP0,CUBE0 /9*0.0,-6.386496963D0,4.28533/ DOUBLE PRECISION XX(MOP),ETA(3,3) DATA ETA /9*0.d0/ CALL INITIALIZE_CONSTANTS C ========================================================================== C Read instructions from input file, such as "con_defect": READ (*,'(A70)') buf C "To relax 0) one pt-defect 1) all pt-defects 2) user-defined configuration:" READ *, OPTION READ (*,'(/,A70)') buf C "CONST_H or CONST_S mode:" READ (*,'(A70)') BUF IF (BUF(1:7).EQ."CONST_H") THEN CONST_H = .TRUE. CONST_S = .FALSE. ELSEIF (BUF(1:7).EQ."CONST_S") THEN CONST_H = .FALSE. CONST_S = .TRUE. ELSE STOP 'Must input a running mode: "CONST_H" or "CONST_S".' ENDIF READ (*,'(/,A70)') buf C "Default number of atoms in corresponding perfect(3C) state:" READ *, NOLD IF(NOLD.GT.NM) STOP 'Error(defect.fh): too many atoms, revise NM.' READ (*,'(/,A70)') buf C "Default lattice parameter(3C) in angstrom:" READ *, CUBE NC = NINT((NOLD/8.)**(1./3.)) H0(1,1) = CUBE*NC H0(2,2) = CUBE*NC H0(3,3) = CUBE*NC CALL NEW_H (ETA) READ (*,'(/,A70)') buf C "Default stress (11,22,33,12,13,23), in MBar:" READ *,SOUT(1,1),SOUT(2,2),SOUT(3,3),SOUT(1,2),SOUT(1,3),SOUT(2,3) SOUT(1,1) = SOUT(1,1) / USTRESS_IN_MBAR SOUT(2,2) = SOUT(2,2) / USTRESS_IN_MBAR SOUT(3,3) = SOUT(3,3) / USTRESS_IN_MBAR SOUT(1,2) = SOUT(1,2) / USTRESS_IN_MBAR SOUT(1,3) = SOUT(1,3) / USTRESS_IN_MBAR SOUT(2,3) = SOUT(2,3) / USTRESS_IN_MBAR SOUT(2,1) = SOUT(1,2) SOUT(3,1) = SOUT(1,3) SOUT(3,2) = SOUT(2,3) READ (*,'(/,A70)') buf C "Default name of the configurational file that could be read in:" READ (*,'(A70)') NAME_READ READ (*,'(/,A70)') buf C "Default name of user defined defect:" READ (*,'(A70)') NAME READ (*,'(/,A70,/,A70)') buf, buf C "Default index of the point defect that could be relaxed:" C 0.perfect 1.Vc 2.VSi 3.C_si 4.Si_c 5.CTc 6.CTsi 7.SiTsi 8.SiTc 9.antisite pair READ *, IW READ (*,'(/,A70)') buf C "Write on configurational file:" READ (*,'(A70)') NAME_WRITE READ (*,'(/,A70)') buf C "Give the initial configuration a perturbation:" READ *, PERTURBATION READ (*,'(/,A70)') buf C "Random number generator seed:" READ *, SEED READ (*,'(/,A70)') buf C "Minimization tolerance:" READ *, TOL READ (*,'(/,A70,/,A70)') buf, buf C Print out the position of atoms around particle i (-1=NO,0=default) within C R and save to XYZ file before (*.xyz.bef) and after optimization (*.xyz.opt): READ *, IATOM, RPRINT PRINT_ATOM = IATOM.NE.-1 USE_DEFAULT = IATOM.EQ.0 C ========================================================================== IF (OPTION.EQ.1) IW = 1 154 IF (OPTION.EQ.2) THEN OPEN (UNIT=LP_CONFIG, STATUS='UNKNOWN', FORM='UNFORMATTED', A FILE=NAME_READ) READ (LP_CONFIG) N IF (N.GT.NM) STOP 'Error(defect.fh): too many atoms, revise NM.' READ (LP_CONFIG) (NTYPE(I),I=1,N) READ (LP_CONFIG) ((H0(I,J),J=1,3),I=1,3) READ (LP_CONFIG) (SX(I),SY(I),SZ(I),I=1,N) C The velocity information is of no use. CLOSE (LP_CONFIG) CALL NEW_H (ETA) ELSE C -------------------------------------------------- C DIAMOND CUBIC STRUCTURE N = NOLD A = 1./FLOAT(NC) X(1) = 0. Y(1) = 0. Z(1) = 0. X(3) = 1./2. Y(3) = 1./2. Z(3) = 0. X(5) = 0. Y(5) = 1./2. Z(5) = 1./2. X(7) = 1./2. Y(7) = 0. Z(7) = 1./2. 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./8.)*A SY(IP+L) = (Y(L)+J-7./8.)*A SZ(IP+L) = (Z(L)+K-7./8.)*A NTYPE (IP+L) = 1 SX(IP+L+1) = (X(L)+I-5./8.)*A SY(IP+L+1) = (Y(L)+J-5./8.)*A SZ(IP+L+1) = (Z(L)+K-5./8.)*A NTYPE (IP+L+1) = 2 218 CONTINUE IP = IP+8 217 CONTINUE IF (CONST_H.AND.(.NOT.((OPTION.EQ.1).AND.(IW.GT.1)))) THEN C Calculate the perfect crystal energy at this volume for benchmark: CALL NEW_H (ETA) CALL TERSOFF_VALUE EPH = POTE/N ENDIF C Cohesive energy curve: C EMAX = 0.2 C MESH = 20 C DO I = -MESH,MESH C ETA(1,1) = EMAX/MESH*I C ETA(2,2) = ETA(1,1) C ETA(3,3) = ETA(1,1) C CALL NEW_H (ETA) C CALL TERSOFF_VALUE C CALL TERSOFF_FORCE C EPH = POTE/N C PRINT *,H(1,1)**3/N,EPH,STRESS(1,1) C ENDDO C STOP C ENDIF IF (IW.EQ.0) THEN c perfect crystal name = 'perfect crystal' eab = 0. if (use_default) iatom = 1 ELSEIF (IW.EQ.1) THEN c carbon vacancy sx(2) = sx(n) sy(2) = sy(n) sz(2) = sz(n) ntype(2) = ntype(n) n = n-1 name = 'Vc' c ab initio result: eab = 5.9 if (use_default) iatom = 1 ELSEIF (IW.EQ.2) THEN c silicon vacancy sx(1) = sx(n) sy(1) = sy(n) sz(1) = sz(n) ntype(1) = ntype(n) n = n-1 name = 'VSi' eab = 6.8 if (use_default) iatom = 2 ELSEIF (IW.EQ.3) THEN c carbon antisite ntype(1) = 2 name = 'C_si' eab = 1.1 if (use_default) iatom = 1 ELSEIF (IW.EQ.4) THEN c silicon antisite ntype(2) = 1 name = 'Si_c' eab = 7.3 if (use_default) iatom = 2 ELSEIF (IW.EQ.5) THEN c interstitial C at Tc site n = n + 1 ntype(n) = 2 sx(n) = sx(1)+A/2. sy(n) = sy(1)+A/2. sz(n) = sz(1)+A/2. name = 'CTc' eab = 11.0 if (use_default) iatom = n ELSEIF (IW.EQ.6) THEN c interstitial C at TSi site n = n + 1 ntype(n) = 2 sx(n) = sx(2)+A/2. sy(n) = sy(2)+A/2. sz(n) = sz(2)+A/2. name = 'CTsi' eab = 8.6 if (use_default) iatom = n ELSEIF (IW.EQ.7) THEN c interstitial Si at TSi site n = n + 1 ntype(n) = 1 sx(n) = sx(2)+A/2. sy(n) = sy(2)+A/2. sz(n) = sz(2)+A/2. name = 'SiTsi' eab = 15.0 if (use_default) iatom = n ELSEIF (IW.EQ.8) THEN c interstitial Si at Tc site n = n + 1 ntype(n) = 1 sx(n) = sx(1)+A/2. sy(n) = sy(1)+A/2. sz(n) = sz(1)+A/2. name = 'SiTc' eab = 14.7 if (use_default) iatom = n ELSEIF (IW.EQ.9) THEN c nearest neighbour antisite pair ntype(1) = 2 ntype(2) = 1 name = 'antisite pair' eab = 5.9 if (use_default) iatom = 2 ELSEIF (IW.EQ.10) THEN c di-vacancy sx(1) = sx(n-1) sy(1) = sy(n-1) sz(1) = sz(n-1) ntype(1) = ntype(n-1) sx(2) = sx(n) sy(2) = sy(n) sz(2) = sz(n) ntype(2) = ntype(n) n = n-2 name = 'divacancy' eab = 8.1 if (use_default) iatom = 3 ENDIF C -------------------------------------------------- ENDIF C ======================================================================== C Write header: IF (.NOT.((OPTION.EQ.1).AND.(IW.GT.1))) PRINT *, ' ' IF (OPTION.EQ.0) THEN print *, 'Relaxing pt-defect ',name(1:index(name,' ')-1) ELSEIF (OPTION.EQ.1) THEN if (iw.eq.1) print *, 'Relaxing all point defects in store ' ELSEIF (OPTION.EQ.2) THEN print *, 'Relaxing the configuration from ', A name_read(1:index(name_read,' ')-1) ELSE stop 'input valid option 1-3.' ENDIF if (option.ne.1) print *, 'of',N,' atoms in cell with PBC' IF (.NOT.((OPTION.EQ.1).AND.(IW.GT.1))) THEN if (const_h) then print *, 'under constant H (angstrom) ' write (*,422) ((H0(I,J),J=1,3),I=1,3) else print *, 'under constant stress (MBar) ' write (*,422) ((SOUT(I,J)*USTRESS_IN_MBAR,J=1,3),I=1,3) endif ENDIF 422 format(3('|',3(1x,f11.5,1x),'|',/)) if (option.ne.1) then print *, A 'and save the relaxed configuration to file ', A name_write(1:index(name_write,' ')-1) // '.' print *, ' ' endif IF (.NOT.((OPTION.EQ.1).AND.(IW.GT.1))) print *, A 'Perfect crystal energy/atom at zero pressure =',ep0,'eV.' IF (CONST_H.AND.(.NOT.((OPTION.EQ.1).AND.(IW.GT.1)))) A WRITE (6,408) EPH 408 FORMAT (' Perfect crystal energy/atom at this volume = ',F10.6, A ' eV.') IF (.NOT.((OPTION.EQ.1).AND.(IW.GT.1))) THEN print *, ' ' call tersoff_value print*,'Before relaxation the total (free) energy = ',pote,'eV,' print*,' (free) energy/atom = ',pote/n,'eV.' print*,' ' volume_before_relax = volume print*,'Before relaxation the total volume = ',volume,'A^3' print*,' ' ENDIF buf = name(1:index(name,' ')-1) // '.xyz.bef' if (IATOM.ne.-1) call print_atoms(buf) C ======================================================================== C A little perturbation: do i = 1,n call random sx(i) = sx(i) + perturbation*(seed-0.5)/((N/8.)**(1./3.)) call random sy(i) = sy(i) + perturbation*(seed-0.5)/((N/8.)**(1./3.)) call random sz(i) = sz(i) + perturbation*(seed-0.5)/((N/8.)**(1./3.)) enddo IP = 1 DO I=1,N XX(IP) = SX(I) XX(IP+1) = SY(I) XX(IP+2) = SZ(I) IP = IP+3 ENDDO IF (CONST_S) THEN XX(IP) = ETA(1,1) XX(IP+1) = ETA(2,2) XX(IP+2) = ETA(3,3) XX(IP+3) = ETA(1,2) XX(IP+4) = ETA(1,3) XX(IP+5) = ETA(2,3) IP = IP+6 ENDIF CALL FRPRMN (XX,IP-1,TOL,ITER,FMIN) IP = 1 DO I=1,N SX(I) = XX(IP) SY(I) = XX(IP+1) SZ(I) = XX(IP+2) IP = IP+3 ENDDO IF (CONST_S) THEN ETA(1,1) = XX(IP) ETA(2,2) = XX(IP+1) ETA(3,3) = XX(IP+2) ETA(1,2) = XX(IP+3) ETA(1,3) = XX(IP+4) ETA(2,3) = XX(IP+5) IP = IP+6 ETA(2,1) = ETA(1,2) ETA(3,1) = ETA(1,3) ETA(3,2) = ETA(2,3) ENDIF IF (OPTION.NE.1) THEN write(*,'(/," Minimization converged after",I5," iterations.",/)') A ITER write(*, '(" After relaxation the total (free) energy ="' // A ',F15.7," eV.",/)') FMIN if (option.ne.2.and.const_h) then write (*,542) fmin,n,eph,fmin-n*eph print *, A '** Potential energy reference is perfect crystal at this H. **' else write (*,542) fmin,n,ep0,fmin-n*ep0 print*, A '** Free energy reference is perfect crystal at zero stress. **' endif 542 format(' The defect formation energy =',F14.6,' -',I5,' *',F10.6, A /, ' =',F14.6,' eV.',/) print*,' ' print*,'After relaxation the total volume = ',volume,'A^3,' print*,'and has increased by', volume-volume_before_relax,'A^3.' print*,' ' ELSE if (iw.eq.1) write A (*,'(" Name Our formation energy(eV) Ab initio(eV)")') if (const_h) then write (*,438) name(1:index(name,' ')-1),fmin-n*eph,eab else write (*,438) name(1:index(name,' ')-1),fmin-n*ep0,eab endif 438 format(1x,a4,6x,f12.5,8x,f13.5) ENDIF if (option.eq.1.and.iw.eq.10) then if (const_h) then write (*, '(/, "** Potential energy reference is perfect ' // A 'crystal at this H **")') else write (*, '(/, "** Free energy reference is perfect ' // A 'crystal at zero stress **")') endif endif IF (.NOT.((OPTION.EQ.1).AND.(IW.GE.1))) print *,' ' buf = name(1:index(name,' ')-1) // '.xyz.opt' if (print_atom) call print_atoms(buf) IF (OPTION.EQ.1.AND.IW.NE.10) THEN IW = IW+1 DO I=1,3 DO K=1,3 ETA(I,K) = 0.D0 ENDDO ENDDO GOTO 154 ENDIF IF (OPTION.NE.1) THEN IF (CONST_H) THEN call tersoff_force write (*,'("After relaxation, the stress in cell is (MBar)")') write (*,422) ((STRESS(I,J)*USTRESS_IN_MBAR,J=1,3),I=1,3) ELSE write (*,'(/,"After relaxation, the strain in cell is")') write (*,422) ((ETA(I,J),J=1,3),I=1,3) write (*,'("the H-matrix is")') write (*,422) ((H(I,J),J=1,3),I=1,3) ENDIF ENDIF IF (OPTION.NE.1) THEN C Write on configurational file: OPEN (UNIT=LP_CONFIG, STATUS='UNKNOWN', FORM='UNFORMATTED', a FILE=NAME_WRITE) WRITE (LP_CONFIG) N WRITE (LP_CONFIG) (NTYPE(I),I=1,N) WRITE (LP_CONFIG) ((H(I,J),J=1,3),I=1,3) DO I = 1,N IF (SX(I).LT.0.D0) SX(I)=SX(I)+1.D0 IF (SX(I).GT.1.D0) SX(I)=SX(I)-1.D0 IF (SY(I).LT.0.D0) SY(I)=SY(I)+1.D0 IF (SY(I).GT.1.D0) SY(I)=SY(I)-1.D0 IF (SZ(I).LT.0.D0) SZ(I)=SZ(I)+1.D0 IF (SZ(I).GT.1.D0) SZ(I)=SZ(I)-1.D0 ENDDO WRITE (LP_CONFIG) (SX(I),SY(I),SZ(I),I=1,N) WRITE (LP_CONFIG) (0.D0,0.D0,0.D0,I=1,N) CLOSE (LP_CONFIG) ENDIF STOP END FUNCTION VALUE(XX) INCLUDE 'defect.fh' DOUBLE PRECISION VALUE,XX(MOP),ETA(3,3) EXTERNAL TERSOFF_VALUE IP = 1 DO I=1,N SX(I) = XX(IP) SY(I) = XX(IP+1) SZ(I) = XX(IP+2) IP = IP+3 ENDDO IF (CONST_S) THEN ETA(1,1) = XX(IP) ETA(2,2) = XX(IP+1) ETA(3,3) = XX(IP+2) ETA(1,2) = XX(IP+3) ETA(1,3) = XX(IP+4) ETA(2,3) = XX(IP+5) IP = IP+6 ETA(2,1) = ETA(1,2) ETA(3,1) = ETA(1,3) ETA(3,2) = ETA(2,3) CALL NEW_H (ETA) ENDIF CALL TERSOFF_VALUE VALUE = POTE - WORK_FUNCTION(ETA) IF (OPTION.NE.1) PRINT *, 'value =',VALUE c PRINT *, 'value =',VALUE RETURN END SUBROUTINE FORCE(XX,F) C Should have called VALUE before calling FORCE INCLUDE 'defect.fh' DOUBLE PRECISION XX(MOP),F(MOP),TMP(3,3),TMPP(3,3),MATNORM COMMON /J/ J(3,3),JI(3,3) DOUBLE PRECISION J,JI IF (OPTION.NE.1) PRINT *, '** force evaluation... **' CALL TERSOFF_FORCE IP = 1 DO I=1,N F(IP) = HI(1,1)*FX(I)+HI(2,1)*FY(I)+HI(3,1)*FZ(I) F(IP+1) = HI(1,2)*FX(I)+HI(2,2)*FY(I)+HI(3,2)*FZ(I) F(IP+2) = HI(1,3)*FX(I)+HI(2,3)*FY(I)+HI(3,3)*FZ(I) IP = IP+3 ENDDO IF (CONST_S) THEN CALL MATINV (J,JI,CHANGE) CALL MATADD (SOUT,STRESS,TMP) CALL MATMUL (JI,TMP,TMPP) CALL MATMUL (TMPP,JI,TMP) IF (OPTION.NE.1) A PRINT *,'** stress discrepancy <',MATNORM(TMP),' MBar **' F(IP) = TMP(1,1) * VOLUME F(IP+1) = TMP(2,2) * VOLUME F(IP+2) = TMP(3,3) * VOLUME F(IP+3) = TMP(1,2) * VOLUME * 2.D0 F(IP+4) = TMP(1,3) * VOLUME * 2.D0 F(IP+5) = TMP(2,3) * VOLUME * 2.D0 IP = IP + 6 ENDIF RETURN END SUBROUTINE TERSOFF_VALUE INCLUDE 'defect.fh' 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,N EP(I) = 0. ENDDO C --------------------------------------------------------- C ALWAYS UPDATE THE VERLET LIST DO IPT = 1,N IDX(IPT) = 0 ENDDO DO 170 IPT = 1,N N1 = NTYPE(IPT) SXI = SX(IPT) SYI = SY(IPT) SZI = SZ(IPT) DO 270 JPT = IPT+1,N N2 = NTYPE(JPT) SXJI = SX(JPT) - SXI SYJI = SY(JPT) - SYI SZJI = SZ(JPT) - SZI IF(SXJI.GT.0.5D0) SXJI=SXJI-1.D0 IF(SXJI.LT.-0.5D0) SXJI=SXJI+1.D0 IF(SYJI.GT.0.5D0) SYJI=SYJI-1.D0 IF(SYJI.LT.-0.50) SYJI=SYJI+1.D0 IF(SZJI.GT.0.5D0) SZJI=SZJI-1.D0 IF(SZJI.LT.-0.5D0) SZJI=SZJI+1.D0 RXX = H(1,1)*SXJI + H(2,1)*SYJI + H(3,1)*SZJI RYY = H(1,2)*SXJI + H(2,2)*SYJI + H(3,2)*SZJI RZZ = 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 c IF (IDX(IPT).EQ.NX.OR.IDX(JPT).EQ.NX) print *, c A 'Warning (defect.fh): at H(1,1) =',H(1,1),', NX out of bound.' 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,N N1 = NTYPE(IPT) SXI = SX(IPT) SYI = SY(IPT) SZI = SZ(IPT) NEG = 0 DO 600 J = 1,IDX(IPT) JPT = LIST(IPT,J) N2 = NTYPE(JPT) RSMAX = PS(N1)*PS(N2) SXJI = SX(JPT) - SXI SYJI = SY(JPT) - SYI SZJI = SZ(JPT) - SZI IF(SXJI.GT.0.5D0) SXJI=SXJI-1.D0 IF(SXJI.LT.-0.5D0) SXJI=SXJI+1.D0 IF(SYJI.GT.0.5D0) SYJI=SYJI-1.D0 IF(SYJI.LT.-0.50) SYJI=SYJI+1.D0 IF(SZJI.GT.0.5D0) SZJI=SZJI-1.D0 IF(SZJI.LT.-0.5D0) SZJI=SZJI+1.D0 RXX = H(1,1)*SXJI + H(2,1)*SYJI + H(3,1)*SZJI RYY = H(1,2)*SXJI + H(2,2)*SYJI + H(3,2)*SZJI RZZ = 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,N POTE = POTE+EP(I) ENDDO RETURN END SUBROUTINE TERSOFF_FORCE INCLUDE 'defect.fh' DIMENSION ZET2(NX,NX),ZET3(NX,NX),FC(NM),FCP(NM), A RX(NX),RY(NX),RZ(NX),RXT(NX),RYT(NX),RZT(NX), A R(NX), RJK2(NX,NX),RJK(NX,NX),RJKX(NX,NX),RJKY(NX,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 FRST(NX),FF1(NX),FF2(NX),FF3(NX),R2(NX),STRS21(NX),STRS22(NX), A STRS23(NX),STRS24(NX),STRS25(NX),STRS34(NX),STRS26(NX), A STRS31(NX),STRS32(NX),STRS33(NX),STRS35(NX),STRS36(NX) DO I=1,N FX(I) = 0. FY(I) = 0. FZ(I) = 0. EP(I) = 0. ENDDO STRESS(1,1) = 0.0D0 STRESS(2,2) = 0.0D0 STRESS(3,3) = 0.0D0 STRESS(1,2) = 0.0D0 STRESS(1,3) = 0.0D0 STRESS(2,3) = 0.0D0 DO 500 IPT = 1,N N1 = NTYPE(IPT) SXI = SX(IPT) SYI = SY(IPT) SZI = SZ(IPT) NEG = 0 DO 600 J = 1,IDX(IPT) JPT = LIST(IPT,J) N2 = NTYPE(JPT) RSMAX = PS(N1)*PS(N2) SXJI = SX(JPT) - SXI SYJI = SY(JPT) - SYI SZJI = SZ(JPT) - SZI IF(SXJI.GT.0.5D0) SXJI=SXJI-1.D0 IF(SXJI.LT.-0.5D0) SXJI=SXJI+1.D0 IF(SYJI.GT.0.5D0) SYJI=SYJI-1.D0 IF(SYJI.LT.-0.50) SYJI=SYJI+1.D0 IF(SZJI.GT.0.5D0) SZJI=SZJI-1.D0 IF(SZJI.LT.-0.5D0) SZJI=SZJI+1.D0 RXX = H(1,1)*SXJI + H(1,2)*SYJI + H(1,3)*SZJI RYY = H(2,1)*SXJI + H(2,2)*SYJI + H(2,3)*SZJI RZZ = H(3,1)*SXJI + H(3,2)*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. FCP(NEG) = 0. ELSE FC(NEG)=0.5*(1.+COS(PI*(R1-CR(NEG)) A /(CS(NEG)-CR(NEG)))) FCP(NEG)=-0.5*PI*SIN(PI*(R1-CR(NEG)) A /(CS(NEG)-CR(NEG)))/(CS(NEG)-CR(NEG)) ENDIF 600 CONTINUE DO 700 J = 1,NEG JPT = LIST2(J) N22=NTYPE(JPT) ZET(J) = 0. SZET1(J) = 0. SZET11(J) = 0. SFTM2X(J) = 0. SFTM2Y(J) = 0. SFTM2Z(J) = 0. SFTM3X(J) = 0. SFTM3Y(J) = 0. SFTM3Z(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) RJK2(J,K) = R2(J)+R2(K)-2.*RMU*R(J)*R(K) RJK(J,K) = SQRT(RJK2(J,K)) RJKX(J,K) = RX(J)-RX(K) RJKY(J,K) = RY(J)-RY(K) RJKZ(J,K) = RZ(J)-RZ(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.*RMU1*PC(N1)/RMU3**2 COM1 = FCP(K)*GMU COM2 = FC(K)*GMUP ZET(J) = ZET(J)+FC(K)*GMU TIP1 = RMU/R(J)-1./R(K) TIP2 = RMU/R(K)-1./R(J) TIP3 = RJK(J,K)/R(J)/R(K) ZET1 = COM2*(-TIP1) SZET1(J) = SZET1(J)+ZET1 ZET2(J,K) = COM1+COM2*(-TIP2) ZET3(J,K) = COM2*(-TIP3) SFTM2X(J) = SFTM2X(J) + ZET2(J,K)*RXT(K) SFTM2Y(J) = SFTM2Y(J) + ZET2(J,K)*RYT(K) SFTM2Z(J) = SFTM2Z(J) + ZET2(J,K)*RZT(K) SFTM3X(J) = SFTM3X(J) + ZET3(J,K)*PRJKX(J,K) SFTM3Y(J) = SFTM3Y(J) + ZET3(J,K)*PRJKY(J,K) SFTM3Z(J) = SFTM3Z(J) + ZET3(J,K)*PRJKZ(J,K) 800 CONTINUE 700 CONTINUE DO 704 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) FCD1 = FCP(J)-RLM1(J)*FC(J) FCD2 = FCP(J)-RLM2(J)*FC(J) IF(ZET(J).NE.0.) THEN BEZET = 1.+ PBT(N1)*(ZET(J)**PN(N1)) BIJ = (BEZET**PWR(N1))*CHIJ(J) BIJ1 = 0.5*PBT(N1)*(ZET(J)**PWR2(N1)) BIJP = BIJ1*(BEZET**PWR1(N1))*CHIJ(J) PF(J) = FC(J)*ELM2B*BIJP ELSE BIJ = 1. PF(J) = 0. ENDIF FRST(J) = FCD1*ELM1A - BIJ*FCD2*ELM2B FF1(J) = FRST(J) + PF(J)*SZET1(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. EP(JPT) = EP(JPT) + VIJ/4. STRS21(J) = 0.D0 STRS22(J) = 0.D0 STRS23(J) = 0.D0 STRS24(J) = 0.D0 STRS25(J) = 0.D0 STRS26(J) = 0.D0 STRS31(J) = 0.D0 STRS32(J) = 0.D0 STRS33(J) = 0.D0 STRS34(J) = 0.D0 STRS35(J) = 0.D0 STRS36(J) = 0.D0 C THREE BODY FORCES AND THREE BODY STRESSES DO 840 K=1,NEG KPT = LIST2(K) IF (K.EQ.J) GOTO 840 FF2(J) = PF(J)*ZET2(J,K) FF3(J) = PF(J)*ZET3(J,K) STRS21(J) = STRS21(J)+RXT(K)*RX(K)*FF2(J) STRS22(J) = STRS22(J)+RXT(K)*RY(K)*FF2(J) STRS23(J) = STRS23(J)+RXT(K)*RZ(K)*FF2(J) STRS24(J) = STRS24(J)+RYT(K)*RY(K)*FF2(J) STRS25(J) = STRS25(J)+RYT(K)*RZ(K)*FF2(J) STRS26(J) = STRS26(J)+RZT(K)*RZ(K)*FF2(J) STRS31(J) = STRS31(J)+PRJKX(J,K)*RJKX(J,K)*FF3(J) STRS32(J) = STRS32(J)+PRJKX(J,K)*RJKY(J,K)*FF3(J) STRS33(J) = STRS33(J)+PRJKX(J,K)*RJKZ(J,K)*FF3(J) STRS34(J) = STRS34(J)+PRJKY(J,K)*RJKY(J,K)*FF3(J) STRS35(J) = STRS35(J)+PRJKY(J,K)*RJKZ(J,K)*FF3(J) STRS36(J) = STRS36(J)+PRJKZ(J,K)*RJKZ(J,K)*FF3(J) C CALC FORCE ON PARTICLE K FKX = -(FF2(J)*RX(K)/R(K)+FF3(J)*(-PRJKX(J,K)))*0.5 FKY = -(FF2(J)*RY(K)/R(K)+FF3(J)*(-PRJKY(J,K)))*0.5 FKZ = -(FF2(J)*RZ(K)/R(K)+FF3(J)*(-PRJKZ(J,K)))*0.5 FX(KPT) = FX(KPT) + FKX FY(KPT) = FY(KPT) + FKY FZ(KPT) = FZ(KPT) + FKZ 840 CONTINUE C -------------------------------------------------------------- C CALC FORCE ON PARTICLE I FIX = ((FRST(J)+PF(J)*SZET1(J))*RXT(J)+PF(J)*SFTM2X(J))*0.5 FIY = ((FRST(J)+PF(J)*SZET1(J))*RYT(J)+PF(J)*SFTM2Y(J))*0.5 FIZ = ((FRST(J)+PF(J)*SZET1(J))*RZT(J)+PF(J)*SFTM2Z(J))*0.5 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.5 FJY=-((FRST(J)+PF(J)*SZET1(J))*RYT(J)+SFTM3Y(J)*PF(J))*0.5 FJZ=-((FRST(J)+PF(J)*SZET1(J))*RZT(J)+SFTM3Z(J)*PF(J))*0.5 FX(JPT) = FX(JPT)+FJX FY(JPT) = FY(JPT)+FJY FZ(JPT) = FZ(JPT)+FJZ C CALC STRESS TENSOR OF THE SYSTEM. STRESS(1,1) = STRESS(1,1) - A (FF1(J)*RXT(J)*RX(J)+STRS21(J)+STRS31(J))/2.D0 STRESS(2,2) = STRESS(2,2) - A (FF1(J)*RYT(J)*RY(J)+STRS24(J)+STRS34(J))/2.D0 STRESS(3,3) = STRESS(3,3) - A (FF1(J)*RZT(J)*RZ(J)+STRS26(J)+STRS36(J))/2.D0 STRESS(1,2) = STRESS(1,2) - A (FF1(J)*RXT(J)*RY(J)+STRS22(J)+STRS32(J))/2.D0 STRESS(1,3) = STRESS(1,3) - A (FF1(J)*RXT(J)*RZ(J)+STRS23(J)+STRS33(J))/2.D0 STRESS(2,3) = STRESS(2,3) - A (FF1(J)*RYT(J)*RZ(J)+STRS25(J)+STRS35(J))/2.D0 704 CONTINUE 500 CONTINUE POTE = 0. DO I = 1,N POTE = POTE + EP(I) ENDDO STRESS(1,1) = STRESS(1,1) / VOLUME STRESS(2,2) = STRESS(2,2) / VOLUME STRESS(3,3) = STRESS(3,3) / VOLUME STRESS(1,2) = STRESS(1,2) / VOLUME STRESS(1,3) = STRESS(1,3) / VOLUME STRESS(2,3) = STRESS(2,3) / VOLUME STRESS(2,1) = STRESS(1,2) STRESS(3,1) = STRESS(1,3) STRESS(3,2) = STRESS(2,3) RETURN END SUBROUTINE NEW_H (ETA) INCLUDE 'defect.fh' DIMENSION ETA(3,3) COMMON /J/ J(3,3), JI(3,3) DOUBLE PRECISION J,JI CALL ETA_TO_J (ETA,J) CALL MATMUL (H0, J, H) CALL MATINV (H, HI, VOLUME) RETURN END SUBROUTINE ETA_TO_J (ETA,J) DOUBLE PRECISION ETA(3,3),J(3,3),TMP(3,3),TMPP(3,3),MATNORM DO I=1,3 DO K=1,3 J(I,K) = ETA(I,K) IF (I.EQ.K) J(I,K) = J(I,K) + 1.D0 ENDDO ENDDO c 2nd order CALL MATMUL (ETA,ETA,TMP) CALL MATEQV (TMP,TMPP) CALL MATADD_MUL (J,-0.5D0,TMP,J) c 3rd order CALL MATMUL (TMPP,ETA,TMP) CALL MATEQV (TMP,TMPP) CALL MATADD_MUL (J,+0.5D0,TMP,J) c 4th order CALL MATMUL (TMPP,ETA,TMP) CALL MATEQV (TMP,TMPP) CALL MATADD_MUL (J,-5.D0/8.D0,TMP,J) c 5th order CALL MATMUL (TMPP,ETA,TMP) CALL MATEQV (TMP,TMPP) CALL MATADD_MUL (J,+7.D0/8.D0,TMP,J) c 6th order c CALL MATMUL (TMPP,ETA,TMP) c CALL MATEQV (TMP,TMPP) c CALL MATADD_MUL (J,-21.D0/16.D0,TMP,J) c 7th order c CALL MATMUL (TMPP,ETA,TMP) c CALL MATEQV (TMP,TMPP) c CALL MATADD_MUL (J,33.D0/16.D0,TMP,J) c 8th order c CALL MATMUL (TMPP,ETA,TMP) c CALL MATEQV (TMP,TMPP) c CALL MATADD_MUL (J,-429.D0/128.D0,TMP,J) c 9th order c CALL MATMUL (TMPP,ETA,TMP) c CALL MATEQV (TMP,TMPP) c CALL MATADD_MUL (J,715.D0/128.D0,TMP,J) RETURN END DOUBLE PRECISION FUNCTION DRIVING_FORCE (ETA,T) INCLUDE 'defect.fh' DOUBLE PRECISION ETA(3,3),ETAT(3,3),J(3,3), A JI(3,3),TMP(3,3),TMPP(3,3),MATNORM DO I=1,3 DO K=1,3 ETAT(I,K) = ETA(I,K)*T ENDDO ENDDO CALL ETA_TO_J (ETAT,J) CALL MATMUL (H0,J,TMP) CALL MATINV (TMP,TMPP,VOLUMET) CALL MATINV (J,JI,CHANGE) CALL MATMUL (JI,SOUT,TMPP) CALL MATMUL (TMPP,JI,TMP) DRIVING_FORCE = 0.D0 DO I=1,3 DO K=1,3 DRIVING_FORCE = DRIVING_FORCE+TMP(I,K)*ETA(I,K)*VOLUMET ENDDO ENDDO RETURN END DOUBLE PRECISION FUNCTION WORK_FUNCTION (ETA) C ************************************************************** C Calculate the work done by constant stress in symmetric C deformation space by straight path Romberg integration: C EPS is the desired accuracy. C ************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (EPS=1.d-10,JMAX=25,JMAXP=JMAX+1,K=13,KM=K-1) DOUBLE PRECISION ETA(3,3),MATNORM DOUBLE PRECISION h(JMAXP),s(JMAXP) if (matnorm(eta).gt.0.5d0) then work_function = -matnorm(eta)**2 return endif h(1) = 1.d0 s(1) = 0.5d0*(driving_force(eta,0.d0)+driving_force(eta,1.d0)) do 12 j = 2,JMAX C Integration by trapzoidal rule with progressively decreasing h h(j)=0.25d0*h(j-1) it=2**(j-2) tnm=dble(it) del=1.d0/tnm t=0.5d0*del sum=0.d0 do jj=1,it sum=sum+driving_force(eta,t) t=t+del enddo s(j)=0.5d0*(s(j-1)+sum/tnm) if (j.ge.K) then call polint(h(j-KM),s(j-KM),K,0.d0,WORK_FUNCTION,error) c print *, abs(error) / abs(WORK_FUNCTION) if (abs(error).le.EPS*abs(WORK_FUNCTION)) then return endif endif 12 continue STOP A 'Error(work_function): can not converge in Romberg integration.' END SUBROUTINE INITIALIZE_CONSTANTS INCLUDE 'defect.fh' UENERGY_IN_JOULE = 1.602192D-19 ULENGTH_IN_METER = 1.D-10 UMASS_IN_KG = 1.66053D-27 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 AVO = 6.0221367D+23 BOLZ = 1.380658D-23 HBAR = 1.054589D-34 PM(1) = 28.086 PM(2) = 12.011 NAME_ATOM(1) = 'Si' PA(1) = 1.8308E3 PB(1) = 4.7118E2 PL1(1) = 2.4799 PL2(1) = 1.7322 PN(1) = 7.8734D-1 PBT(1) = (1.1D-6)**PN(1) PC(1) = (1.0039E5)**2. PD(1) = (1.6217E1)**2. PH(1) = -5.9825D-1 PR(1) = 2.36 PS(1) = 2.56 C2D2(1) = 1.+PC(1)/PD(1) PWR(1) = -1./(2.*PN(1)) PWR1(1) = PWR(1)-1. PWR2(1) = PN(1)-1. NAME_ATOM(2) = 'C ' PA(2) = 1.5448E3 PB(2) = 3.8963E2 PL1(2) = 3.4653 PL2(2) = 2.3064 PN(2) = 9.9054D-1 PBT(2) = (4.1612D-6)**PN(2) PC(2) = (1.9981E4)**2. PD(2) = (7.0340)**2. PH(2) = -3.9953D-1 PR(2) = 2.36 PS(2) = 2.56 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 = 1.0086 C Meijie's parameter: PRSIC = 2.36D0 PSSIC = 2.56D0 RCUT = MAX(PS(1),PS(2)) RETURN END SUBROUTINE PRINT_ATOMS (NAME_XYZ) INCLUDE 'defect.fh' CHARACTER *70 NAME_XYZ DIMENSION RXX(NM),RYY(NM),RZZ(NM),IX(NM) CHARACTER*2 NAME_OF_THEM(NM) IF (OPTION.NE.1) THEN PRINT *,'-------------- atomic positions print-out --------------' PRINT*,'The atom',iatom,' at origin is a ',NAME_ATOM(ntype(iatom)) WRITE (*,'(" of (X,Y,Z) =",3F11.5)') X(IATOM), Y(IATOM), Z(IATOM) PRINT *,'--------------------------------------------------------' PRINT *,' Type Seq. DX DY DZ R ' ENDIF IP = 0 DO J = 1,N SXJP = SX(J) - SX(IATOM) SYJP = SY(J) - SY(IATOM) SZJP = SZ(J) - SZ(IATOM) IF(SXJP.GT.0.5D0) SXJP=SXJP-1.D0 IF(SXJP.LT.-0.5D0) SXJP=SXJP+1.D0 IF(SYJP.GT.0.5D0) SYJP=SYJP-1.D0 IF(SYJP.LT.-0.5D0) SYJP=SYJP+1.D0 IF(SZJP.GT.0.5D0) SZJP=SZJP-1.D0 IF(SZJP.LT.-0.5D0) SZJP=SZJP+1.D0 RX = H(1,1)*SXJP + H(2,1)*SYJP + H(3,1)*SZJP RY = H(1,2)*SXJP + H(2,2)*SYJP + H(3,2)*SZJP RZ = H(1,3)*SXJP + H(2,3)*SYJP + H(3,3)*SZJP R2 = RX*RX+RY*RY+RZ*RZ IF (R2.LT.RPRINT**2) THEN R1 = SQRT(R2) IF (OPTION.NE.1) WRITE (*,'(3x,A2,I7,4F11.6)') A NAME_ATOM(NTYPE(J)),J,RX,RY,RZ,R1 IP = IP+1 IX(IP) = J RXX(IP) = RX RYY(IP) = RY RZZ(IP) = RZ ENDIF ENDDO OPEN (UNIT=LP_XYZ, STATUS='UNKNOWN', FORM='FORMATTED', A FILE=NAME_XYZ) WRITE (LP_XYZ,*) IP WRITE (LP_XYZ,'(" Defect ",A6)') NAME(1:6) DO 321 I=1,IP 321 WRITE (LP_XYZ,'(3X,A2,3(F12.6))') A NAME_ATOM(NTYPE(IX(I))),RXX(I),RYY(I),RZZ(I) CLOSE (LP_XYZ) IF (OPTION.NE.1) THEN PRINT *,'--------------------------------------------------------' PRINT *, ' ' ENDIF RETURN END SUBROUTINE MATADD (A,B,C) DOUBLE PRECISION A(3,3), B(3,3), C(3,3) C(1,1) = A(1,1) + B(1,1) C(1,2) = A(1,2) + B(1,2) C(1,3) = A(1,3) + B(1,3) C(2,1) = A(2,1) + B(2,1) C(2,2) = A(2,2) + B(2,2) C(2,3) = A(2,3) + B(2,3) C(3,1) = A(3,1) + B(3,1) C(3,2) = A(3,2) + B(3,2) C(3,3) = A(3,3) + B(3,3) RETURN END SUBROUTINE MATADD_MUL (A,CC,B,C) DOUBLE PRECISION A(3,3), B(3,3), C(3,3), CC C(1,1) = A(1,1) + CC * B(1,1) C(1,2) = A(1,2) + CC * B(1,2) C(1,3) = A(1,3) + CC * B(1,3) C(2,1) = A(2,1) + CC * B(2,1) C(2,2) = A(2,2) + CC * B(2,2) C(2,3) = A(2,3) + CC * B(2,3) C(3,1) = A(3,1) + CC * B(3,1) C(3,2) = A(3,2) + CC * B(3,2) C(3,3) = A(3,3) + CC * B(3,3) RETURN END SUBROUTINE MATINV (A,B,C) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(3,3),B(3,3),C D11=A(2,2)*A(3,3)-A(2,3)*A(3,2) D22=A(3,3)*A(1,1)-A(3,1)*A(1,3) D33=A(1,1)*A(2,2)-A(1,2)*A(2,1) D12=A(2,3)*A(3,1)-A(2,1)*A(3,3) D23=A(3,1)*A(1,2)-A(3,2)*A(1,1) D31=A(1,2)*A(2,3)-A(1,3)*A(2,2) D13=A(2,1)*A(3,2)-A(3,1)*A(2,2) D21=A(3,2)*A(1,3)-A(1,2)*A(3,3) D32=A(1,3)*A(2,1)-A(2,3)*A(1,1) C=A(1,1)*D11+A(1,2)*D12+A(1,3)*D13 B(1,1)=D11/C B(2,2)=D22/C B(3,3)=D33/C B(1,2)=D21/C B(2,3)=D32/C B(3,1)=D13/C B(2,1)=D12/C B(3,2)=D23/C B(1,3)=D31/C RETURN END SUBROUTINE MATMUL (A,B,C) DOUBLE PRECISION A(3,3), B(3,3), C(3,3) C(1,1) = A(1,1)*B(1,1) + A(1,2)*B(2,1) + A(1,3)*B(3,1) C(1,2) = A(1,1)*B(1,2) + A(1,2)*B(2,2) + A(1,3)*B(3,2) C(1,3) = A(1,1)*B(1,3) + A(1,2)*B(2,3) + A(1,3)*B(3,3) C(2,1) = A(2,1)*B(1,1) + A(2,2)*B(2,1) + A(2,3)*B(3,1) C(2,2) = A(2,1)*B(1,2) + A(2,2)*B(2,2) + A(2,3)*B(3,2) C(2,3) = A(2,1)*B(1,3) + A(2,2)*B(2,3) + A(2,3)*B(3,3) C(3,1) = A(3,1)*B(1,1) + A(3,2)*B(2,1) + A(3,3)*B(3,1) C(3,2) = A(3,1)*B(1,2) + A(3,2)*B(2,2) + A(3,3)*B(3,2) C(3,3) = A(3,1)*B(1,3) + A(3,2)*B(2,3) + A(3,3)*B(3,3) RETURN END SUBROUTINE MATEQV (A,B) DOUBLE PRECISION A(3,3),B(3,3) B(1,1) = A(1,1) B(1,2) = A(1,2) B(1,3) = A(1,3) B(2,1) = A(2,1) B(2,2) = A(2,2) B(2,3) = A(2,3) B(3,1) = A(3,1) B(3,2) = A(3,2) B(3,3) = A(3,3) RETURN END SUBROUTINE MATSUB (A,B,C) DOUBLE PRECISION A(3,3), B(3,3), C(3,3) C(1,1) = A(1,1) - B(1,1) C(1,2) = A(1,2) - B(1,2) C(1,3) = A(1,3) - B(1,3) C(2,1) = A(2,1) - B(2,1) C(2,2) = A(2,2) - B(2,2) C(2,3) = A(2,3) - B(2,3) C(3,1) = A(3,1) - B(3,1) C(3,2) = A(3,2) - B(3,2) C(3,3) = A(3,3) - B(3,3) RETURN END DOUBLE PRECISION FUNCTION MATNORM(A) DOUBLE PRECISION A(3,3) INTEGER I,K MATNORM = 0.D0 DO I=1,3 DO K=1,3 MATNORM = MATNORM + A(I,K)**2 ENDDO ENDDO MATNORM = SQRT(MATNORM) RETURN END SUBROUTINE RANDOM INCLUDE 'defect.fh' PARAMETER (K=5701,J=3612,M=566927,RM=566927.) IX = INT(SEED*RM) KR = J*IX + K IRAND = MOD(KR,M) SEED = (FLOAT(IRAND)+0.5)/RM RETURN END SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) EXTERNAL func DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY PARAMETER (GOLD=1.618034, GLIMIT=1000.d0, TINY=1.d-20) DOUBLE PRECISION dum,fu,q,r,u,ulim fa=func(ax) fb=func(bx) if(fb.gt.fa)then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif cx=bx+GOLD*(bx-ax) fc=func(cx) 1 if(fb.ge.fc)then r=(bx-ax)*(fb-fc) q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) if((bx-u)*(u-cx).gt.0.)then fu=func(u) if(fu.lt.fc)then ax=bx fa=fb bx=u fb=fu return else if(fu.gt.fb)then cx=u fc=fu return endif u=cx+GOLD*(cx-bx) fu=func(u) else if((cx-u)*(u-ulim).gt.0.)then fu=func(u) if(fu.lt.fc)then bx=cx cx=u u=cx+GOLD*(cx-bx) fb=fc fc=fu fu=func(u) endif else if((u-ulim)*(ulim-cx).ge.0.)then u=ulim fu=func(u) else u=cx+GOLD*(cx-bx) fu=func(u) endif ax=bx bx=cx cx=u fa=fb fb=fc fc=fu goto 1 endif return END FUNCTION brent(ax,bx,cx,f,tol,xmin) INTEGER ITMAX EXTERNAL f DOUBLE PRECISION brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS PARAMETER (ITMAX=1000,CGOLD=.3819660,ZEPS=1.0d-20) INTEGER iter DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1, A tol2,u,v,w,x,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0. fx=f(x) fv=fx fw=fx do 11 iter=1,ITMAX xm=0.5*(a+b) tol1=tol*abs(x)+ZEPS tol2=2.*tol1 if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3 if(abs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.*(q-r) if(q.gt.0.) p=-p q=abs(q) etemp=e e=d if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) *goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(abs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=f(u) if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif 11 continue pause 'brent exceed maximum iterations' 3 xmin=x brent=fx return END FUNCTION f1dim(x) INTEGER NM,NMAX DOUBLE PRECISION f1dim,VALUE,x PARAMETER (NMAX=100000) CU USES VALUE INTEGER j,ncom DOUBLE PRECISION pcom(NMAX),xicom(NMAX),xt(NMAX) COMMON /f1com/ pcom,xicom,ncom do 11 j=1,ncom xt(j)=pcom(j)+x*xicom(j) 11 continue f1dim=VALUE(xt) return END SUBROUTINE linmin(p,xi,n,fret) INTEGER n,NM,NMAX DOUBLE PRECISION fret,p(n),xi(n),TOL PARAMETER (NMAX=100000,TOL=2.d-4) CU USES brent,f1dim,mnbrak INTEGER j,ncom DOUBLE PRECISION ax,bx,fa,fb,fx,xmin,xx,pcom(NMAX), A xicom(NMAX),brent COMMON /f1com/ pcom,xicom,ncom EXTERNAL f1dim DOUBLE PRECISION f1dim ncom=n do 11 j=1,n pcom(j)=p(j) xicom(j)=xi(j) 11 continue ax=0. xx=1. call mnbrak(ax,xx,bx,fa,fx,fb,f1dim) fret=brent(ax,xx,bx,f1dim,TOL,xmin) do 12 j=1,n xi(j)=xmin*xi(j) p(j)=p(j)+xi(j) 12 continue return END SUBROUTINE frprmn (p,n,ftol,iter,fret) INTEGER iter,n,NM,NMAX,ITMAX EXTERNAL VALUE DOUBLE PRECISION fret,ftol,p(n),EPS,VALUE PARAMETER (NMAX=100000,ITMAX=10000,EPS=1.d-20) CU USES FORCE,VALUE,linmin INTEGER its,j DOUBLE PRECISION dgg,fp,gam,gg,g(NMAX),h(NMAX),xi(NMAX) if (n.gt.NMAX) then print *, 'error (frprmn): too many variables, revise NMAX >',n stop endif fp=VALUE(p) call FORCE(p,xi) do 11 j=1,n g(j)=-xi(j) h(j)=g(j) xi(j)=h(j) 11 continue do 14 its=1,ITMAX iter=its call linmin(p,xi,n,fret) if(2.*abs(fret-fp).le.ftol*(abs(fret)+abs(fp)+EPS)) return fp=VALUE(p) call FORCE(p,xi) gg=0. dgg=0. do 12 j=1,n gg=gg+g(j)**2 C dgg=dgg+xi(j)**2 dgg=dgg+(xi(j)+g(j))*xi(j) 12 continue if(gg.eq.0.)return gam=dgg/gg do 13 j=1,n g(j)=-xi(j) h(j)=g(j)+gam*h(j) xi(j)=h(j) 13 continue 14 continue pause 'frprmn maximum iterations exceeded' return END SUBROUTINE polint(xa,ya,n,x,y,dy) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER n,NMAX DOUBLE PRECISION dy,x,y,xa(n),ya(n) PARAMETER (NMAX=20) INTEGER i,m,ns DOUBLE PRECISION den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) ns=1 dif=dabs(x-xa(1)) do 11 i=1,n dift=dabs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.) pause 'failure in polint' den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END