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
#     include "fitall.fh"
      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     ---------------------------------------------------------
#     include "fitall.fh"      
      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:
#     include "fitall.fh"
      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     ------------------------------------------------------
#     include "fitall.fh"
      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     ------------------------------------------------------
#     include "fitall.fh"
      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     ------------------------------------------
#     include "fitall.fh"
      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     -----------------------------------------
#     include "fitall.fh"
      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     ----------------------------------------------------------------
#     include "fitall.fh"
      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)
#     include "fitall.fh" 
C     ---------------------------------------------------------------
C     SUBROUTINE TO DIAGONALIZE THE HAMILTONIAN FOR A SPECIFIC K
C     THE INPUT ARE NA(NUMBER OF ATOMS), AND FKX,FKY,FKZ. THE OVERLAP 
C     MATRIX TM(MMA,MMA,MAC,4,4) IN STACKED FORM MUST BE READY. 
C     THE OUTPUT EIGENVALUES ARE WRITTEN IN W(MMD), WHICH IS ORDERD.
C     ----------------------------------------------------------------
      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     --------------------------------------------------------------
#     include "fitall.fh"
      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     ----------------------------------------------------------
#     include "fitall.fh"
      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     ========================================================================