# make a general spread binary (any or all parameters)
# see two seminal EAM-X papers by M. Daw & M. Chandross, Acta Mat (2023)
# MDaw & MChandross 13 Dec 2022

def MakeSpreadBinary():
    
    from Make_LAMMPS_eamalloy_v2 import Make_LAMMPS_eamalloy_v2
    import math
    import numpy as np

    # collect base values (usually Averagium)
    baseparams = np.loadtxt('baseparams.set')
    baser1nn = baseparams[0]
    baseEc = baseparams[1]
    baseB = baseparams[2]
    basephi0 = baseparams[3]
    basebeta = baseparams[4]
    basercut = baseparams[5]
    
    basemass = 64
    
    ntypes = 2
    
    # set both types to base values
    ones=np.ones(ntypes)
    r1nns=ones*baser1nn
    Ecs=ones*baseEc
    Bs=ones*baseB
    betas=ones*basebeta
    phi0s=ones*basephi0
    masses=ones*basemass
    names=np.array([str(itype) for itype in range(ntypes)])
    names = ["A","B"]
    
    rho0s=ones*1.
    
    # cross values
    chis=np.ones((ntypes,ntypes))
    
    # make "spread" in parameters [r1, Ec, B, beta, phi0, rho0, chi]
    dr1nn, dEc, dB, dbeta, dphi0, drho0, chi = np.loadtxt('deltas.set')
#    print(" deltas = ",dr1nn,dEc,dB,dbeta,dphi0,drho0,chi)

#   here's a summary of how the "deltas" work (see papers for complete details)
#   a "spread binary" is made in reference to a reference metal (see the
#   portion above that reads in the file "baseparams.set")
#   the is then a delta defined for each parameter (see list above read
#   in from "deltas.set")
#   for each parameter then the value for elements A and B are modified
#   from the reference metal
#   for example, the parameter r1nn sets the nearest neighbor distance
#   dr1nn is the "delta" for that parameter
#   the delta is the fractional change
#   so the nearest neighbor distances for elements A and B are
#   r1nnA = (1+dr1nn)*r1nn
#   r1nnB = (1-dr1nn)*r1nn
#   similarly for the other parameters Ec, B, beta, phi0, and rho0
#   the parameters chi are not strictly "delta"-style because they
#   involves pairs of types, but they are included in the "deltas.set"
#   just for simplicity

    factorsL = np.array([1+dr1nn,1-dr1nn])  # length scales
    factorsE = np.array([1+dEc,1-dEc]) # energy scales

#   the Length and Energy scales are set by r1nn and Ec
#   the other parameters B, beta, and phi0 are "scaled" then
#   according to the prescription below (see papers on scaling)
#   if scaled in this way, the resulting SCALED properties will be
#   the same
#   for instance, the ratio of vacancy formation energy (VFE) to
#   cohesive energy (Ec) will remain the same
    
    r1nns = r1nns*factorsL
    Ecs = Ecs*factorsE
    Bs = Bs*factorsE/factorsL**3
    betas = betas/factorsL
    phi0s = phi0s*factorsE
    
    Bs = Bs*np.array([1+dB ,1-dB])
    betas = betas*np.array([1+dbeta,1-dbeta])
    phi0s = phi0s*np.array([1+dphi0,1-dphi0])
    rho0s = rho0s*np.array([1+drho0,1-drho0])
    chis[0,1] = 1+chi
    chis[1,0] = 1+chi

#   setting the cutoff distance to 1.95*r1nn (r1nn=nearest neighbor distance)
#   should give a maximum stability for FCC relative to HCP in each
#   pure metal
#   another choice might be preferred to make HCP stable relative to FCC
#   (see first paper for discussion)
    
    rcuts = 1.95*r1nns
    
    header1 = "AB binary from Av, deltas: dr1,dEc,dB,dbeta,dphi0,drho0,chi "
    header2 = " "+str(dr1nn)+" "+str(dEc)+" "+str(dB)+" "+str(dbeta)+" "+str(dphi0)+" "+str(drho0)+" "+str(chi)
    headers = (header1,header2)
    
    Make_LAMMPS_eamalloy_v2(headers,ntypes,betas,r1nns,rho0s,phi0s,Ecs,Bs,masses,names,rcuts,chis)
    
#    print(" done ")

if __name__ == '__main__':
    import sys
    results = MakeSpreadBinary()