C ------------------------------------------------------------ C PROGRAM FitAll V1.1 C AN OPTIMIZATION PROGRAM TO FIT THE ELECTRONIC BAND STRUCTURE C AND COHESIVE ENERGIES OF BINARY SI-C SYSTEMS USING EMPIRICAL C TIGHT-BINDING MODEL C 1. ABLE TO HANDLE A DATABASE OF ARBITRARY STRUCTURES C 2. ENVIRONMENT-DEPENDENT TB MODEL (PRB 53, 979, 1996) C 3. GENERIC PARAMETER TABLE FOR THE BINARY SYSTEM C C AUG.7 1998 C JU LI (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 ------------------------------------------------------------ PROGRAM FitAll # include "fitall.fh" C MINA and ASA module DOUBLE PRECISION LOWER_BOUND(NOP),UPPER_BOUND(NOP),GUESS(NOP) EXTERNAL ASA_MAIN,FE ITOT = 1 WRAPPING_UP = .FALSE. DO_EBS = .TRUE. C ---------------------------------------------------------------- C These are going to be optimized: READ (*,'("Si-C Parameters (Upper_bound Guess Lower_bound):")') IP = 1 DO I = 1,13 READ (*,'(A70)') NAME_PARA(I) READ *,BOUND(IP,2),BOUND(IP+1,2),BOUND(IP+2,2),BOUND(IP+3,2) READ *,GUESS(IP), GUESS(IP+1), GUESS(IP+2), GUESS(IP+3) READ *,BOUND(IP,1),BOUND(IP+1,1),BOUND(IP+2,1),BOUND(IP+3,1) IP = IP+4 ENDDO C ---------------------------------------------------------------- C ---------------------------------------------------------------- C These are not going to be optimized: READ (*,'(/,"These parameters are not going to be optimized:")') READ (*,'("For Si:")') READ (*,'("a1 a2 a3 a4")') READ *, (A1(IOP),A2(IOP),A3(IOP),A4(IOP),DELTA(IOP,1,1), A B1(IOP,1),B2(IOP,1),B3(IOP,1), IOP = 2,7) READ (*,'(/,"For C:")') READ (*,'("a1 a2 a3 a4")') READ *, (A1(IOP),A2(IOP),A3(IOP),A4(IOP),DELTA(IOP,2,2), A B1(IOP,2),B2(IOP,2),B3(IOP,2), IOP = 9,14) READ (*,'(/,"Polynomial coefficients for Si:")') READ *, C0, C1, C2, C3, C4 READ (*,'(/,"Polynomial coefficients for C:")') READ *, D0, D1, D2, D3, D4 READ (*,'(/,"The Interaction Cutoffs are:")') READ *, RCUTSI_L,RCUTSI_H,RCUTC_L,RCUTC_H,RCUTSIC_L,RCUTSIC_H RCUT(1,1,1) = RCUTSI_L RCUT(1,2,1) = RCUTSIC_L RCUT(2,1,1) = RCUTSIC_L RCUT(2,2,1) = RCUTC_L RCUT(1,1,2) = RCUTSI_H RCUT(1,2,2) = RCUTSIC_H RCUT(2,1,2) = RCUTSIC_H RCUT(2,2,2) = RCUTC_H C ---------------------------------------------------------------- READ (*,'(/,"Number of Object Groups: ",I)') NOBJ NPHASE = 0 TOT_COH_WEIGHT = 0. TOT_BAND_WEIGHT = 0. DO 10 IOBJ = 1,NOBJ READ (*,'(/,"Object Group Name: ",A70)') NAME(IOBJ) 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,*) ((A(IOBJ,I,J),J=1,3),I=1,3) C Invert the matrix: D11 = A(IOBJ,2,2)*A(IOBJ,3,3)-A(IOBJ,2,3)*A(IOBJ,3,2) D22 = A(IOBJ,3,3)*A(IOBJ,1,1)-A(IOBJ,3,1)*A(IOBJ,1,3) D33 = A(IOBJ,1,1)*A(IOBJ,2,2)-A(IOBJ,1,2)*A(IOBJ,2,1) D12 = A(IOBJ,2,3)*A(IOBJ,3,1)-A(IOBJ,2,1)*A(IOBJ,3,3) D23 = A(IOBJ,3,1)*A(IOBJ,1,2)-A(IOBJ,3,2)*A(IOBJ,1,1) D31 = A(IOBJ,1,2)*A(IOBJ,2,3)-A(IOBJ,1,3)*A(IOBJ,2,2) D13 = A(IOBJ,2,1)*A(IOBJ,3,2)-A(IOBJ,3,1)*A(IOBJ,2,2) D21 = A(IOBJ,3,2)*A(IOBJ,1,3)-A(IOBJ,1,2)*A(IOBJ,3,3) D32 = A(IOBJ,1,3)*A(IOBJ,2,1)-A(IOBJ,2,3)*A(IOBJ,1,1) DET(IOBJ) = A(IOBJ,1,1)*D11+A(IOBJ,1,2)*D12+A(IOBJ,1,3)*D13 P(IOBJ,1,1) = D11/DET(IOBJ) P(IOBJ,2,2) = D22/DET(IOBJ) P(IOBJ,3,3) = D33/DET(IOBJ) P(IOBJ,1,2) = D21/DET(IOBJ) P(IOBJ,2,3) = D32/DET(IOBJ) P(IOBJ,3,1) = D13/DET(IOBJ) P(IOBJ,2,1) = D12/DET(IOBJ) P(IOBJ,3,2) = D23/DET(IOBJ) P(IOBJ,1,3) = D31/DET(IOBJ) C "Number of atoms in the unit cell:" READ (LP,'(A70,/,I,/,A70)') BUF,NPA(IOBJ),BUF C "Type X Y Z" READ (LP,*) (NTYPE(IOBJ,N),(TAU(IOBJ,N,I),I=1,3),N=1,NPA(IOBJ)) CLOSE (LP) READ (*,'("IBZ k-point file: ",A70)') BUF PRINT *,'Loading IBZ k-points from ', BUF(1:INDEX(BUF,' ')-1) OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED',FILE=BUF) READ (LP,*) NZKP(IOBJ) READ (LP,*) ((ZKP(IOBJ,NK,I),I=1,4),NK=1,NZKP(IOBJ)) CLOSE (LP) READ (*,'("Ground weight: ",F)') GROUND_WEIGHT(IOBJ) READ (*,'("Curvature weight: ",F)') CURV_WEIGHT(IOBJ) TOT_COH_WEIGHT = TOT_COH_WEIGHT + A GROUND_WEIGHT(IOBJ) + CURV_WEIGHT(IOBJ) READ (*,'("Number of Characteristic Lengths: ",I)') N NPHASE_OBJ(IOBJ,1) = NPHASE+1 NPHASE_OBJ(IOBJ,2) = NPHASE+N NPHASE = NPHASE_OBJ(IOBJ,2) READ (*,'("Length Band Band wgt.")') DO 20 IPHASE = NPHASE_OBJ(IOBJ,1),NPHASE_OBJ(IOBJ,2) READ (*,'(F16,A16,F)') A AA(IPHASE), NAME_BAND(IPHASE), WB(IPHASE) MYOBJ(IPHASE) = IOBJ C Initialize the CPS cohesive energy target (see note.1): BUF = 'Coh_Database/' // A NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1) // '.cps' OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED',FILE=BUF) READ (LP, *) NDATA, NDATA, NDATA DO N = 1,NDATA READ (LP, *) AA_LDA, VOLUME_LDA, CPS(IPHASE) IF (ABS(AA_LDA-AA(IPHASE)).LT.EPS) GOTO 30 ENDDO PRINT *,'CPS energy at AA =',AA(IPHASE),'not found in ', A BUF(1:INDEX(BUF,' ')-1) STOP 30 PRINT *,'CPS energy found in ', BUF(1:INDEX(BUF,' ')-1) WRITE (*, '(A8," No.", I2, ": Length = ",F9.6, A " A, CPS = ",F9.6," eV.")') A NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1), A IPHASE-NPHASE_OBJ(IOBJ,1)+1, AA(IPHASE), CPS(IPHASE) CLOSE (LP) C Initialize the band structure target: IF (NAME_BAND(IPHASE).EQ.'NULL') WB(IPHASE)=0. TOT_BAND_WEIGHT = TOT_BAND_WEIGHT+WB(IPHASE) IF (WB(IPHASE).GT.0) THEN PRINT *,'Loading band structure targets from Band_Database/' A //NAME_BAND(IPHASE)(1:index(NAME_BAND(IPHASE),' ')-1) OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED', A FILE='Band_Database/'//NAME_BAND(IPHASE)) READ (LP,'("Constant of length: ",F)') CC IF (ABS(CC-AA(IPHASE)).GT.EPS) THEN PRINT *,'The description of band ', NAME_BAND(IPHASE) PRINT *,'is different from Band_Database/'//NAME_BAND(IPHASE) STOP ENDIF READ (LP, '("Number of target k-points: ",I)') NBKP(IPHASE) READ (LP, '("Number of links in band structure plot: ",I)') A LINK(IPHASE) READ (LP, '("From To Sections")') READ (LP, *) A (IA(IPHASE,I),IB(IPHASE,I),NS(IPHASE,I),I=1,LINK(IPHASE)) READ (LP, '("LDA gamma point energy: ",F)') GAMMA_LDA(IPHASE) READ (LP, '("LDA Fermi surface: ",F)') FERMI_LDA(IPHASE) DO NK = 1,NBKP(IPHASE) READ (LP,'(A70,/,A70)') BUF,BUF READ (LP, *) NBA(IPHASE,NK),(BKP(IPHASE,NK,I),I=1,3) BKP(IPHASE,NK,1) = BKP(IPHASE,NK,1)*PI/AA(IPHASE) BKP(IPHASE,NK,2) = BKP(IPHASE,NK,2)*PI/AA(IPHASE) BKP(IPHASE,NK,3) = BKP(IPHASE,NK,3)*PI/AA(IPHASE) READ (LP, '(A70)') BUF READ (LP, *) (NTH(IPHASE,NK,I),BTARGET(IPHASE,NK,I), A BWEIGHT(IPHASE,NK,I),I=1,NBA(IPHASE,NK)) ENDDO CLOSE(LP) ENDIF C Initialize the structure: CALL INITIALIZE_STRUCTURE(IPHASE) 20 CONTINUE PRINT *,' ' 10 CONTINUE READ (*,'(/,"Number of Deformation Modes: ",I)') NDEF DO IDEF = 1,NDEF READ (*,'(/,"Mode name: ",A70)') NAME_DEF(IDEF) READ (*,'("Mode weight: ",F)') DEF_WEIGHT(IDEF) READ (*,'("Deformation elements: ",I)') NELE(IDEF) TOT_COH_WEIGHT = TOT_COH_WEIGHT + DEF_WEIGHT(IDEF) DO 40 I = 1,NELE(IDEF) READ (*,'(A70)') BUF C Look for this element in all object groups: DO IOBJ = 1,NOBJ IF (NAME(IOBJ).EQ.BUF) THEN C Only the first structure in each object group IELE(IDEF,I) = NPHASE_OBJ(IOBJ,1) C Make sure that the program calculates all ground C states if DEF_WEIGHT(IDEF) != 0 IF (DEF_WEIGHT(IDEF).GT.0.AND. A GROUND_WEIGHT(IOBJ).LE.0.AND. A CURV_WEIGHT(IOBJ).LE.0) GROUND_WEIGHT(IOBJ) = EPS GOTO 40 ENDIF ENDDO PRINT *, 'Cannot find ', BUF(1:INDEX(BUF,' ')-1), A ' in all above defined object groups.' STOP 40 CONTINUE ENDDO C Renormalize the IBZ k-point sampling weights: DO IOBJ = 1,NOBJ SUMWEIGHT = 0.D0 DO I = 1,NZKP(IOBJ) SUMWEIGHT = SUMWEIGHT + ZKP(IOBJ,I,4) ENDDO DO I = 1,NZKP(IOBJ) ZKP(IOBJ,I,4) = ZKP(IOBJ,I,4) / SUMWEIGHT / 4 / NPA(IOBJ) ENDDO ENDDO C Renormalize the cohesive energy weights: DO IOBJ = 1,NOBJ GROUND_WEIGHT(IOBJ) = GROUND_WEIGHT(IOBJ)/TOT_COH_WEIGHT CURV_WEIGHT(IOBJ) = CURV_WEIGHT(IOBJ)/TOT_COH_WEIGHT ENDDO C Renormalize the deformation mode weights: DO IDEF = 1,NDEF DEF_WEIGHT(IDEF) = DEF_WEIGHT(IDEF)/TOT_COH_WEIGHT ENDDO C Renormalize the band structure weights: DO IPHASE = 1,NPHASE IF (WB(IPHASE).GT.0) THEN THIS_BAND_TOTAL = 0. DO NK = 1,NBKP(IPHASE) DO I = 1,NBA(IPHASE,NK) THIS_BAND_TOTAL = THIS_BAND_TOTAL + BWEIGHT(IPHASE,NK,I) ENDDO ENDDO DO NK = 1,NBKP(IPHASE) DO I = 1,NBA(IPHASE,NK) BWEIGHT(IPHASE,NK,I) = BWEIGHT(IPHASE,NK,I) / THIS_BAND_TOTAL A * WB(IPHASE) / TOT_BAND_WEIGHT ENDDO ENDDO ENDIF ENDDO READ (*,'(/,"Calculate EBS only initially (y/n): ",A70)') BUF FIX_EBS = (BUF.EQ.'y'.OR.BUF.EQ.'Y') READ (*,'("Band error ratio after divided by eV: ",F)') A BAND_ERROR_RATIO READ (*,'("Minimizer (MINA=1; ASA=2): ",I)') MINIMIZER READ (*,'("DEL: ",F)') DEL READ (*,'("NDIV: ",I)') NDIV READ (*,'("Annealing rate: ",F)') ANNEAL_RATE READ (*,'("Maximum number of iterations:",I)') MAXITER READ (*,'("Frequency of saving para: ",I)') IOFTEN OPEN (UNIT=LP_PARA, STATUS='UNKNOWN', FORM='FORMATTED', A FILE="para") IF (MAXITER.NE.0) THEN IF (MINIMIZER.EQ.1) THEN C The MINA optimizer CALL MINA (fe, nop, ndiv, del, bound, guess, guess, emin, ierr) ELSE C Try the current parameter set: CC = FE(GUESS) C The ASA optimizer DO I = 1,NOP LOWER_BOUND(I) = BOUND(I,1) UPPER_BOUND(I) = BOUND(I,2) ENDDO CALL ASA_MAIN (emin, nop, lower_bound, upper_bound, guess, A anneal_rate, main_exit_code) ENDIF ENDIF WRAPPING_UP = .TRUE. CC = FE(GUESS) PRINT *,'The average error for this parameter set is', CC PRINT *,' ' CALL WRAPUP STOP END SUBROUTINE WRAPUP C ------------------------------------------------------------- C Calculate the properties using the current parameter set. C ------------------------------------------------------------- # include "fitall.fh" PRINT *,'Start property calculations... ' PRINT *,' ' DO 10 IPHASE = 1,NPHASE C Plot out the band structure using optimized parameters IF (NAME_BAND(IPHASE).NE.'NULL') THEN BUF = NAME_BAND(IPHASE)(1:INDEX(NAME_BAND(IPHASE),' ')-1) A // '.band' PRINT *,'Plotting out the band structure for ', A NAME_BAND(IPHASE)(1:INDEX(NAME_BAND(IPHASE),' ')-1), A ' on ', BUF(1:INDEX(BUF,' ')-1), ' ..' CALL OVERLAP (IPHASE) OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED',FILE=BUF) RKTOT = 0.D0 DO I=1,LINK(IPHASE) IKA = IA(IPHASE,I) IKB = IB(IPHASE,I) DKX = (BKP(IPHASE,IKB,1)-BKP(IPHASE,IKA,1))/NS(IPHASE,I) DKY = (BKP(IPHASE,IKB,2)-BKP(IPHASE,IKA,2))/NS(IPHASE,I) DKZ = (BKP(IPHASE,IKB,3)-BKP(IPHASE,IKA,3))/NS(IPHASE,I) DK = DSQRT(DKX**2+DKY**2+DKZ**2) NSI = NS(IPHASE,I) IF (I.EQ.LINK(IPHASE)) NSI=NSI+1 DO J=1,NSI FKX = BKP(IPHASE,IKA,1)+(J-1)*DKX FKY = BKP(IPHASE,IKA,2)+(J-1)*DKY FKZ = BKP(IPHASE,IKA,3)+(J-1)*DKZ FKXA = FKX*AA(IPHASE)/PI FKYA = FKY*AA(IPHASE)/PI FKZA = FKZ*AA(IPHASE)/PI IF (J.EQ.1) THEN ISPEC = IKA ELSE ISPEC = 0 ENDIF IF (I.EQ.LINK(IPHASE).AND.J.EQ.NSI) ISPEC = IKB CALL DIAG_K (IPHASE,FKX,FKY,FKZ) WRITE (LP,'(3(F9.6,1X),I2,1X,F7.4,30(1X,F10.6))') A FKXA,FKYA,FKZA,ISPEC,RKTOT,(W(II),II=1,NPA(MYOBJ(IPHASE))*4) RKTOT = RKTOT+DK ENDDO ENDDO CLOSE(LP) ENDIF 10 CONTINUE PRINT *,' ' C Cohesive energy curves of object groups: PRINT *,'Cohesive energy curves of object groups: ' PRINT *,' ' DO IOBJ = 1,NOBJ OPEN (UNIT=LP, STATUS='UNKNOWN', FORM='FORMATTED', A FILE=NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1)//'.out') PRINT *,'Calculating cohesive energy for ', A NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1) PRINT *,' ' WRITE (*,'( A " Length(A) Volume(A^3) CPS(eV) ETB(eV) Ef(eV) ", A "Ebs(eV) Erep(eV) Ec(eV) Mulliken charges ->")') CURV_ERROR = 0.D0 ABS_CURV_ERROR = 0.D0 DO IPHASE = NPHASE_OBJ(IOBJ,1),NPHASE_OBJ(IOBJ,2) WRITE (*,'(F9.6,11(F11.6))') A AA(IPHASE), ABS(DET(IOBJ))*AA(IPHASE)**3/NPA(IOBJ), A CPS(IPHASE),ETB(IPHASE),EFERMI(IPHASE),EBS(IPHASE),EREP(IPHASE), A EC(IPHASE), (CHARGE(IPHASE,N),N=1,NPA(IOBJ)) WRITE (LP,'(F9.6,11(F11.6))') A AA(IPHASE), ABS(DET(IOBJ))*AA(IPHASE)**3/NPA(IOBJ), A CPS(IPHASE),ETB(IPHASE),EFERMI(IPHASE),EBS(IPHASE),EREP(IPHASE), A EC(IPHASE), (CHARGE(IPHASE,N),N=1,NPA(IOBJ)) IF (IPHASE.GT.NPHASE_OBJ(IOBJ,1)) THEN CC = ( ETB(IPHASE)-ETB(NPHASE_OBJ(IOBJ,1)) A -CPS(IPHASE)+CPS(NPHASE_OBJ(IOBJ,1)) ) / A (CPS(IPHASE)-CPS(NPHASE_OBJ(IOBJ,1))) / A (NPHASE_OBJ(IOBJ,2)-NPHASE_OBJ(IOBJ,1)) CURV_ERROR = CURV_ERROR + CC ABS_CURV_ERROR = ABS_CURV_ERROR + ABS(CC) ENDIF ENDDO CLOSE (LP) WRITE (*,'(" Relative ground error = ",f5.1, A "%, curv error = ",f5.1,"% (",f4.1,"%)",/)') A -(ETB(NPHASE_OBJ(IOBJ,1))-CPS(NPHASE_OBJ(IOBJ,1))) A /CPS(NPHASE_OBJ(IOBJ,1))*100, CURV_ERROR*100, A ABS_CURV_ERROR*100 ENDDO C Deformation mode energetics: DO IDEF = 1,NDEF PRINT *,'Deformation mode ', A NAME_DEF(IDEF)(1:INDEX(NAME_DEF(IDEF),' ')-1), A ' energetics:' PRINT *,' ' OPEN (UNIT=LP, STATUS='UNKNOWN', FORM='FORMATTED', A FILE=NAME_DEF(IDEF)(1:INDEX(NAME_DEF(IDEF),' ')-1)//'.out') DEF_ERROR = 0.D0 ABS_DEF_ERROR = 0.D0 DO I = 1,NELE(IDEF) IPHASE = IELE(IDEF,I) WRITE (*,'(6(F11.6))') CPS(IPHASE),ETB(IPHASE), A EFERMI(IPHASE),EBS(IPHASE),EREP(IPHASE),EC(IPHASE) WRITE (LP,'(6(F11.6))') CPS(IPHASE),ETB(IPHASE), A EFERMI(IPHASE),EBS(IPHASE),EREP(IPHASE),EC(IPHASE) IF (I.GT.1) THEN CC = ( ETB(IPHASE)-ETB(IELE(IDEF,1))- A CPS(IPHASE)+CPS(IELE(IDEF,1)) ) / A (CPS(IPHASE)-CPS(IELE(IDEF,1))) / A ( NELE(IDEF)-1 ) DEF_ERROR = DEF_ERROR + CC ABS_DEF_ERROR = ABS_DEF_ERROR + ABS(CC) ENDIF ENDDO CLOSE (LP) WRITE (*,'(" Relative error = ",f5.1,"% (",f4.1,"%)",/)') A DEF_ERROR*100, ABS_DEF_ERROR*100 ENDDO CLOSE(LP_PARA) STOP END SUBROUTINE INITIALIZE_STRUCTURE (IPHASE) C -------------------------------------------------------------- C ASSIGN ATOMIC POSITIONS AND CALCULATE THE COORDINATION NUMBERS C -------------------------------------------------------------- # include "fitall.fh" PARAMETER (MML=7) IOBJ = MYOBJ(IPHASE) C Step 1: unit cell in the middle IP = 0 DO N = 1,NPA(IOBJ) IP = IP+1 X(IPHASE,IP,1) = AA(IPHASE)*TAU(IOBJ,N,1) X(IPHASE,IP,2) = AA(IPHASE)*TAU(IOBJ,N,2) X(IPHASE,IP,3) = AA(IPHASE)*TAU(IOBJ,N,3) IPT(IPHASE,IP) = N ENDDO C Step 2: all those which can interact with center cell DO 10 II = -MML,MML DO 10 JJ = -MML,MML DO 10 KK = -MML,MML IF (II.EQ.0.AND.JJ.EQ.0.AND.KK.EQ.0) GOTO 10 DO 10 N = 1,NPA(IOBJ) IP = IP+1 IF (IP.GT.MMP) THEN PRINT *, 'initialize_structure(): MMP reached.' STOP ENDIF X(IPHASE,IP,1)=AA(IPHASE)*(II*A(IOBJ,1,1)+JJ*A(IOBJ,2,1)+ A KK*A(IOBJ,3,1)+ TAU(IOBJ,N,1)) X(IPHASE,IP,2)=AA(IPHASE)*(II*A(IOBJ,1,2)+JJ*A(IOBJ,2,2)+ A KK*A(IOBJ,3,2)+ TAU(IOBJ,N,2)) X(IPHASE,IP,3)=AA(IPHASE)*(II*A(IOBJ,1,3)+JJ*A(IOBJ,2,3)+ A KK*A(IOBJ,3,3)+ TAU(IOBJ,N,3)) IPT(IPHASE,IP) = N DO M = 1,NPA(IOBJ) IF ( (X(IPHASE,IP,1)-X(IPHASE,M,1))**2+ A (X(IPHASE,IP,2)-X(IPHASE,M,2))**2+ A (X(IPHASE,IP,3)-X(IPHASE,M,3))**2.LE. A RCUT(NTYPE(IOBJ,N),NTYPE(IOBJ,M),2)**2 ) GOTO 10 ENDDO IP = IP-1 10 CONTINUE C It happens that they are also those which can screen C core-outside interactions: NPP(IPHASE) = IP C Calculate the distances between atoms: DO I = 1,NPA(IOBJ) DO JP = 1,NPP(IPHASE) R(IPHASE,I,JP) = DSQRT((X(IPHASE,I,1)-X(IPHASE,JP,1))**2+ A (X(IPHASE,I,2)-X(IPHASE,JP,2))**2+ A (X(IPHASE,I,3)-X(IPHASE,JP,3))**2) ENDDO ENDDO C Calculate the effective coordination number according to (5) IF (ITOT.EQ.1) WRITE (*, '(" NPP = ", I3, A ", and effective coordination number of atoms are")') A NPP(IPHASE) DO I=1,NPA(IOBJ) G(IPHASE,I,1) = 0.D0 G(IPHASE,I,2) = 0.D0 I_TYPE = NTYPE(IOBJ,I) DO 30 JP=1,NPP(IPHASE) J_TYPE = NTYPE(IOBJ,IPT(IPHASE,JP)) ETAIJ = 0.D0 RIJ = R(IPHASE,I,JP) IF (RIJ.GT.RCUT(I_TYPE,J_TYPE,2).OR.I.EQ.JP) GOTO 30 DO 40 LP=1,NPP(IPHASE) L_TYPE=NTYPE(IOBJ,IPT(IPHASE,LP)) RIL = R(IPHASE,I,LP) RJL = DSQRT((X(IPHASE,JP,1)-X(IPHASE,LP,1))**2 A +(X(IPHASE,JP,2)-X(IPHASE,LP,2))**2 A +(X(IPHASE,JP,3)-X(IPHASE,LP,3))**2) IF (RIL.GT.RCUT(I_TYPE,L_TYPE,2).OR. A RJL.GT.RCUT(J_TYPE,L_TYPE,2).OR. A LP.EQ.I.OR.LP.EQ.JP) GOTO 40 C Wired-in screening formula for effective coordination number ETAIJ = ETAIJ+2.0*DEXP(-0.0478*((RIL+RJL)/RIJ)**7.16) 40 CONTINUE XX = DEXP(ETAIJ) CC = 1.D0-(XX-1./XX)/(XX+1./XX) C Cutoff smoother from CZ (June 21, 97): IF (RIJ.GT.RCUT(I_TYPE,J_TYPE,1)) CC = CC* A DCOS( PI/2 * (RIJ-RCUT(I_TYPE,J_TYPE,1)) / A (RCUT(I_TYPE,J_TYPE,2)-RCUT(I_TYPE,J_TYPE,1)) )**2 G(IPHASE,I,J_TYPE) = G(IPHASE,I,J_TYPE) + CC 30 CONTINUE IF (ITOT.EQ.1) WRITE (*, '(I2,3F14.10)') A I, G(IPHASE,I,1), G(IPHASE,I,2), G(IPHASE,I,1)+G(IPHASE,I,2) ENDDO IF (ITOT.EQ.1) PRINT *,' ' RETURN END DOUBLE PRECISION FUNCTION FE (GUESS) C ----------------------------------------- C GET THE ERROR OF THE GUESS PARAMETER SET C ----------------------------------------- # include "fitall.fh" DOUBLE PRECISION GUESS(NOP) C Feed GUESS parameters into the Generic Table: IP = 1 DO IOP = 17,25 A1(IOP) = GUESS(IP) A2(IOP) = GUESS(IP+1) A3(IOP) = GUESS(IP+2) A4(IOP) = GUESS(IP+3) IP = IP+4 ENDDO DO IOP = 2,7 DELTA(IOP,1,2) = GUESS(IP) IP = IP+1 ENDDO DO IOP = 9,14 DELTA(IOP,2,1) = GUESS(IP) IP = IP+1 ENDDO RIGID_SHIFT = GUESS(IP) S_SI = GUESS(IP+1) S_C = GUESS(IP+2) DO IOP = 1,25 D(IOP,1) = 1. D(IOP,2) = 1. F(IOP,1) = 0.0 F(IOP,2) = 0.0 ENDDO C Simple screening rules: DO IOP = 2,7 B1(IOP,2) = B1(IOP+7,2) B2(IOP,2) = B2(IOP+7,2) B3(IOP,2) = B3(IOP+7,2) C B1(IOP,2) = B1(IOP+7,2)**GUESS(IP+3) * B1(IOP,1)**(1-GUESS(IP+3)) C B2(IOP,2) = B2(IOP+7,2)**GUESS(IP+3) * B2(IOP,1)**(1-GUESS(IP+3)) C B3(IOP,2) = B3(IOP+7,2)**GUESS(IP+3) * B3(IOP,1)**(1-GUESS(IP+3)) ENDDO DO IOP = 9,14 B1(IOP,1) = B1(IOP-7,1) B2(IOP,1) = B2(IOP-7,1) B3(IOP,1) = B3(IOP-7,1) C B1(IOP,1) = B1(IOP-7,1)**GUESS(IP+3) * B1(IOP,2)**(1-GUESS(IP+3)) C B2(IOP,1) = B2(IOP-7,1)**GUESS(IP+3) * B2(IOP,2)**(1-GUESS(IP+3)) C B3(IOP,1) = B3(IOP-7,1)**GUESS(IP+3) * B3(IOP,2)**(1-GUESS(IP+3)) ENDDO DO IOP = 18,23 B1(IOP,1) = DSQRT(B1(IOP-16,1)*B1(IOP-9,1)) B2(IOP,1) = DSQRT(B2(IOP-16,1)*B2(IOP-9,1)) B3(IOP,1) = DSQRT(B3(IOP-16,1)*B3(IOP-9,1)) B1(IOP,2) = DSQRT(B1(IOP-16,2)*B1(IOP-9,2)) B2(IOP,2) = DSQRT(B2(IOP-16,2)*B2(IOP-9,2)) B3(IOP,2) = DSQRT(B3(IOP-16,2)*B3(IOP-9,2)) DELTA(IOP,1,1) = DELTA(IOP-16,1,1) DELTA(IOP,2,2) = DELTA(IOP-9,2,2) C DELTA(IOP,1,2) = DELTA(IOP-16,1,2) C DELTA(IOP,2,1) = DELTA(IOP-9,2,1) DELTA(IOP,1,2) = ( DELTA(IOP-16,1,2)+DELTA(IOP-9,2,1) ) / 2. A * GUESS(IP+3) DELTA(IOP,2,1) = DELTA(IOP,1,2) ENDDO C screening to Es,Ep(Si) = Es,Ep(C) B1(17,1) = B1(18,1) B2(17,1) = B2(18,1) B3(17,1) = B3(18,1) DELTA(17,1,1) = DELTA(18,1,1) DELTA(17,1,2) = DELTA(18,1,2) B1(17,2) = B1(18,2) B2(17,2) = B2(18,2) B3(17,2) = B3(18,2) DELTA(17,2,1) = DELTA(18,2,1) DELTA(17,2,2) = DELTA(18,2,2) C screening to Phi(Si->C) = Phi(C->Si) B1(24,1) = B1(23,1) B2(24,1) = B2(23,1) B3(24,1) = B3(23,1) DELTA(24,1,1) = DELTA(23,1,1) DELTA(24,1,2) = DELTA(23,1,2) B1(24,2) = B1(23,2) B2(24,2) = B2(23,2) B3(24,2) = B3(23,2) DELTA(24,2,1) = DELTA(23,2,1) DELTA(24,2,2) = DELTA(23,2,2) C screening to ps_sigma = sp_sigma B1(25,1) = B1(20,1) B2(25,1) = B2(20,1) B3(25,1) = B3(20,1) DELTA(25,1,1) = DELTA(20,1,1) DELTA(25,1,2) = DELTA(20,1,2) B1(25,2) = B1(20,2) B2(25,2) = B2(20,2) B3(25,2) = B3(20,2) DELTA(25,2,1) = DELTA(20,2,1) DELTA(25,2,2) = DELTA(20,2,2) IF (ITOT.EQ.MAXITER) WRAPPING_UP = .TRUE. COH_ERROR = 0.D0 IF (DO_EBS) BAND_ERROR = 0.D0 DO IPHASE = 1,NPHASE IOBJ = MYOBJ(IPHASE) IF ( CURV_WEIGHT(IOBJ).GT.0.OR. A (GROUND_WEIGHT(IOBJ).GT.0.AND. A IPHASE.EQ.NPHASE_OBJ(IOBJ,1)).OR. A WB(IPHASE).GT.0.OR.WRAPPING_UP) THEN C Construct the overlap matrix and calculate Erep CALL OVERLAP (IPHASE) IF (DO_EBS) CALL CALC_EBS(IPHASE) CALL CALC_EC (IPHASE) ETB(IPHASE) = EREP(IPHASE) + EBS(IPHASE) + EC(IPHASE) ENDIF IF (DO_EBS.AND.WB(IPHASE).GT.0) THEN C First, calculate the lowest gamma point energy: CALL DIAG_K (IPHASE,0.D0,0.D0,0.D0) EGAMMA = W(1) C Calculate the band structure: DO NK = 1,NBKP(IPHASE) CALL DIAG_K (IPHASE, BKP(IPHASE,NK,1), BKP(IPHASE,NK,2), A BKP(IPHASE,NK,3)) DO I = 1,NBA(IPHASE,NK) BAND_ERROR = BAND_ERROR + BWEIGHT(IPHASE,NK,I) * A ABS(W(NTH(IPHASE,NK,I))-EGAMMA-BTARGET(IPHASE,NK,I)) ENDDO ENDDO ENDIF ENDDO IF (FIX_EBS) DO_EBS = .FALSE. C Object group errors: DO IOBJ = 1,NOBJ DO IPHASE = NPHASE_OBJ(IOBJ,1),NPHASE_OBJ(IOBJ,2) IF (IPHASE.EQ.NPHASE_OBJ(IOBJ,1)) THEN COH_ERROR = COH_ERROR + GROUND_WEIGHT(IOBJ) * A ABS((ETB(IPHASE)-CPS(IPHASE))/CPS(IPHASE)) ELSE COH_ERROR = COH_ERROR + A CURV_WEIGHT(IOBJ)/(NPHASE_OBJ(IOBJ,2)-NPHASE_OBJ(IOBJ,1)) A * ABS( (ETB(IPHASE)-ETB(NPHASE_OBJ(IOBJ,1)) A -CPS(IPHASE)+CPS(NPHASE_OBJ(IOBJ,1)))/ A (CPS(IPHASE)-CPS(NPHASE_OBJ(IOBJ,1))) ) ENDIF ENDDO ENDDO C Deformation mode errors: DO IDEF = 1,NDEF DO I = 2,NELE(IDEF) COH_ERROR = COH_ERROR + DEF_WEIGHT(IDEF)/(NELE(IDEF)-1) A * ABS( (ETB(IELE(IDEF,I))-ETB(IELE(IDEF,1))- A CPS(IELE(IDEF,I))+CPS(IELE(IDEF,1)))/ A (CPS(IELE(IDEF,I))-CPS(IELE(IDEF,1))) ) ENDDO ENDDO FE = BAND_ERROR * BAND_ERROR_RATIO + A COH_ERROR * (1-BAND_ERROR_RATIO) WRITE (*,'(" step=",I6,", band=",F6.3," eV, coh=", F6.3, A ", fe=", F6.3)') ITOT, BAND_ERROR, COH_ERROR, FE IF (MOD(ITOT,IOFTEN).EQ.0) THEN C ------------------------------------------------------------------ C PRINT OUT THE PARAMETERS WRITE (LP_PARA,'("step=",I6,", band=",F6.3," eV, coh=", F6.3, A ", fe=", F6.3)') ITOT, BAND_ERROR, COH_ERROR, FE IP = 1 DO I=1,13 WRITE (LP_PARA,'(A70)') NAME_PARA(I) WRITE (LP_PARA,'(4F11.6,1X)') A BOUND(IP,2),BOUND(IP+1,2),BOUND(IP+2,2),BOUND(IP+3,2) WRITE (LP_PARA,'(4F11.6,1x)') A GUESS(IP), GUESS(IP+1), GUESS(IP+2), GUESS(IP+3) WRITE (LP_PARA,'(4F11.6,1x)') A BOUND(IP,1),BOUND(IP+1,1),BOUND(IP+2,1),BOUND(IP+3,1) IP = IP+4 ENDDO ENDIF C ------------------------------------------------------------------ IF (ITOT.EQ.MAXITER) THEN PRINT *,'Maximum Iteration =',MAXITER,' reached;' PRINT *,'Optimization stop.' PRINT *,' ' CALL WRAPUP ENDIF ITOT = ITOT+1 RETURN END SUBROUTINE OVERLAP (IPHASE) C --------------------------------------------------------------- C SUBROUTINE TO CALCULATE THE OVERLAP MATRIX TM(MMA,MMA,MAC,4,4) C AND DX(MMA,MMA,MAC,3) IN STACKED FORM. THE INPUTS ARE THE C REAL COORDINATES X(MMS,MMP,3) AND EFFECTIVE COORDINATION C NUMBER G(MMS,MMA,MMT). THE GENERIC PARAMETER TABLE MUST BE C FILLED OUT BEFOREHAND. EREP WILL ALSO BE EVALUATED. C ---------------------------------------------------------------- # include "fitall.fh" DIMENSION H(NTO), ETA(NTO) LOGICAL SWAP SWAP = .FALSE. IOBJ = MYOBJ(IPHASE) C On-site energy: ES0(1) = -2.071 + RIGID_SHIFT EP0(1) = 4.709 + RIGID_SHIFT ES0(2) = -6.041 EP0(2) = 1.024 DO I = 1,NPA(IOBJ) ES(I) = ES0(NTYPE(IOBJ,I)) EP(I) = EP0(NTYPE(IOBJ,I)) PHI(I,1) = 0.0 PHI(I,2) = 0.0 DO J = 1,NPA(IOBJ) IDX(I,J) = 0 ENDDO ENDDO C First fill in off-diagonal hopping integrals and C calculate influences to the on-site energies. DO 10 III = 1,NPA(IOBJ) I = III IP = I I_TYPE = NTYPE(IOBJ,I) DO 10 JJJP = 1,NPP(IPHASE) C We can document all the interactions by only C counting those that has j >= i IF (IPT(IPHASE,JJJP).LT.III.OR.JJJP.EQ.IP) GOTO 10 JP = JJJP J = IPT(IPHASE,JP) J_TYPE = NTYPE(IOBJ,J) RIJ = R(IPHASE,I,JP) IF (RIJ.GT.RCUT(I_TYPE,J_TYPE,2)) GOTO 10 IF (I_TYPE.EQ.1.AND.J_TYPE.EQ.1) THEN MIN = 2 IF (.NOT.DO_EBS) MIN = 7 MAX = 7 ELSEIF (I_TYPE.EQ.2.AND.J_TYPE.EQ.2) THEN MIN = 9 IF (.NOT.DO_EBS) MIN = 14 MAX = 14 ELSE MIN = 17 IF (.NOT.DO_EBS) MIN = 23 MAX = 25 IF (.NOT.DO_EBS) MIN = 24 C Always make i Si IF (I_TYPE.EQ.2) THEN JJ = J JJP = JP J = I JP = IP I = JJ IP = JJP I_TYPE = 1 J_TYPE = 2 SWAP = .TRUE. ENDIF ENDIF XJ = X(IPHASE,JP,1) - X(IPHASE,IP,1) YJ = X(IPHASE,JP,2) - X(IPHASE,IP,2) ZJ = X(IPHASE,JP,3) - X(IPHASE,IP,3) CX = XJ/RIJ CY = YJ/RIJ CZ = ZJ/RIJ C First, calculate the interaction strengths: DO IOP = MIN,MAX C Rescale the distances T0(1) = 4.000467 T0(2) = G0-T0(1) RR = RIJ * ABS( 1 + 0.5*DELTA(IOP,I_TYPE,1)*G(IPHASE,I,1)/G0 A + 0.5*DELTA(IOP,I_TYPE,2)*G(IPHASE,I,2)/G0 A - 0.5*DELTA(IOP,I_TYPE,J_TYPE)*T0(1)/G0 A - 0.5*DELTA(IOP,I_TYPE,I_TYPE)*T0(2)/G0 A + 0.5*DELTA(IOP,J_TYPE,1)*G(IPHASE,J,1)/G0 A + 0.5*DELTA(IOP,J_TYPE,2)*G(IPHASE,J,2)/G0 A - 0.5*DELTA(IOP,J_TYPE,I_TYPE)*T0(1)/G0 A - 0.5*DELTA(IOP,J_TYPE,J_TYPE)*T0(2)/G0 ) H(IOP) = A1(IOP)*RR**(-A2(IOP))*DEXP(-A3(IOP)*RR**A4(IOP)) IF (RIJ.GT.RCUT(I_TYPE,J_TYPE,1)) H(IOP) = H(IOP) * A DCOS( PI/2 * (RIJ-RCUT(I_TYPE,J_TYPE,1)) / A (RCUT(I_TYPE,J_TYPE,2)-RCUT(I_TYPE,J_TYPE,1)) )**2 ETA(IOP) = 0. ENDDO C Then, calculate the screening by looping over l's DO 20 LP = 1,NPP(IPHASE) IF (LP.EQ.IP.OR.LP.EQ.JP) GOTO 20 L_TYPE = NTYPE(IOBJ,IPT(IPHASE,LP)) XL = X(IPHASE,LP,1)-X(IPHASE,IP,1) YL = X(IPHASE,LP,2)-X(IPHASE,IP,2) ZL = X(IPHASE,LP,3)-X(IPHASE,IP,3) RIL2 = XL*XL + YL*YL + ZL*ZL RJL2 = (XJ-XL)*(XJ-XL) + (YJ-YL)*(YJ-YL) + (ZJ-ZL)*(ZJ-ZL) IF (RIL2.GT.RCUT(I_TYPE,L_TYPE,2)**2.OR. A RJL2.GT.RCUT(J_TYPE,L_TYPE,2)**2) GOTO 20 RIL = DSQRT(RIL2) RJL = DSQRT(RJL2) DO IOP = MIN,MAX C Distance to the additional force center between i-j RIF = DSQRT((XL-0.5*XJ)**2+(YL-0.5*YJ)**2+(ZL-0.5*ZJ)**2) CC = ( D(IOP,L_TYPE) * RIL + A 1./D(IOP,L_TYPE) * RJL + A F(IOP,L_TYPE) * RIF ) / RIJ ETA(IOP) = ETA(IOP) + B1(IOP,L_TYPE) * A DEXP( -B2(IOP,L_TYPE) * CC**B3(IOP,L_TYPE) ) ENDDO 20 CONTINUE C Screen the interactions DO IOP = MIN,MAX XX = DEXP(ETA(IOP)) H(IOP) = H(IOP)*2./(XX+1./XX)/XX C To avoid numerical problems: IF (H(IOP).GT.1/EPS.OR.H(IOP).LT.-1/EPS) H(IOP)=0. ENDDO IF (DO_EBS) THEN ID = IDX(I,J)+1 IF (ID.GT.MAC) THEN PRINT *, 'overlap(): MAC reached.' STOP ENDIF IDX(I,J) = ID DX(I,J,ID,1) = -XJ DX(I,J,ID,2) = -YJ DX(I,J,ID,3) = -ZJ ENDIF IF (I_TYPE.EQ.J_TYPE) THEN IF (I.EQ.J) THEN IF (DO_EBS) THEN ES(I) = ES(I) + H(MIN) EP(I) = EP(I) + H(MIN) ENDIF PHI(I,J_TYPE) = PHI(I,J_TYPE) + H(MAX) ELSE IF (DO_EBS) THEN ES(I) = ES(I) + H(MIN) EP(I) = EP(I) + H(MIN) ES(J) = ES(J) + H(MIN) EP(J) = EP(J) + H(MIN) ENDIF PHI(I,J_TYPE) = PHI(I,J_TYPE) + H(MAX) PHI(J,I_TYPE) = PHI(J,I_TYPE) + H(MAX) ENDIF IF (DO_EBS) THEN TM(I,J,ID,1,1) = H(MIN+1) TM(I,J,ID,1,2) = H(MIN+2)*CX TM(I,J,ID,1,3) = H(MIN+2)*CY TM(I,J,ID,1,4) = H(MIN+2)*CZ TM(I,J,ID,2,1) = -TM(I,J,ID,1,2) TM(I,J,ID,3,1) = -TM(I,J,ID,1,3) TM(I,J,ID,4,1) = -TM(I,J,ID,1,4) C px-px TM(I,J,ID,2,2) = CX*CX*H(MIN+3)+(1-CX*CX)*H(MIN+4) C px-py TM(I,J,ID,2,3) = CX*CY*(H(MIN+3)-H(MIN+4)) TM(I,J,ID,3,2) = TM(I,J,ID,2,3) C px-pz TM(I,J,ID,2,4) = CX*CZ*(H(MIN+3)-H(MIN+4)) TM(I,J,ID,4,2) = TM(I,J,ID,2,4) C py-py TM(I,J,ID,3,3) = CY*CY*H(MIN+3)+(1-CY*CY)*H(MIN+4) C py-pz TM(I,J,ID,3,4) = CY*CZ*(H(MIN+3)-H(MIN+4)) TM(I,J,ID,4,3) = TM(I,J,ID,3,4) C pz-pz TM(I,J,ID,4,4) = CZ*CZ*H(MIN+3)+(1-CZ*CZ)*H(MIN+4) ENDIF ELSE IF (DO_EBS) THEN ES(I) = ES(I) + H(MIN) EP(I) = EP(I) + H(MIN) ES(J) = ES(J) + H(MIN+1) EP(J) = EP(J) + H(MIN+1) PHI(I,J_TYPE) = PHI(I,J_TYPE) + H(MAX-2) PHI(J,I_TYPE) = PHI(J,I_TYPE) + H(MAX-1) TM(I,J,ID,1,1) = H(MIN+2) TM(I,J,ID,1,2) = H(MIN+3) * CX TM(I,J,ID,1,3) = H(MIN+3) * CY TM(I,J,ID,1,4) = H(MIN+3) * CZ TM(I,J,ID,2,1) = -H(MAX) * CX TM(I,J,ID,3,1) = -H(MAX) * CY TM(I,J,ID,4,1) = -H(MAX) * CZ C px-px TM(I,J,ID,2,2) = CX*CX*H(MIN+4)+(1-CX*CX)*H(MIN+5) C px-py TM(I,J,ID,2,3) = CX*CY*(H(MIN+4)-H(MIN+5)) TM(I,J,ID,3,2) = TM(I,J,ID,2,3) C px-pz TM(I,J,ID,2,4) = CX*CZ*(H(MIN+4)-H(MIN+5)) TM(I,J,ID,4,2) = TM(I,J,ID,2,4) C py-py TM(I,J,ID,3,3) = CY*CY*H(MIN+4)+(1-CY*CY)*H(MIN+5) C py-pz TM(I,J,ID,3,4) = CY*CZ*(H(MIN+4)-H(MIN+5)) TM(I,J,ID,4,3) = TM(I,J,ID,3,4) C pz-pz TM(I,J,ID,4,4) = CZ*CZ*H(MIN+4)+(1-CZ*CZ)*H(MIN+5) ELSE PHI(I,J_TYPE) = PHI(I,J_TYPE) + H(MAX-1) PHI(J,I_TYPE) = PHI(J,I_TYPE) + H(MAX) ENDIF ENDIF C Change back if swapped: IF (SWAP) THEN JJ = J JJP = JP J = I JP = IP I = JJ IP = JJP I_TYPE = 2 J_TYPE = 1 SWAP = .FALSE. ENDIF 10 CONTINUE IF (DO_EBS) THEN C Fill in the on-site terms DO I=1,NPA(IOBJ) ID = IDX(I,I)+1 IF (ID.GT.MAC) THEN PRINT *, 'overlap(): MAC reached.' STOP ENDIF IDX(I,I) = ID DO J=1,4 DO K=1,4 TM(I,I,ID,J,K)=0. ENDDO ENDDO TM(I,I,ID,1,1)=ES(I) TM(I,I,ID,2,2)=EP(I) TM(I,I,ID,3,3)=EP(I) TM(I,I,ID,4,4)=EP(I) DX(I,I,ID,1)=0. DX(I,I,ID,2)=0. DX(I,I,ID,3)=0. ENDDO ENDIF X0 = 0. X1 = 0. X2 = 0. X3 = 0. X4 = 0. Y0 = 0. Y1 = 0. Y2 = 0. Y3 = 0. Y4 = 0. DO K = 1,NPA(IOBJ) IF (NTYPE(IOBJ,K).EQ.1) THEN X0 = X0 + 1. X1 = X1 + (PHI(K,1)/S_SI + PHI(K,2)) X2 = X2 + (PHI(K,1)/S_SI + PHI(K,2))**2. X3 = X3 + (PHI(K,1)/S_SI + PHI(K,2))**3. X4 = X4 + (PHI(K,1)/S_SI + PHI(K,2))**4. ELSE Y0 = Y0 + 1. Y1 = Y1 + (PHI(K,2)/S_C + PHI(K,1)) Y2 = Y2 + (PHI(K,2)/S_C + PHI(K,1))**2. Y3 = Y3 + (PHI(K,2)/S_C + PHI(K,1))**3. Y4 = Y4 + (PHI(K,2)/S_C + PHI(K,1))**4. ENDIF ENDDO C From pure Si fit: FOFSI = (c0-4*rigid_shift)*x0 + c1*x1*s_si + A c2*x2*s_si**2. + c3*x3*s_si**3. + c4*x4*s_si**4. C From pure C fit: FOFC = d0*y0 + d1*y1*s_c + d2*y2*s_c**2. + A d3*y3*s_c**3. + d4*y4*s_c**4. EREP(IPHASE) = (FOFSI+FOFC)/NPA(IOBJ) RETURN END SUBROUTINE DIAG_K (IPHASE,FKX,FKY,FKZ) # include "fitall.fh" C --------------------------------------------------------------- C SUBROUTINE TO DIAGONALIZE THE HAMILTONIAN FOR A SPECIFIC K C THE INPUT ARE NA(NUMBER OF ATOMS), AND FKX,FKY,FKZ. THE OVERLAP C MATRIX TM(MMA,MMA,MAC,4,4) IN STACKED FORM MUST BE READY. C THE OUTPUT EIGENVALUES ARE WRITTEN IN W(MMD), WHICH IS ORDERD. C ---------------------------------------------------------------- COMPLEX*16 TMAT(MMD,MMD),FACTOR,AP(MMD*(MMD+1)/2),ADD,ZAUZ(2*MMD) DOUBLE PRECISION AUZ(6*MMD) NA = NPA(MYOBJ(IPHASE)) NOB = 4*NA DO II = 1,NOB DO JJ = II,NOB TMAT(II,JJ) = DCMPLX(0.,0.) ENDDO ENDDO DO I = 1,NA II = 4*(I-1) DO J = 1,NA JJ = 4*(J-1) DO ID = 1,IDX(I,J) PHASE = FKX*DX(I,J,ID,1)+FKY*DX(I,J,ID,2)+FKZ*DX(I,J,ID,3) FACTOR = DCMPLX (DCOS(PHASE),DSIN(PHASE)) DO IP=1,4 DO JP=1,4 ADD = TM(I,J,ID,IP,JP)*FACTOR C Because we may have (Si-C) exchanged cases: IF (I.LE.J) THEN TMAT(II+IP,JJ+JP)=TMAT(II+IP,JJ+JP)+ADD ELSE TMAT(JJ+JP,II+IP)=TMAT(JJ+JP,II+IP)+CONJG(ADD) ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO LENAI = 0 DO JJ = 1,NOB DO II = 1,JJ LENAI = LENAI+1 AP(LENAI) = TMAT(II,JJ) ENDDO ENDDO C LAPACK driver INFO = 0 CALL ZHPEV('V', 'U', NOB, AP, W, Z, MMD, ZAUZ, AUZ, INFO) RETURN END SUBROUTINE SEA (IPHASE,EF,FVALUE,DERIVATIVE) C -------------------------------------------------------------- C CALCULATE THE DOS INTEGRAL AND FIRST DERIVATIVE ON THE LEFT C OF CUTOFF ENERGY EF, FOR A GIVEN DISCRETE SPECTRUM STORED IN C ELIST(MMD*MZK) WITH GAUSSIAN BROADENING EDEL AND NORMALIZED C WEIGHT ZKP(MZK,4) SUCH THAT ALL ORBITALS ADD UP TO 1. C -------------------------------------------------------------- # include "fitall.fh" PARAMETER (CUTOFF=30.) IOBJ = MYOBJ(IPHASE) NORB = NZKP(IOBJ) * NPA(IOBJ) * 4 FVALUE = 0.D0 DERIVATIVE = 0.D0 DO I=1,NORB DE = (EF - ELIST(I)) / EDEL IF (DE.GT.CUTOFF) THEN EWEIGHT(I) = ZKP(IOBJ,(I-1)/(4*NPA(IOBJ))+1,4) ELSEIF (DE.GT.-CUTOFF) THEN EWEIGHT(I) = ZKP(IOBJ,(I-1)/(4*NPA(IOBJ))+1,4) * A (1. + DERF(DE/DSQRT(2.D0))) / 2. DERIVATIVE = DERIVATIVE + ZKP(IOBJ,(I-1)/(4*NPA(IOBJ))+1,4) A * DEXP(-DE**2/2.D0) / DSQRT(2*PI) / EDEL ELSE EWEIGHT(I) = 0. ENDIF FVALUE = FVALUE + EWEIGHT(I) ENDDO IF (DERIVATIVE.EQ.0.) DERIVATIVE = 1000000. RETURN END SUBROUTINE CALC_EBS (IPHASE) C ---------------------------------------------------------- C 1. CALCULATE THE FERMI ENERGY. C 2. CALCULATE THE ELECTRONIC BAND STRUCTURE ENERGY EBS. C 3. CALCULATE THE PARTIAL CHARGE ON EACH ATOM. C ---------------------------------------------------------- # include "fitall.fh" DIMENSION EORDER(MMD*MZK) IOBJ = MYOBJ(IPHASE) C Sample the IBZ k-points: IE = 0 DO NK = 1,NZKP(IOBJ) FKX = 2 * PI / AA(IPHASE) * ( ZKP(IOBJ,NK,1)*P(IOBJ,1,1) A + ZKP(IOBJ,NK,2)*P(IOBJ,1,2) + ZKP(IOBJ,NK,3)*P(IOBJ,1,3) ) FKY = 2 * PI / AA(IPHASE) * ( ZKP(IOBJ,NK,1)*P(IOBJ,2,1) A + ZKP(IOBJ,NK,2)*P(IOBJ,2,2) + ZKP(IOBJ,NK,3)*P(IOBJ,2,3) ) FKZ = 2 * PI / AA(IPHASE) * ( ZKP(IOBJ,NK,1)*P(IOBJ,3,1) A + ZKP(IOBJ,NK,2)*P(IOBJ,3,2) + ZKP(IOBJ,NK,3)*P(IOBJ,3,3) ) CALL DIAG_K (IPHASE,FKX,FKY,FKZ) DO I = 1, 4*NPA(IOBJ) IE = IE + 1 ELIST(IE) = W(I) DO J = 1, 4*NPA(IOBJ) ZLIST(J,IE) = Z(J,I) ENDDO ENDDO ENDDO EDEL = 0.05 * DSQRT(2.D0) NORB = 4 * NPA(IOBJ) * NZKP(IOBJ) C Order the eigenvalues. DO I = 1, NORB EORDER(I) = ELIST(I) ENDDO DO I = 1, NORB-1 DO J = NORB, I+1, -1 IF (EORDER(J) .LT. EORDER(J-1)) THEN EE = EORDER(J-1) EORDER(J-1) = EORDER(J) EORDER(J) = EE ENDIF ENDDO ENDDO C Bracket the Fermi energy between EORDER(IL) C and EORDER(IL+1) IL = NORB/2 EL = EORDER(IL) IH = IL+1 EH = EORDER(IH) CALL SEA (IPHASE, EL, FL, DERIVATIVE) CALL SEA (IPHASE, EH, FH, DERIVATIVE) 10 IF (FH.GT.0.5.AND.FL.GT.0.5) THEN IL = IL-1 IH = IH-1 EH = EL EL = EORDER(IL) FH = FL CALL SEA (IPHASE, EL, FL, DERIVATIVE) GOTO 10 ELSEIF (FH.LT.0.5.AND.FL.LT.0.5) THEN IL = IL+1 IH = IH+1 EL = EH EH = EORDER(IH) FL = FH CALL SEA (IPHASE, EH, FH, DERIVATIVE) GOTO 10 ENDIF C Guess the fermi energy from the low end where C derivative is big. EF = EL C Newton's method to determine Fermi energy. ITER = 0 20 CALL SEA (IPHASE, EF, FVALUE, DERIVATIVE) IF (ABS(FVALUE-0.5).LT.EPS) GOTO 30 EF1 = EF-(FVALUE-0.5)/DERIVATIVE IF (EF1.GT.EH) THEN EF1 = (EF + EH)/2. ELSEIF (EF1.LT.EL) THEN EF1 = (EF + EL)/2. ENDIF IF (FVALUE.GT.0.5) THEN EH = EF ELSE EL = EF ENDIF EF = EF1 ITER = ITER + 1 IF (ITER.LE.50) GOTO 20 PRINT *, A 'calc_ebs(): Newtons method out of maximum iteration.' STOP 30 EFERMI(IPHASE) = EF C Calculate the band structure energy. EBS(IPHASE) = 0. DO I = 1,NORB DE = (EF - ELIST(I)) / EDEL EBS(IPHASE) = EBS(IPHASE) + ZKP(IOBJ,(I-1)/(4*NPA(IOBJ))+1,4) A * ( ELIST(I) * (1 + DERF(DE/DSQRT(2.d0))) / 2. - A DEXP(-DE**2/2) / DSQRT(2*PI) * EDEL ) ENDDO C Energy per atom EBS(IPHASE) = 2 * EBS(IPHASE) * 4 C Calculate the charge occupation on each atomic site by eigenstates: DO IP = 1,NPA(IOBJ) CHARGE (IPHASE,IP) = 0. ENDDO DO I = 1,4*NPA(IOBJ) IP = (I-1)/4+1 DO J = 1,NORB CHARGE (IPHASE,IP) = CHARGE (IPHASE,IP) + A (REAL(ZLIST(I,J))**2 + IMAG(ZLIST(I,J))**2) * EWEIGHT(J) ENDDO ENDDO DO IP = 1,NPA(IOBJ) CHARGE (IPHASE,IP) = CHARGE (IPHASE,IP) * 2 * NPA(IOBJ) * 4 ENDDO RETURN END SUBROUTINE CALC_EC (IPHASE) C ----------------------------------------------------- C CALCULATE THE COULOMB ENERGY USING EWALD SUMMATION C ----------------------------------------------------- # include "fitall.fh" DIMENSION CFX(MMA),CFY(MMA),CFZ(MMA),S1(MMA),S2(MMA), A S3(MMA),Q(MMA), CSTRESS(3,3),HH(3,3),T(MMA,3) EXTERNAL INIT_EWALD, EWALD, EXIT_EWALD, EWALD_DYNAMICAL EC(IPHASE) = 0.D0 RETURN END C ======================================================================== C Notes: C C 1. CPS cohesive energy data should be store in Coh_Database/3C.cps C etc, the format should be: C lengthscale in A volume/atom(A^3) energy/atom(eV). C C Here I am using the data which is the smooth universal binding C curve fitted from raw LDA data. CZ says that since the raw LDA data C points come with some noise, it's better to use the smooth curve as the C database. I agree with him because I remembered in fitting the curve C a small change of parameters such as the domain width can induce big C error in the bulk modulus. C ========================================================================