C ------------------------------------------------------------
C PROGRAM FitAll V3.0
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 SYSTEMS
C 4. OPTIONAL U-TERM SELF CONSISTENT CHARGE ITERATION
C 5. COULOMB POTENTIAL AND FIRST-NEIGHBOR SOFTENING
C
C SEP.4 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
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IPHASE,I,J,N,NK,IP,IOP,IOBJ,NDATA,IDEF,MINIMIZER,
A NDIV,IERR,MAIN_EXIT_CODE
DOUBLE PRECISION CC,SUPPORT_BOUND,D11,D22,D33,D12,D23,D31,
A D13,D21,D32,TOT_COH_WEIGHT,TOT_BAND_WEIGHT,AA_LDA,
A C_INITIAL_CHARGE,SUMWEIGHT,THIS_BAND_TOTAL,DEL,
A ANNEAL_RATE,EMIN,LOWER_BOUND(NOP),UPPER_BOUND(NOP),
A GUESS(NOP),FE
EXTERNAL ASA_MAIN,FE
ITOT = 1
WRAPPING_UP = .FALSE.
NAME_ATOM(1) = 'Si'
NAME_ATOM(2) = 'C'
C ----------------------------------------------------------------
C These are going to be optimized:
READ (*,
A '("Fitting Parameters (Upper_bound Guess Lower_bound):")')
IP = 1
DO I = 1,16
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 (*,'(/,"Incremental Screening Cutoff:")')
READ *, CUT_INC(1), CUT_INC(2)
READ (*,'(/,"Screening Cutoff:")')
READ *, CUT_ETA(1), CUT_ETA(2)
READ (*,'(/,"Activate real-space cutoffs (y/n): "
A ,A70)') BUF
REAL_SPACE_CUTOFF = BUF.EQ.'y'.OR.BUF.EQ.'Y'
READ (*,'("They are (Si-Si, C-C, Si-C):")')
READ *,RCUT(1,1,1),RCUT(1,1,2),RCUT(2,2,1),RCUT(2,2,2),
A RCUT(1,2,1),RCUT(1,2,2)
RCUT(2,1,1) = RCUT(1,2,1)
RCUT(2,1,2) = RCUT(1,2,2)
READ (*,'(/,"Charge iteration tolerance: ",F)') CHARGE_TOLERANCE
READ (*,'("Use charges from last step (y/n): ",A70)') BUF
USE_LAST_STEP_CHARGES = BUF.EQ.'y'.OR.BUF.EQ.'Y'
READ (*,'("Charge transfer protection bound, exponent:")')
READ *, SUPPORT_BOUND, SUPPORT_EXPONENT
C Make the exponent even:
SUPPORT_EXPONENT = NINT(SUPPORT_EXPONENT/2)*2
C At DQ = support_bound, there is an elevation of
C potential by 1 volt.
SUPPORT(1) = 1./SUPPORT_BOUND**(SUPPORT_EXPONENT-1)
SUPPORT(2) = SUPPORT(1)
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)
C Count how many Si and C atoms in each cell:
NATOM(IOBJ,1) = 0
NATOM(IOBJ,2) = 0
DO N = 1,NPA(IOBJ)
NATOM(IOBJ,NTYPE(IOBJ,N)) = NATOM(IOBJ,NTYPE(IOBJ,N)) + 1
ENDDO
IF (NATOM(IOBJ,1)+NATOM(IOBJ,2).NE.NPA(IOBJ)) THEN
PRINT *,'fitall: Si/C checksum wrong.'
STOP
ENDIF
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.",
A " Initial C charge")')
DO 20 IPHASE = NPHASE_OBJ(IOBJ,1),NPHASE_OBJ(IOBJ,2)
READ (*,'(F16,A16,F16,F)') AA(IPHASE),
A NAME_BAND(IPHASE), WB(IPHASE), C_INITIAL_CHARGE
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(IPHASE), 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 AA(IPHASE) = AA_LDA
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 charges:
CALL INITIALIZE_CHARGES (IPHASE,C_INITIAL_CHARGE)
C Initialize the structure:
CALL INITIALIZE_STRUCTURE (IPHASE)
20 CONTINUE
PRINT *,' '
C Check if all volumes in this object group
C have the same class structure: they should.
DO IPHASE = NPHASE_OBJ(IOBJ,1),NPHASE_OBJ(IOBJ,2)
DO N = 1,NPA(IOBJ)
IF (MYCLASS(IPHASE,N).NE.MYCLASS(NPHASE_OBJ(IOBJ,1),N)) THEN
PRINT *, 'fitall: different class structure in a group.'
STOP
ENDIF
ENDDO
ENDDO
C Calculate the dimensionless bare Coulomb
C potential coefficients for this group:
DO N = 1,NPA(IOBJ)
S1(N) = P(IOBJ,1,1)*TAU(IOBJ,N,1) +
A P(IOBJ,2,1)*TAU(IOBJ,N,2) +
A P(IOBJ,3,1)*TAU(IOBJ,N,3)
S2(N) = P(IOBJ,1,2)*TAU(IOBJ,N,1) +
A P(IOBJ,2,2)*TAU(IOBJ,N,2) +
A P(IOBJ,3,2)*TAU(IOBJ,N,3)
S3(N) = P(IOBJ,1,3)*TAU(IOBJ,N,1) +
A P(IOBJ,2,3)*TAU(IOBJ,N,2) +
A P(IOBJ,3,3)*TAU(IOBJ,N,3)
IF (S1(N).GT.0.5) S1(N) = S1(N) - 1.
IF (S1(N).LT.-0.5) S1(N) = S1(N) + 1.
IF (S2(N).GT.0.5) S2(N) = S2(N) - 1.
IF (S2(N).LT.-0.5) S2(N) = S2(N) + 1.
IF (S3(N).GT.0.5) S3(N) = S3(N) - 1.
IF (S3(N).LT.-0.5) S3(N) = S3(N) + 1.
ENDDO
CALL SIMPLE_EWALD (A(IOBJ,1,1),A(IOBJ,1,2),A(IOBJ,1,3),
A A(IOBJ,2,1),A(IOBJ,2,2),A(IOBJ,2,3),
A A(IOBJ,3,1),A(IOBJ,3,2),A(IOBJ,3,3),
A NPA(IOBJ),S1,S2,S3,1D-10,POTE_MATRIX)
DO N = 1,NPA(IOBJ)*NPA(IOBJ)
OBJ_POTE_MATRIX(IOBJ,N) = POTE_MATRIX(N) * ENERGY_TO_EV
ENDDO
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 the 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 first:
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 ---------------------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IPHASE,I,J,IKA,IKB,NSI,ISPEC,II,IOBJ,N,IDEF
DOUBLE PRECISION CC,DKX,DKY,DKZ,DK,FKX,FKY,FKZ,CURV_ERROR,
A ABS_CURV_ERROR,DEF_ERROR,ABS_DEF_ERROR
PRINT *,'Starts property calculations... '
PRINT *,' '
DO 10 IPHASE = 1,NPHASE
C Plot out the band structure using optimized parameters
C and converged Mulliken charges from last FE
IF (NAME_BAND(IPHASE)(1:INDEX(NAME_BAND(IPHASE),' ')-1).NE.'NULL')
A 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,1)
OPEN (UNIT=LP,STATUS='UNKNOWN',FORM='FORMATTED',FILE=BUF)
CC = 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
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 FKX*AA(IPHASE)/PI,FKY*AA(IPHASE)/PI,FKZ*AA(IPHASE)/PI,
A ISPEC,CC,(W(II),II=1,NPA(MYOBJ(IPHASE))*4)
CC = CC+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) Eu(eV) Mulliken charges ->")')
CURV_ERROR = 0.D0
ABS_CURV_ERROR = 0.D0
DO IPHASE = NPHASE_OBJ(IOBJ,1),NPHASE_OBJ(IOBJ,2)
C VOLUME_LDA can be any identifier on the 2nd column of
C the .cps file, not just "volume" -- like the bond length.
WRITE (*,'(F9.6,11(F11.6))')
A AA(IPHASE), VOLUME_LDA(IPHASE), CPS(IPHASE),
A ETB(IPHASE), EFERMI(IPHASE), EBS(IPHASE), EREP(IPHASE),
A EU(IPHASE), (CHARGE(IPHASE,N),N=1,NPA(IOBJ))
WRITE (LP,'(F9.6,11(F11.6))')
A AA(IPHASE), VOLUME_LDA(IPHASE), CPS(IPHASE),
A ETB(IPHASE), EFERMI(IPHASE), EBS(IPHASE), EREP(IPHASE),
A EU(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)
C Negative means lower ground state or softer curvature
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), EU(IPHASE)
WRITE (LP,'(6(F11.6))') CPS(IPHASE), ETB(IPHASE),
A EFERMI(IPHASE), EBS(IPHASE), EREP(IPHASE), EU(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_CHARGES(IPHASE, C_INITIAL_CHARGE)
C Initialize charges for self-consistent iteration:
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IPHASE,IOBJ,N
DOUBLE PRECISION C_INITIAL_CHARGE,BB,CC,SUMCHARGE,
A SIMPLE_METHOD_CHARGE(MMT)
IOBJ = MYOBJ(IPHASE)
IF (C_INITIAL_CHARGE.GT.0) THEN
SIMPLE_METHOD_CHARGE(2) = C_INITIAL_CHARGE
SIMPLE_METHOD_CHARGE(1) = (4.*NPA(IOBJ) -
A SIMPLE_METHOD_CHARGE(2)*NATOM(IOBJ,2)) / NATOM(IOBJ,1)
DO N = 1,NPA(IOBJ)
CHARGE_INI(IPHASE,N) = SIMPLE_METHOD_CHARGE(NTYPE(IOBJ,N))
ENDDO
ELSE
C Load charges from .out file
BUF = NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1)//'.out'
PRINT *, 'Loading initial charges from ', BUF(1:INDEX(BUF,' ')-1)
OPEN (UNIT=LP, STATUS='UNKNOWN', FORM='FORMATTED', FILE=BUF)
C read on until error happens: a little brutal
10 READ (LP,'(F9.6,11(F11.6))')
A BB,CC,CC,CC,CC,CC,CC,CC,(CHARGE_INI(IPHASE,N),N=1,NPA(IOBJ))
IF (ABS(BB-AA(IPHASE)).LT.EPS) GOTO 20
GOTO 10
20 CLOSE (LP)
ENDIF
C Renormalize the charges:
SUMCHARGE = 0.
DO N = 1,NPA(IOBJ)
SUMCHARGE = SUMCHARGE + CHARGE_INI(IPHASE,N)
ENDDO
DO N = 1,NPA(IOBJ)
CHARGE_INI(IPHASE,N) = CHARGE_INI(IPHASE,N) -
A (SUMCHARGE/NPA(IOBJ) - 4.)
ENDDO
RETURN
END
SUBROUTINE INITIALIZE_STRUCTURE (IPHASE)
C ------------------------------------------------------
C CALCULATE THE COORDINATION NUMBERS AND DETERMINE THE
C INTERACTION AND SCREENING LIST, EQUIVALENT CLASSES
C ------------------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER MML,MMH,MMP
PARAMETER (MML=2,MMH=(2*MML+1)**3,MMP=MMH*MMA)
INTEGER IPHASE,IOBJ,IP,N,II,JJ,KK,I,I_TYPE,J,J_TYPE,
A JP,L_TYPE
DOUBLE PRECISION CC,ETAIJ,X(3,MMP)
IOBJ = MYOBJ(IPHASE)
NCLASS(IPHASE) = 0
C Calculate the G(r) if it is the first volume:
IF (IPHASE.EQ.NPHASE_OBJ(IOBJ,1)) OPEN
A (UNIT=LP_GR,STATUS='UNKNOWN',FORM='FORMATTED', FILE='GR/'//
A NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1)//'.gr')
C Step 1: unit cell in the middle
IP = 0
DO N = 1,NPA(IOBJ)
IP = IP+1
X(1,IP) = AA(IPHASE) * TAU(IOBJ,N,1)
X(2,IP) = AA(IPHASE) * TAU(IOBJ,N,2)
X(3,IP) = AA(IPHASE) * TAU(IOBJ,N,3)
ENDDO
C Step 2: those which enclose the 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
X(1,IP) = AA(IPHASE)*( II*A(IOBJ,1,1)+JJ*A(IOBJ,2,1)+
A KK*A(IOBJ,3,1)+ TAU(IOBJ,N,1) )
X(2,IP) = AA(IPHASE)*( II*A(IOBJ,1,2)+JJ*A(IOBJ,2,2)+
A KK*A(IOBJ,3,2)+ TAU(IOBJ,N,2) )
X(3,IP) = AA(IPHASE)*( II*A(IOBJ,1,3)+JJ*A(IOBJ,2,3)+
A KK*A(IOBJ,3,3)+ TAU(IOBJ,N,3) )
10 CONTINUE
PRINT *,
A '# Class Si-coord. C-coord. Tot-coord. Initial Q'
NIT(IPHASE) = 0
C Two chemically distinct types of coordination numbers:
DO I = 1,NPA(IOBJ)
G(IPHASE,I,1) = 0.D0
G(IPHASE,I,2) = 0.D0
ENDDO
DO 20 I = 1,NPA(IOBJ)
I_TYPE = NTYPE(IOBJ,I)
C We can enlist all the interactions by
C just recording those with j >= i:
DO 40 J = I,NPA(IOBJ)
J_TYPE = NTYPE(IOBJ,J)
DO 40 JP = J,MMH*NPA(IOBJ),NPA(IOBJ)
C Exclude self-interaction:
IF (JP.EQ.I) GOTO 40
NIT(IPHASE) = NIT(IPHASE) + 1
IF (NIT(IPHASE).GT.MMI) THEN
PRINT *, 'initialize_structure(): MMI reached.'
STOP
ENDIF
IT(1,NIT(IPHASE),IPHASE) = I
IT(2,NIT(IPHASE),IPHASE) = J
C IT(3,.) is 1 + the total number of screening
C atoms for this interaction
IT(3,NIT(IPHASE),IPHASE) = 1
C First record is i-j relative coordinate and distance:
R(1,1,NIT(IPHASE),IPHASE) = X(1,JP)-X(1,I)
R(2,1,NIT(IPHASE),IPHASE) = X(2,JP)-X(2,I)
R(3,1,NIT(IPHASE),IPHASE) = X(3,JP)-X(3,I)
R(4,1,NIT(IPHASE),IPHASE) = DSQRT( R(1,1,NIT(IPHASE),IPHASE)**2
A + R(2,1,NIT(IPHASE),IPHASE)**2 + R(3,1,NIT(IPHASE),IPHASE)**2 )
C ***** for C-C interaction ******
if (real_space_cutoff
C A .and.i_type.eq.j_type
A .and.i_type.eq.2.and.j_type.eq.2
A .and.R(4,1,NIT(IPHASE),IPHASE).gt.rcut(i_type,j_type,2)) goto 45
C ***** for C-C interaction ******
C For Si-C interaction, always make the first record Si:
IF (J_TYPE.EQ.1.AND.I_TYPE.EQ.2) THEN
IT(1,NIT(IPHASE),IPHASE) = J
IT(2,NIT(IPHASE),IPHASE) = I
R(1,1,NIT(IPHASE),IPHASE) = -R(1,1,NIT(IPHASE),IPHASE)
R(2,1,NIT(IPHASE),IPHASE) = -R(2,1,NIT(IPHASE),IPHASE)
R(3,1,NIT(IPHASE),IPHASE) = -R(3,1,NIT(IPHASE),IPHASE)
ENDIF
C Starts counting screening atoms:
ETAIJ = 0.
DO 50 LP = 1,MMH*NPA(IOBJ)
IF (LP.EQ.I.OR.LP.EQ.JP) GOTO 50
C Record number + 1
IT(3,NIT(IPHASE),IPHASE) = IT(3,NIT(IPHASE),IPHASE) + 1
IF (IT(3,NIT(IPHASE),IPHASE)-1.GT.MSC) THEN
PRINT *, 'initialize_structure(): MSC reached.'
STOP
ENDIF
C Starting from IT(4,.) are screening atom indices:
IT(IT(3,NIT(IPHASE),IPHASE)+2,NIT(IPHASE),IPHASE) =
A MOD(LP-1,NPA(IOBJ)) + 1
L_TYPE = NTYPE( IOBJ,IT(IT(3,NIT(IPHASE),IPHASE)+2,
A NIT(IPHASE),IPHASE) )
C RIL
R(1,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) =
A DSQRT( (X(1,LP)-X(1,I))**2 + (X(2,LP)-X(2,I))**2
A + (X(3,LP)-X(3,I))**2 )
C RJL
R(2,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) =
A DSQRT( (X(1,JP)-X(1,LP))**2 + (X(2,JP)-X(2,LP))**2
A + (X(3,JP)-X(3,LP))**2 )
IF (J_TYPE.EQ.1.AND.I_TYPE.EQ.2) THEN
CC = R(1,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE)
R(1,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) =
A R(2,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE)
R(2,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) = CC
ENDIF
C RIF: Distance to the additional force center between i-j
R(3,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) =
A DSQRT( (X(1,LP)-0.5*X(1,JP)-0.5*X(1,I))**2+
A (X(2,LP)-0.5*X(2,JP)-0.5*X(2,I))**2+
A (X(3,LP)-0.5*X(3,JP)-0.5*X(3,I))**2 )
C Wired-in formula for effective coordination number:
CC = 2.0 * DEXP( -0.0478 * (
A (R(1,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE)
A + R(2,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE))
A / R(4,1,NIT(IPHASE),IPHASE) ) ** 7.16 )
R(4,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) = CC
C ***** for C-C interaction ******
if (real_space_cutoff
C A .and.i_type.eq.j_type
A .and.i_type.eq.2.and.j_type.eq.2
A ) then
IF (R(1,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE).LT.
A RCUT(NTYPE(IOBJ,IT(1,NIT(IPHASE),IPHASE)),L_TYPE,2).AND.
A R(2,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE).LT.
A RCUT(NTYPE(IOBJ,IT(2,NIT(IPHASE),IPHASE)),L_TYPE,2)) THEN
R(5,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) = 1.0
ETAIJ = ETAIJ + CC
GOTO 50
ELSE
GOTO 55
ENDIF
endif
C ***** for C-C interaction ******
C If the increment is too small, it does not
C qualify for a screening atom -> reject it
IF (CC.LT.CUT_INC(1)) THEN
GOTO 55
ELSEIF (CC.LT.CUT_INC(2)) THEN
C Smoother for the increment rejection:
R(5,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) =
A DCOS( PI/2 * (CC-CUT_INC(2)) /
A (CUT_INC(1)-CUT_INC(2)) ) ** EXPONENT
ELSE
R(5,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE) = 1.0
ENDIF
ETAIJ = ETAIJ + R(5,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE)
A * R(4,IT(3,NIT(IPHASE),IPHASE),NIT(IPHASE),IPHASE)
C If the screening is already too strong, reject the interaction:
IF (ETAIJ.GT.CUT_ETA(2)) GOTO 45
GOTO 50
55 IT(3,NIT(IPHASE),IPHASE) = IT(3,NIT(IPHASE),IPHASE) - 1
50 CONTINUE
CC = 1.D0-DTANH(ETAIJ)
C ***** for C-C interaction ******
if (real_space_cutoff
C A .and.i_type.eq.j_type
A .and.i_type.eq.2.and.j_type.eq.2
A ) then
IF (R(4,1,NIT(IPHASE),IPHASE).GT.RCUT(I_TYPE,J_TYPE,1)) THEN
R(5,1,NIT(IPHASE),IPHASE) = DCOS ( PI/2 *
A (R(4,1,NIT(IPHASE),IPHASE)-RCUT(I_TYPE,J_TYPE,1)) /
A (RCUT(I_TYPE,J_TYPE,2)-RCUT(I_TYPE,J_TYPE,1)) ) ** 2
ELSE
R(5,1,NIT(IPHASE),IPHASE) = 1.
ENDIF
GOTO 60
endif
C ***** for C-C interaction ******
C Smoother for the interaction rejection:
IF (ETAIJ.GT.CUT_ETA(1)) THEN
C fifth field of the first record is the smoother:
R(5,1,NIT(IPHASE),IPHASE) = DCOS( PI/2 * (ETAIJ-CUT_ETA(1))
A / (CUT_ETA(2)-CUT_ETA(1)) ) ** EXPONENT
ELSE
R(5,1,NIT(IPHASE),IPHASE) = 1.
ENDIF
60 IF (I.EQ.J) THEN
G(IPHASE,I,J_TYPE) = G(IPHASE,I,J_TYPE)
A + CC * R(5,1,NIT(IPHASE),IPHASE)
ELSE
G(IPHASE,I,J_TYPE) = G(IPHASE,I,J_TYPE)
A + CC * R(5,1,NIT(IPHASE),IPHASE)
G(IPHASE,J,I_TYPE) = G(IPHASE,J,I_TYPE)
A + CC * R(5,1,NIT(IPHASE),IPHASE)
ENDIF
C ***
C IF (IPHASE.EQ.37.AND.(I.EQ.1.OR.J.EQ.1)) PRINT *,
C A ETAIJ, CC, R(5,1,NIT(IPHASE),IPHASE)
C ***
GOTO 40
45 NIT(IPHASE) = NIT(IPHASE) - 1
40 CONTINUE
C Sort into classes:
DO J = 1,I-1
IF ( ABS(G(IPHASE,I,1)-G(IPHASE,J,1)).LT.1E-5.AND.
A ABS(G(IPHASE,I,2)-G(IPHASE,J,2)).LT.1E-5.AND.
A I_TYPE.EQ.NTYPE(IOBJ,J) ) THEN
MYCLASS(IPHASE,I) = MYCLASS(IPHASE,J)
MEMCLASS(IPHASE,MYCLASS(IPHASE,I)) =
A MEMCLASS(IPHASE,MYCLASS(IPHASE,I)) + 1
GOTO 25
ENDIF
ENDDO
C A new class is found:
NCLASS(IPHASE) = NCLASS(IPHASE) + 1
MYCLASS(IPHASE,I) = NCLASS(IPHASE)
MEMCLASS(IPHASE,MYCLASS(IPHASE,I)) = 1
C Study G(r) of this class:
IF (IPHASE.EQ.NPHASE_OBJ(IOBJ,1)) CALL GR(IPHASE,I)
25 IF (ITOT.EQ.1) WRITE (*, '(1X,A2,I3,1X,4F14.10)')
A NAME_ATOM(I_TYPE), MYCLASS(IPHASE,I),
A G(IPHASE,I,1), G(IPHASE,I,2), G(IPHASE,I,1)+G(IPHASE,I,2),
A CHARGE_INI(IPHASE,I)
20 CONTINUE
WRITE (*,'(" Total of",I5," unique interactions.",/)')
A NIT(IPHASE)
IF (IPHASE.EQ.NPHASE_OBJ(IOBJ,1)) CLOSE(LP_GR)
RETURN
END
LOGICAL FUNCTION EQUIVALENT_BOND (IOBJ,I,X1,Y1,Z1,X2,Y2,Z2)
C ------------------------------------------------------
C TEST IF TWO BONDS OF THE SAME LENGTH ARE EQUIVALENT,
C I.E., THEY CAN BE RELATED BY A POINT GROUP OPERATION:
C LET ME DO A STRANGE CHECKSUM ON SURROUNDING ATOMS.
C ------------------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IOBJ,I,MAX_N1,MAX_N2,MAX_N3,II,JJ,KK,N
DOUBLE PRECISION X1,Y1,Z1,X2,Y2,Z2,ACUT,ACUT2,
A COEFFICIENT(MMT),SUM1,SUM2,XX,YY,ZZ
ACUT = 1.84542391
ACUT2 = ACUT*ACUT
MAX_N1 =
A ACUT*DSQRT(P(IOBJ,1,1)**2+P(IOBJ,2,1)**2+P(IOBJ,3,1)**2)+1
MAX_N2 =
A ACUT*DSQRT(P(IOBJ,1,2)**2+P(IOBJ,2,2)**2+P(IOBJ,3,2)**2)+1
MAX_N3 =
A ACUT*DSQRT(P(IOBJ,1,3)**2+P(IOBJ,2,3)**2+P(IOBJ,3,3)**2)+1
SUM1 = 0.
SUM2 = 0.
COEFFICIENT(1) = 0.97541
COEFFICIENT(2) = 2.29377
DO 10 II = -MAX_N1,MAX_N1
DO 10 JJ = -MAX_N2,MAX_N2
DO 10 KK = -MAX_N3,MAX_N3
DO 10 N = 1,NPA(IOBJ)
C A sphere of atoms which encloses the two bonds
XX = II*A(IOBJ,1,1) + JJ*A(IOBJ,2,1) + KK*A(IOBJ,3,1)
A + TAU(IOBJ,N,1) - TAU(IOBJ,I,1)
YY = II*A(IOBJ,1,2) + JJ*A(IOBJ,2,2) + KK*A(IOBJ,3,2)
A + TAU(IOBJ,N,2) - TAU(IOBJ,I,2)
ZZ = II*A(IOBJ,1,3) + JJ*A(IOBJ,2,3) + KK*A(IOBJ,3,3)
A + TAU(IOBJ,N,3) - TAU(IOBJ,I,3)
IF (XX*XX+YY*YY+ZZ*ZZ.LT.ACUT2) THEN
SUM1 = SUM1 + COEFFICIENT(NTYPE(IOBJ,N))
A / ( 1.3213 + XX*XX + YY*YY + ZZ*ZZ +
A COEFFICIENT(3-NTYPE(IOBJ,N)) * (XX*X1 + YY*Y1 + ZZ*Z1) )
SUM2 = SUM2 + COEFFICIENT(NTYPE(IOBJ,N))
A / ( 1.3213 + XX*XX + YY*YY + ZZ*ZZ +
A COEFFICIENT(3-NTYPE(IOBJ,N)) * (XX*X2 + YY*Y2 + ZZ*Z2) )
ENDIF
10 CONTINUE
EQUIVALENT_BOND = ABS((SUM1-SUM2)/SUM1).LT.1D-7
RETURN
END
SUBROUTINE GR(IPHASE,I)
C ------------------------------------------
C ANALYZE THE NEIGHBORS OF ATOM I UP TO MNG
C ------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER MNG,MLG
DOUBLE PRECISION TINY,RSCREEN
PARAMETER (MNG=10,MLG=MNG,TINY=1D-7,RSCREEN=2.D0)
INTEGER IPHASE,I,J,K,L,IOBJ,MAXCOUNT,II,JJ,KK,N,NT,
A MIT,J_TYPE,ICOUNT(MNG),NCOUNT(MNG)
DOUBLE PRECISION XX,YY,ZZ,CC,ETAIJ,RX(3,MNG),R2(MNG),
A GCHECK(MMT)
LOGICAL EQUIVALENT_BOND
EXTERNAL EQUIVALENT_BOND
IOBJ = MYOBJ(IPHASE)
MAXCOUNT = 0
DO 10 II = MLG,-MLG,-1
DO 10 JJ = MLG,-MLG,-1
DO 10 KK = MLG,-MLG,-1
DO 10 N = 1,NPA(IOBJ)
IF (II.EQ.0.AND.JJ.EQ.0.AND.KK.EQ.0.AND.N.EQ.I) GOTO 10
NT = NTYPE(IOBJ,N)
XX = II*A(IOBJ,1,1) + JJ*A(IOBJ,2,1) + KK*A(IOBJ,3,1)
A + TAU(IOBJ,N,1) - TAU(IOBJ,I,1)
YY = II*A(IOBJ,1,2) + JJ*A(IOBJ,2,2) + KK*A(IOBJ,3,2)
A + TAU(IOBJ,N,2) - TAU(IOBJ,I,2)
ZZ = II*A(IOBJ,1,3) + JJ*A(IOBJ,2,3) + KK*A(IOBJ,3,3)
A + TAU(IOBJ,N,3) - TAU(IOBJ,I,3)
CC = XX*XX + YY*YY + ZZ*ZZ
IF (CC.GT.RSCREEN*RSCREEN) GOTO 10
IF ( MAXCOUNT.EQ.0.OR.
A (MAXCOUNT.GE.1.AND.CC-R2(MAXCOUNT).GT.TINY) ) THEN
IF (MAXCOUNT.LT.MNG) THEN
MAXCOUNT = MAXCOUNT + 1
RX(1,MAXCOUNT) = XX
RX(2,MAXCOUNT) = YY
RX(3,MAXCOUNT) = ZZ
R2(MAXCOUNT) = CC
ICOUNT(MAXCOUNT) = 1
NCOUNT(MAXCOUNT) = NT
ENDIF
GOTO 10
ENDIF
C Something will happen: newcomers have higher priority
DO K = MAXCOUNT,0,-1
IF (K.EQ.0.OR.(K.GE.1.AND.CC-R2(K).GT.TINY)) THEN
C Move the higher ones
IF (MAXCOUNT.LT.MNG) MAXCOUNT = MAXCOUNT+1
DO J = MAXCOUNT,K+2,-1
RX(1,J) = RX(1,J-1)
RX(2,J) = RX(2,J-1)
RX(3,J) = RX(3,J-1)
R2(J) = R2(J-1)
ICOUNT(J) = ICOUNT(J-1)
NCOUNT(J) = NCOUNT(J-1)
ENDDO
RX(1,K+1) = XX
RX(2,K+1) = YY
RX(3,K+1) = ZZ
R2(K+1) = CC
ICOUNT(K+1) = 1
NCOUNT(K+1) = NT
GOTO 10
ELSEIF ( ABS(CC-R2(K)).LT.TINY.AND.
A NT.EQ.NCOUNT(K) ) THEN
C See note. 3
IF (EQUIVALENT_BOND(IOBJ,I,XX,YY,ZZ,RX(1,K),RX(2,K),RX(3,K))) THEN
ICOUNT(K) = ICOUNT(K)+1
GOTO 10
ENDIF
ENDIF
ENDDO
10 CONTINUE
WRITE (LP_GR,'(/,"Class",I3,", Atom No.",I4,", ",A2,
A ", g(Si) =",F9.6,", g(C) =",F9.6)'), MYCLASS(IPHASE,I),I,
A NAME_ATOM(NTYPE(IOBJ,I)), G(IPHASE,I,1), G(IPHASE,I,2)
WRITE (LP_GR,
A '("Neig. Reduced_R Count Type DX DY DZ",
A " ETA Coord.")')
GCHECK(1) = 0.
GCHECK(2) = 0.
DO K = 1,MAXCOUNT
C Go through lists to see if "K" interacts with
C I and contributes to I's coordination number:
DO MIT = 1,NIT(IPHASE)
IF ( (IT(1,MIT,IPHASE).EQ.I.AND.
A ABS(RX(1,K)-R(1,1,MIT,IPHASE)/AA(IPHASE)).LT.TINY.AND.
A ABS(RX(2,K)-R(2,1,MIT,IPHASE)/AA(IPHASE)).LT.TINY.AND.
A ABS(RX(3,K)-R(3,1,MIT,IPHASE)/AA(IPHASE)).LT.TINY).OR.
A (IT(2,MIT,IPHASE).EQ.I.AND.
A ABS(RX(1,K)+R(1,1,MIT,IPHASE)/AA(IPHASE)).LT.TINY.AND.
A ABS(RX(2,K)+R(2,1,MIT,IPHASE)/AA(IPHASE)).LT.TINY.AND.
A ABS(RX(3,K)+R(3,1,MIT,IPHASE)/AA(IPHASE)).LT.TINY) ) THEN
ETAIJ = 0.D0
DO L = 2,IT(3,MIT,IPHASE)
ETAIJ = ETAIJ + R(5,L,MIT,IPHASE) * R(4,L,MIT,IPHASE)
ENDDO
GOTO 20
ENDIF
ENDDO
20 IF (MIT.GT.NIT(IPHASE)) THEN
ETAIJ = -9.9D0
CC = 0.D0
ELSE
CC = R(5,1,MIT,IPHASE) * (1.D0-DTANH(ETAIJ))
ENDIF
GCHECK(NCOUNT(K)) = GCHECK(NCOUNT(K)) + CC*ICOUNT(K)
WRITE (LP_GR,'(I4,3X,F9.7,I6,4X,A2,2X,3F9.5,F9.5,F9.6)')
A K, DSQRT(R2(K)), ICOUNT(K), NAME_ATOM(NCOUNT(K)),
A RX(1,K), RX(2,K), RX(3,K), ETAIJ, CC
ENDDO
C Coordination checksum:
DO J_TYPE = 1,2
IF (ABS(GCHECK(J_TYPE)-G(IPHASE,I,J_TYPE)).GT.TINY) THEN
WRITE (*, '("gr(): coordination checksum failure for atom",
A I4," in IPHASE =",I3)') I,IPHASE
STOP
ENDIF
ENDDO
RETURN
END
DOUBLE PRECISION FUNCTION FE(GUESS)
C -----------------------------------------
C GET THE ERROR OF THE GUESS PARAMETER SET
C -----------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IP,IOP,IOBJ,NK,I,N,IDEF,ISTEP,IPHASE
DOUBLE PRECISION DELTA_CHARGE,RATIO_MIX,EGAMMA,GUESS(NOP)
C Feed GUESS parameters into the Generic Table:
IP = 1
DO IOP = 1,15,7
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 = 16,24
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)
UTC(1) = GUESS(IP+2)
UTC(2) = GUESS(IP+3)
S_SI = 1.
S_C = 1.
DO IOP = 1,NTO
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)**1.*B1(IOP,1)**(1-1.)
B2(IOP,2) = B2(IOP+7,2)**1.*B2(IOP,1)**(1-1.)
B3(IOP,2) = B3(IOP+7,2)**1.*B3(IOP,1)**(1-1.)
ENDDO
DO IOP = 9,14
B1(IOP,1) = B1(IOP-7,1)**1.*B1(IOP,2)**(1-1.)
B2(IOP,1) = B2(IOP-7,1)**1.*B2(IOP,2)**(1-1.)
B3(IOP,1) = B3(IOP-7,1)**1.*B3(IOP,2)**(1-1.)
ENDDO
DO IOP = 17,22
B1(IOP,1) = DSQRT(B1(IOP-15,1)*B1(IOP-8,1))
B2(IOP,1) = DSQRT(B2(IOP-15,1)*B2(IOP-8,1))
B3(IOP,1) = DSQRT(B3(IOP-15,1)*B3(IOP-8,1))
B1(IOP,2) = DSQRT(B1(IOP-15,2)*B1(IOP-8,2))
B2(IOP,2) = DSQRT(B2(IOP-15,2)*B2(IOP-8,2))
B3(IOP,2) = DSQRT(B3(IOP-15,2)*B3(IOP-8,2))
DELTA(IOP,1,1) = DELTA(IOP-15,1,1)*GUESS(IP+1)
A +DELTA(IOP-8,2,1)*(1-GUESS(IP+1))
DELTA(IOP,2,2) = DELTA(IOP-8,2,2)*GUESS(IP+1)
A +DELTA(IOP-15,1,2)*(1-GUESS(IP+1))
DELTA(IOP,1,2) = DELTA(IOP-15,1,2)*GUESS(IP+1)
A +DELTA(IOP-8,2,2)*(1-GUESS(IP+1))
DELTA(IOP,2,1) = DELTA(IOP-8,2,1)*GUESS(IP+1)
A +DELTA(IOP-15,1,1)*(1-GUESS(IP+1))
ENDDO
C Screening to Es,Ep(Si) = Es,Ep(C)
B1(16,1) = B1(17,1)
B2(16,1) = B2(17,1)
B3(16,1) = B3(17,1)
DELTA(16,1,1) = DELTA(17,1,1)
DELTA(16,1,2) = DELTA(17,1,2)
B1(16,2) = B1(17,2)
B2(16,2) = B2(17,2)
B3(16,2) = B3(17,2)
DELTA(16,2,1) = DELTA(17,2,1)
DELTA(16,2,2) = DELTA(17,2,2)
C Screening to Qs,Qp(*)
DO IOP = 1,15,7
C B1(IOP,1) = B1(IOP+1,1)
C B2(IOP,1) = B2(IOP+1,1)
C B3(IOP,1) = B3(IOP+1,1)
C DELTA(IOP,1,1) = DELTA(IOP+1,1,1)
C DELTA(IOP,1,2) = DELTA(IOP+1,1,2)
C B1(IOP,2) = B1(IOP+1,2)
C B2(IOP,2) = B2(IOP+1,2)
C B3(IOP,2) = B3(IOP+1,2)
C DELTA(IOP,2,1) = DELTA(IOP+1,2,1)
C DELTA(IOP,2,2) = DELTA(IOP+1,2,2)
B1(IOP,1) = 2.0
B2(IOP,1) = 0.0478
B3(IOP,1) = 7.16
DELTA(IOP,1,1) = 0.
DELTA(IOP,1,2) = 0.
B1(IOP,2) = 2.0
B2(IOP,2) = 0.0478
B3(IOP,2) = 7.16
DELTA(IOP,2,1) = 0.
DELTA(IOP,2,2) = 0.
ENDDO
C Screening to Phi(Si->C) = Phi(C->Si)
B1(23,1) = B1(22,1)
B2(23,1) = B2(22,1)
B3(23,1) = B3(22,1)
DELTA(23,1,1) = DELTA(22,1,1)
DELTA(23,1,2) = DELTA(22,1,2)
B1(23,2) = B1(22,2)
B2(23,2) = B2(22,2)
B3(23,2) = B3(22,2)
DELTA(23,2,1) = DELTA(22,2,1)
DELTA(23,2,2) = DELTA(22,2,2)
C Screening to ps_sigma = sp_sigma
B1(24,1) = B1(19,1)
B2(24,1) = B2(19,1)
B3(24,1) = B3(19,1)
DELTA(24,1,1) = DELTA(19,1,1)
DELTA(24,1,2) = DELTA(19,1,2)
B1(24,2) = B1(19,2)
B2(24,2) = B2(19,2)
B3(24,2) = B3(19,2)
DELTA(24,2,1) = DELTA(19,2,1)
DELTA(24,2,2) = DELTA(19,2,2)
IF (ITOT.EQ.MAXITER) WRAPPING_UP = .TRUE.
DO_BAND = (.NOT.FIX_EBS).OR.ITOT.EQ.1.OR.BAND_ERROR_RATIO.GT.0
COH_ERROR = 0.D0
IF (DO_BAND) BAND_ERROR = 0.D0
DO IPHASE = 1,NPHASE
IOBJ = MYOBJ(IPHASE)
C Initialize the charges for this run:
IF (ITOT.EQ.1.OR.(.NOT.(USE_LAST_STEP_CHARGES.OR.FIX_EBS))) THEN
DO N = 1,NPA(IOBJ)
CHARGE(IPHASE,N) = CHARGE_INI(IPHASE,N)
ENDDO
ENDIF
DO N = 1,NPA(IOBJ)
CHARGE_NEW(IPHASE,N) = CHARGE(IPHASE,N)
ENDDO
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.OR.ITOT.EQ.1 ) THEN
ISTEP = 1
C Construct the overlap matrix and calculate Erep and Eu:
10 CALL OVERLAP (IPHASE,ISTEP)
IF ((.NOT.FIX_EBS).OR.ITOT.EQ.1) CALL CALC_EBS(IPHASE)
ETB(IPHASE) = EREP(IPHASE) + EBS(IPHASE) + EU(IPHASE)
IF (UTC(1).EQ.0.AND.UTC(2).EQ.0) GOTO 20
DELTA_CHARGE = 0.
C print *, 'old ', (CHARGE(IPHASE,N),N=1,NPA(IOBJ))
C print *, 'new ', (CHARGE_NEW(IPHASE,N),N=1,NPA(IOBJ))
DO N = 1,NPA(IOBJ)
DELTA_CHARGE = DELTA_CHARGE +
A ABS(CHARGE_NEW(IPHASE,N)-CHARGE(IPHASE,N))
ENDDO
DELTA_CHARGE = DELTA_CHARGE / NPA(IOBJ)
IF (DELTA_CHARGE.GT.CHARGE_TOLERANCE) THEN
C If the charge transfer is too large, exit the iteration:
DO N = 1,NPA(IOBJ)
IF (ABS(CHARGE(IPHASE,N)-4.).GT.3.) THEN
PRINT *, 'Warning: charge transfer too large, exit loop.'
GOTO 20
ENDIF
ENDDO
C Re-mix the charges and do a new iteration:
IF (ISTEP.LE.10) THEN
RATIO_MIX = 0.20
ELSEIF (ISTEP.LE.25) THEN
RATIO_MIX = 0.10
ELSEIF (ISTEP.LE.50) THEN
RATIO_MIX = 0.05
ELSE
GOTO 20
ENDIF
DO N = 1,NPA(IOBJ)
CHARGE(IPHASE,N) = RATIO_MIX * CHARGE_NEW(IPHASE,N) +
A (1.-RATIO_MIX) * CHARGE(IPHASE,N)
ENDDO
ISTEP = ISTEP + 1
GOTO 10
ENDIF
20 CONTINUE
C WRITE (*, '(A8," No.", I2, " converged in",I3," steps.")')
C A NAME(IOBJ)(1:INDEX(NAME(IOBJ),' ')-1),
C A IPHASE-NPHASE_OBJ(IOBJ,1)+1, ISTEP
ENDIF
IF (DO_BAND.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
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,16
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,ISTEP)
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 AND EU WILL BE EVALUATED.
C ----------------------------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IPHASE,ISTEP,IOP,IOBJ,N,I,J,K,L,I_TYPE,J_TYPE,
A L_TYPE,MIN,MAX,MIT,ID
DOUBLE PRECISION CC,CX,CY,CZ,H(NTO),ETA(NTO),FOFC,FOFSI,
A X0,X1,X2,X3,X4,Y0,Y1,Y2,Y3,Y4,XX,ETA_H,ETA_L,RR,
A RIJ,RIL,RJL,RIF,XJ,YJ,ZJ,BB
IOBJ = MYOBJ(IPHASE)
C Load in charges for Coulomb potential summation:
DO N = 1,NPA(IOBJ)
TMP_CHARGE(N) = CHARGE(IPHASE,N)
ENDDO
C On first entry, load in the bare
C dimensionless Coulomb coefficients:
IF (ISTEP.EQ.1) THEN
DO N = 1,NPA(IOBJ)*NPA(IOBJ)
POTE_MATRIX(N) = OBJ_POTE_MATRIX(IOBJ,N) / AA(IPHASE)
ENDDO
ENDIF
C IF (IPHASE.EQ.1) THEN
C print *, 'mat = ', (POTE_MATRIX(N),n=1,NPA(IOBJ)*NPA(IOBJ))
C print *, 'q = ', (TMP_CHARGE(N),n=1,NPA(IOBJ))
C print *, TOTAL_EWALD_ENERGY(1.D0,NPA(IOBJ),TMP_CHARGE,POTE_MATRIX)
C ENDIF
C If not the first entry, we only need
C to modify the on-site matrix elements:
IF (ISTEP.GT.1) THEN
EU(IPHASE) =
A -TOTAL_EWALD_ENERGY(1.D0,NPA(IOBJ),TMP_CHARGE,POTE_MATRIX)
DO I = 1,NPA(IOBJ)
I_TYPE = NTYPE(IOBJ,I)
ID = IDX(I,I)
CC = EWALD_POTENTIAL(I,1.D0,NPA(IOBJ),TMP_CHARGE,POTE_MATRIX) +
A UTC(I_TYPE)*(CHARGE(IPHASE,I)-4.) +
A SUPPORT(I_TYPE)*(CHARGE(IPHASE,I)-4.)**(SUPPORT_EXPONENT-1)
TM(I,I,ID,1,1) = ES(I) + CC
TM(I,I,ID,2,2) = EP(I) + CC
TM(I,I,ID,3,3) = EP(I) + CC
TM(I,I,ID,4,4) = EP(I) + CC
EU(IPHASE) = EU(IPHASE) +
A UTC(I_TYPE)*(CHARGE(IPHASE,I)-4.)**2/2. -
A UTC(I_TYPE)*(CHARGE(IPHASE,I)-4.)*CHARGE(IPHASE,I) +
A SUPPORT(I_TYPE)*(CHARGE(IPHASE,I)-4.)**SUPPORT_EXPONENT
A /SUPPORT_EXPONENT -
A SUPPORT(I_TYPE)*(CHARGE(IPHASE,I)-4.)**(SUPPORT_EXPONENT-1)
A *CHARGE(IPHASE,I)
ENDDO
EU(IPHASE) = EU(IPHASE) / NPA(IOBJ)
RETURN
ENDIF
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 non-selfconsistent the onsite energies.
DO 10 MIT = 1,NIT(IPHASE)
I = IT(1,MIT,IPHASE)
I_TYPE = NTYPE(IOBJ,I)
J = IT(2,MIT,IPHASE)
J_TYPE = NTYPE(IOBJ,J)
XJ = R(1,1,MIT,IPHASE)
YJ = R(2,1,MIT,IPHASE)
ZJ = R(3,1,MIT,IPHASE)
RIJ = R(4,1,MIT,IPHASE)
C Command Center:
IF (I_TYPE.EQ.1.AND.J_TYPE.EQ.1) THEN
MIN = 1
IF (.NOT.DO_BAND) MIN = 7
MAX = 7
ELSEIF (I_TYPE.EQ.2.AND.J_TYPE.EQ.2) THEN
MIN = 8
IF (.NOT.DO_BAND) MIN = 14
MAX = 14
ELSE
MIN = 15
MAX = 24
IF (.NOT.DO_BAND) MIN = 22
IF (.NOT.DO_BAND) MAX = 23
ENDIF
C First, calculate the interaction strength:
DO IOP = MIN,MAX
C Rescale the distances
T0(1) = 4.0000000
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))
ETA(IOP) = 0.
ENDDO
C Then, calculate the screening by looping over L's
DO L = 2,IT(3,MIT,IPHASE)
L_TYPE = NTYPE(IOBJ,IT(2+L,MIT,IPHASE))
RIL = R(1,L,MIT,IPHASE)
RJL = R(2,L,MIT,IPHASE)
RIF = R(3,L,MIT,IPHASE)
DO IOP = MIN,MAX
CC = ( D(IOP,L_TYPE) * RIL +
A 1./D(IOP,L_TYPE) * RJL +
A F(IOP,L_TYPE) * RIF ) / RIJ
ETA(IOP) = ETA(IOP) + R(5,L,MIT,IPHASE) *
A B1(IOP,L_TYPE) * DEXP(-B2(IOP,L_TYPE)*CC**B3(IOP,L_TYPE))
ENDDO
ENDDO
DO IOP = MIN,MAX
IF (IOP.EQ.1.OR.IOP.EQ.8.OR.IOP.EQ.15) THEN
C Uij coeff. softening: only for first neighbor
ETA_L = 0.8
ETA_H = 1.2
IF (ETA(IOP).GT.ETA_H) THEN
C No softening, let it be bare:
H(IOP) = 0.D0
ELSEIF (ETA(IOP).LT.ETA_L) THEN
C First neighbor, soften Uij to what I want:
H(IOP) = H(IOP) - ENERGY_TO_EV / R(4,1,MIT,IPHASE)
ELSE
C smoothly link to 1/r behavior:
H(IOP) = (H(IOP) - ENERGY_TO_EV / R(4,1,MIT,IPHASE)) *
A DCOS(PI/2 *(ETA(IOP)-ETA_L)/(ETA_H-ETA_L))**EXPONENT
ENDIF
ELSE
C Screen the interaction:
XX = DEXP(ETA(IOP))
H(IOP) = H(IOP)*2./(XX+1./XX)/XX * R(5,1,MIT,IPHASE)
ENDIF
C To avoid numerical problems:
IF (ABS(H(IOP)).GT.500.) H(IOP) = 0.
ENDDO
IF (DO_BAND) 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
CX = XJ/RIJ
CY = YJ/RIJ
CZ = ZJ/RIJ
IF (I_TYPE.EQ.J_TYPE) THEN
IF (I.EQ.J) THEN
IF (DO_BAND) THEN
C POTE_MATRIX((I-1)*NPA(IOBJ)+J) =
C A POTE_MATRIX((I-1)*NPA(IOBJ)+J) + H(MIN)
ES(I) = ES(I) + H(MIN+1)
EP(I) = EP(I) + H(MIN+1)
ENDIF
PHI(I,J_TYPE) = PHI(I,J_TYPE) + H(MAX)
ELSE
IF (DO_BAND) THEN
C POTE_MATRIX((I-1)*NPA(IOBJ)+J) =
C A POTE_MATRIX((I-1)*NPA(IOBJ)+J) + H(MIN)
C POTE_MATRIX((J-1)*NPA(IOBJ)+I) =
C A POTE_MATRIX((J-1)*NPA(IOBJ)+I) + H(MIN)
ES(I) = ES(I) + H(MIN+1)
EP(I) = EP(I) + H(MIN+1)
ES(J) = ES(J) + H(MIN+1)
EP(J) = EP(J) + H(MIN+1)
ENDIF
PHI(I,J_TYPE) = PHI(I,J_TYPE) + H(MAX)
PHI(J,I_TYPE) = PHI(J,I_TYPE) + H(MAX)
ENDIF
IF (DO_BAND) THEN
TM(I,J,ID,1,1) = H(MIN+1+1)
TM(I,J,ID,1,2) = H(MIN+1+2)*CX
TM(I,J,ID,1,3) = H(MIN+1+2)*CY
TM(I,J,ID,1,4) = H(MIN+1+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+1+3)+(1-CX*CX)*H(MIN+1+4)
C px-py
TM(I,J,ID,2,3) = CX*CY*(H(MIN+1+3)-H(MIN+1+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+1+3)-H(MIN+1+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+1+3)+(1-CY*CY)*H(MIN+1+4)
C py-pz
TM(I,J,ID,3,4) = CY*CZ*(H(MIN+1+3)-H(MIN+1+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+1+3)+(1-CZ*CZ)*H(MIN+1+4)
ENDIF
ELSE
IF (DO_BAND) THEN
C POTE_MATRIX((I-1)*NPA(IOBJ)+J) =
C A POTE_MATRIX((I-1)*NPA(IOBJ)+J) + H(MIN)
C POTE_MATRIX((J-1)*NPA(IOBJ)+I) =
C A POTE_MATRIX((J-1)*NPA(IOBJ)+I) + H(MIN)
ES(I) = ES(I) + H(MIN+1)
EP(I) = EP(I) + H(MIN+1)
ES(J) = ES(J) + H(MIN+2)
EP(J) = EP(J) + H(MIN+2)
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+1+2)
TM(I,J,ID,1,2) = H(MIN+1+3) * CX
TM(I,J,ID,1,3) = H(MIN+1+3) * CY
TM(I,J,ID,1,4) = H(MIN+1+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+1+4)+(1-CX*CX)*H(MIN+1+5)
C px-py
TM(I,J,ID,2,3) = CX*CY*(H(MIN+1+4)-H(MIN+1+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+1+4)-H(MIN+1+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+1+4)+(1-CY*CY)*H(MIN+1+5)
C py-pz
TM(I,J,ID,3,4) = CY*CZ*(H(MIN+1+4)-H(MIN+1+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+1+4)+(1-CZ*CZ)*H(MIN+1+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
10 CONTINUE
IF (DO_BAND) 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
I_TYPE = NTYPE(IOBJ,I)
CC = EWALD_POTENTIAL(I,1.D0,NPA(IOBJ),TMP_CHARGE,POTE_MATRIX) +
A UTC(I_TYPE)*(CHARGE(IPHASE,I)-4.) +
A SUPPORT(I_TYPE)*(CHARGE(IPHASE,I)-4.)**(SUPPORT_EXPONENT-1)
TM(I,I,ID,1,1) = ES(I) + CC
TM(I,I,ID,2,2) = EP(I) + CC
TM(I,I,ID,3,3) = EP(I) + CC
TM(I,I,ID,4,4) = EP(I) + CC
DX(I,I,ID,1) = 0.
DX(I,I,ID,2) = 0.
DX(I,I,ID,3) = 0.
ENDDO
ENDIF
BB = TOTAL_EWALD_ENERGY(1.D0,NPA(IOBJ),TMP_CHARGE,POTE_MATRIX)
EU(IPHASE) = -BB
print *, iphase, BB, 'orig deficit = ', EU(IPHASE)
DO I = 1,NPA(IOBJ)
I_TYPE = NTYPE(IOBJ,I)
BB = BB - EWALD_POTENTIAL(I,1.D0,NPA(IOBJ),TMP_CHARGE,POTE_MATRIX)
A * charge(IPHASE,I)
EU(IPHASE) = EU(IPHASE) +
A UTC(I_TYPE)*(CHARGE(IPHASE,I)-4.)**2/2. -
A UTC(I_TYPE)*(CHARGE(IPHASE,I)-4.)*CHARGE(IPHASE,I) +
A SUPPORT(I_TYPE)*(CHARGE(IPHASE,I)-4.)**SUPPORT_EXPONENT
A /SUPPORT_EXPONENT -
A SUPPORT(I_TYPE)*(CHARGE(IPHASE,I)-4.)**(SUPPORT_EXPONENT-1)
A *CHARGE(IPHASE,I)
ENDDO
EU(IPHASE) = EU(IPHASE)/NPA(IOBJ)
IF (IPHASE.EQ.1) THEN
PRINT *, (POTE_MATRIX(N),N=1,NPA(IOBJ)*NPA(IOBJ))
PRINT *, (tmp_charge(N),N=1,NPA(IOBJ))
PRINT *, (charge(iphase,N),N=1,NPA(IOBJ))
C STOP
ENDIF
print *, 'real deficit = ', BB
stop
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)
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
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 ----------------------------------------------------------------
INTEGER IPHASE,NA,NOB,II,JJ,I,J,ID,IP,JP,INFO,LENAI
DOUBLE PRECISION FKX,FKY,FKZ,PHASE,AUZ(6*MMD)
COMPLEX*16 TMAT(MMD,MMD),FACTOR,AP(MMD*(MMD+1)/2),ADD,
A ZAUZ(2*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 --------------------------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
DOUBLE PRECISION CUTOFF
PARAMETER (CUTOFF=30.D0)
INTEGER IPHASE,IOBJ,I,NORB
DOUBLE PRECISION EF,EF1,FVALUE,DERIVATIVE,DE,DERF
EXTERNAL DERF
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 ----------------------------------------------------------
C ********************************
C GLOBAL VARIABLES OF FitAll V3.0
C ********************************
IMPLICIT NONE
INTEGER MMO,MMS,MMT,MMF,MME,MMA,MMI,MSC,MAC,MMC,MZK,
A MBK,MMB,MSG,MMD,NTO,NOP
PARAMETER (MMO=32,MMS=128,MMT=2,MMF=30,MME=64,MMA=9,
A MMI=128,MSC=64,MAC=128,MMC=MMA,MZK=1024,MBK=8,
A MMB=4*MMA,MSG=8,MMD=4*MMA,NTO=24,NOP=64)
DOUBLE PRECISION G0,EXPONENT,EPS,PI,VACUUM_PM,
A ELECTRON_CHARGE,ENERGY_TO_EV
PARAMETER (G0=4.4165281974,EXPONENT=1.5,EPS=1D-6,
A PI=3.14159265358979323846D0,VACUUM_PM=
A 8.854187817D-12,ELECTRON_CHARGE=1.602176454D-19,
A ENERGY_TO_EV=ELECTRON_CHARGE*1D10/4/PI/VACUUM_PM)
C ------------------------------------------------------------
C MMO: MAXIMUM NUMBER OF OBJECT GROUPS
C MMS: MAXIMUM NUMBER OF REALIZED STRUCTURES
C MMT: MAXIMUM NUMBER OF ATOMIC TYPES
C MMF: MAXIMUM NUMBER OF DEFORMATION MODES
C MME: MAXIMUM NUMBER OF ELEMENTS IN A DEFORMATION MODE
C MMA: MAXIMUM NUMBER OF ATOMS PER UNIT CELL
C MMI: MAXIMUM NUMBER OF UNIQUE INTERACTIONS IN THE SYSTEM
C MSC: MAXIMUM NUMBER OF SCREENING ATOMS PER INTERACTION
C MAC: MAXIMUM NUMBER OF PAIR IMAGE INTERACTIONS
C MMC: MAXIMUM NUMBER OF EQUIVALENT ATOMIC CLASSES
C MZK: MAXIMUM NUMBER OF INTEGRATION K-POINTS FOR OBJECTS
C MBK: MAXIMUN NUMBER OF TARGET K-POINTS FOR BAND STRUCTURES
C MMB: MAXIMUN NUMBER OF BANDS TO FIT FOR EACH TARGET K-POINT
C MSG: MAXIMUM NUMBER OF BAND SEGMENTS TO BE PLOTTED OUT
C NTO: NUMBER OF ROWS ON THE GENERIC PARAMETER TABLE
C NOP: NUMBER OF PARAMETERS FOR OPTIMIZATION
C G0: COORDINATION NUMBER OF REFERENCE STRUCTURE (DIAMOND)
C ------------------------------------------------------------
COMMON /INTEGER/ ITOT,MAXITER,IOFTEN,NPA(MMO),MYOBJ(MMS),
A NTYPE(MMO,MMA),IT(3+MSC,MMI,MMS),NPHASE_OBJ(MMO,2),
A NPHASE,NOBJ,NIT(MMS),NATOM(MMO,MMT),NCLASS(MMS),
A MYCLASS(MMS,MMA),MEMCLASS(MMS,MMC),IDX(MMA,MMA),
A NBKP(MMS),LINK(MMS),IA(MMS,MSG),IB(MMS,MSG),
A NS(MMS,MSG),NBA(MMS,MBK),NTH(MMS,MBK,MMB),
A NZKP(MMO),NDEF,NELE(MMF),IELE(MMF,MME)
COMMON /DOUBLE/ CHARGE_TOLERANCE,RCUT(MMT,MMT,2),ES0(MMT),
A EP0(MMT),A1(NTO),A2(NTO),A3(NTO),A4(NTO),B1(NTO,MMT),
A B2(NTO,MMT),B3(NTO,MMT),DELTA(NTO,MMT,MMT),D(NTO,MMT),
A F(NTO,MMT),CUT_ETA(2),CUT_INC(2),BOUND(NOP,2),T0(MMT),
A C0,C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC(MMT),SUPPORT(MMT),SUPPORT_EXPONENT,
A A(MMO,3,3),P(MMO,3,3),DET(MMO),AA(MMS),TAU(MMO,MMA,3),
A G(MMS,MMA,MMT),R(5,1+MSC,MMI,MMS),QCLASS(MMS,MMC),
A CHARGE_INI(MMS,MMA),ES(MMA),EP(MMA),PHI(MMA,MMT),
A TM(MMA,MMA,MAC,4,4),DX(MMA,MMA,MAC,3),W(MMD),WB(MMS),
A BKP(MMS,MBK,3),GAMMA_LDA(MMS),FERMI_LDA(MMS),
A BWEIGHT(MMS,MBK,MMB),VOLUME_LDA(MMS),BAND_ERROR,
A BTARGET(MMS,MBK,MMB),ZKP(MMO,MZK,4),ELIST(MMD*MZK),
A EDEL,EFERMI(MMS),EWEIGHT(MMD*MZK),CHARGE(MMS,MMA),
A CHARGE_NEW(MMS,MMA),EBS(MMS),EREP(MMS),EU(MMS),
A ETB(MMS),CPS(MMS),GROUND_WEIGHT(MMO),COH_ERROR,
A CURV_WEIGHT(MMO),DEF_WEIGHT(MMF),BAND_ERROR_RATIO,
A S1(MMA),S2(MMA),S3(MMA),POTE_MATRIX(MMA*MMA),
A TMP_CHARGE(MMA),OBJ_POTE_MATRIX(MMO,MMA*MMA)
COMMON /COMPLEX/ Z(MMD,MMD),ZLIST(MMD,MMD*MZK)
COMMON /LOGICAL/ WRAPPING_UP,FIX_EBS,DO_BAND,
A USE_LAST_STEP_CHARGES,REAL_SPACE_CUTOFF
COMMON /CHARACT/ BUF,NAME_PARA(NOP),NAME_ATOM(MMT),NAME(MMO),
A NAME_BAND(MMS),NAME_DEF(MMF)
DATA LP,LP_GR,LP_PARA /19,20,21/
EXTERNAL SIMPLE_EWALD,EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
INTEGER ITOT,MAXITER,IOFTEN,NPA,MYOBJ,NTYPE,IT,NPHASE_OBJ,
A NPHASE,NOBJ,NIT,NATOM,NCLASS,MYCLASS,MEMCLASS,IDX,
A NBKP,LINK,IA,IB,NS,NBA,NTH,NZKP,NDEF,NELE,IELE,
A LP,LP_GR,LP_PARA
DOUBLE PRECISION CHARGE_TOLERANCE,RCUT,ES0,EP0,A1,A2,A3,
A A4,DELTA,B1,B2,B3,D,F,CUT_ETA,CUT_INC,BOUND,T0,C0,
A C1,C2,C3,C4,D0,D1,D2,D3,D4,RIGID_SHIFT,S_SI,S_C,
A C_ES,C_EP,UTC,SUPPORT,SUPPORT_EXPONENT,A,P,DET,AA,
A TAU,G,R,QCLASS,CHARGE_INI,ES,EP,PHI,TM,DX,W,WB,BKP,
A GAMMA_LDA,FERMI_LDA,BWEIGHT,VOLUME_LDA,BTARGET,
A BAND_ERROR,ZKP,ELIST,EDEL,EFERMI,EWEIGHT,CHARGE,
A CHARGE_NEW,EBS,EREP,EU,ETB,CPS,GROUND_WEIGHT,
A CURV_WEIGHT,DEF_WEIGHT,BAND_ERROR_RATIO,COH_ERROR,
A S1,S2,S3,POTE_MATRIX,TMP_CHARGE,OBJ_POTE_MATRIX,
A EWALD_POTENTIAL,TOTAL_EWALD_ENERGY
COMPLEX*16 Z,ZLIST
LOGICAL WRAPPING_UP,FIX_EBS,DO_BAND,USE_LAST_STEP_CHARGES,
A REAL_SPACE_CUTOFF
CHARACTER *70 BUF,NAME_PARA,NAME_ATOM,NAME,NAME_DEF,NAME_BAND
INTEGER IPHASE,IOBJ,IE,NK,I,J,IP,IC,ITER,IH,IL,NORB
DOUBLE PRECISION FKX,FKY,FKZ,TOTAL,DE,EF1,EF,EH,FH,EL,FL,
A DERF,EE,DERIVATIVE,FVALUE,EORDER(MMD*MZK),SUM_CHARGE(MMT)
EXTERNAL DERF
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)
C To avoid numerical problems:
IF (W(I).GT.1/EPS.OR.W(I).LT.-1/EPS) ELIST(IE) = 0.
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.25) GOTO 20
PRINT *,
A 'calc_ebs(): Newtons method out of maximum iteration.'
C 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 atomic sites by
C eigenstates and sum up in equivalent atom classes.
DO IP = 1,NPA(IOBJ)
CHARGE_NEW (IPHASE,IP) = 0.
ENDDO
DO IC = 1,NCLASS(IPHASE)
QCLASS (IPHASE,IC) = 0.
ENDDO
TOTAL = 0.
DO I = 1,4*NPA(IOBJ)
IP = (I-1)/4+1
DO J = 1,NORB
CHARGE_NEW (IPHASE,IP) = CHARGE_NEW (IPHASE,IP) +
A (REAL(ZLIST(I,J))**2 + IMAG(ZLIST(I,J))**2) * EWEIGHT(J)
ENDDO
ENDDO
DO IP = 1,NPA(IOBJ)
CHARGE_NEW (IPHASE,IP) = CHARGE_NEW (IPHASE,IP) * 8 * NPA(IOBJ)
QCLASS(IPHASE,MYCLASS(IPHASE,IP)) =
A QCLASS(IPHASE,MYCLASS(IPHASE,IP)) + CHARGE_NEW(IPHASE,IP)
TOTAL = TOTAL + CHARGE_NEW(IPHASE,IP)
ENDDO
C Symmetrize charges within each equivalent class after IBZ
C integration (see note 2):
DO IP = 1,NPA(IOBJ)
CHARGE_NEW(IPHASE,IP) = QCLASS(IPHASE,MYCLASS(IPHASE,IP)) /
A MEMCLASS(IPHASE,MYCLASS(IPHASE,IP)) - (TOTAL/NPA(IOBJ)-4.)
ENDDO
RETURN
END
C ========================================================================
C Notes:
C
C (1) Composite cohesive energy targets should be store in
C Coh_Database/3C.cps etc. in the format (Length scale in A)
C + (Global identifier such as volume/atom) + (Cohesive
C energy/atom(eV)). Here I use the data which is the smooth
C universal binding curve fitted from raw LDA data. CZ says
C that since the raw LDA data points come with some noise,
C it is better to use the smooth curve as database. I agree
C with him because I remembered in fitting the curve a small
C change of parameters such as the domain width can induce
C big error in the bulk modulus. As a convention, the exact
C energy minimum of this phase should be put in the first record.
C
C 2. Although IBZ is enough for total energy summation, it can
C lead to fictitious results in LDOS calculation if the symmetry
C operations contains non-primitive translations, such as for the
C two Si atoms in 2H SiC, though physically equivalent, may have
C different IBZ summed charges. The reason is simply although
C O\psi(k)=\psi(Ok) has the same energy as \psi(k), the charge
C on atom i in \psi(k) could shift to an equivalent but
C different site after O\psi(k): so although \sum_{BZ}E(k) =
C N\sum_{IBZ}E(k), and \sum_{BZ}\rho(k,r) has the full symmetry
C of the crystal, \sum_{IBZ}\rho(k,r) may not. This problem has a
C simple solution: one can prove that charges belonging to
C one class of atoms (two Si above) never goes other classes for
C any transformation, so one just needs to symmetrize within each
C class to the IBZ calculated charges. In order to determine
C equivalent atomic classes without the trouble of doing space
C group operations (./Kpts/Bin/group.f), let me adopt the "unsafe"
C trick of comparing coordination numbers. An agreement to the
C sixth digit in both Si and C coordination number is a sound
C guarantee (against accident) that the two atoms are actually
C equivalent, and should have the same charge.
C
C 3. We classify an atom's neighbors first by their species
C and then distances to i. But some atoms that have the same
C bond length to i are not equivalent: like 2H-SiC (Si ABAB
C stacking + C shift), a Si has 12 "second-nearest neighbors",
C but only 6 are on the basal plane, the other 6 give a different
C coordination number contribution. Even those 6 are inequivalent
C as upper 3 + lower 3, even when the bonds themselves are
C equally screened for there is i+j symmetry.
C
C The definition of an equivalent neighbor is that if you look
C in the 2nd direction, after an axial rotation it would look
C exactly the same as in the 1st direction, which is to say that
C the two neighbors can be related by a point group operation
C centered on i. A full realization of that procedure can be
C implemented but is time consuming. My checksum test of conserved
C quantities is faster and should work in most situations.
C ========================================================================