Python EquationOfState.fit方法代码示例

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
 def fit(filename):
     configs = read(filename + '@:')
     volumes = [a.get_volume() for a in configs]
     energies = [a.get_potential_energy() for a in configs]
     eos = EquationOfState(volumes, energies)
     v0, e0, B = eos.fit()
     return (4 * v0)**(1 / 3.0)

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
    def bulk_summary(self, plot, a0):
        natoms = len(self.atoms)
        eos = EquationOfState(self.volumes, self.energies)
        v, e, B = eos.fit()
        x = (v / self.atoms.get_volume())**(1.0 / 3)

        self.log('Fit using %d points:' % len(self.energies))
        self.log('Volume per atom: %.3f Ang^3' % (v / natoms))
        if a0:
            a = a0 * x
            self.log('Lattice constant: %.3f Ang' % a)
        else:
            a = None
        self.log('Bulk modulus: %.1f GPa' % (B * 1e24 / units.kJ))
        self.log('Total energy: %.3f eV (%d atom%s)' %
                 (e, natoms, ' s'[1:natoms]))

        if plot:
            import pylab as plt
            plt.plot(self.volumes, self.energies, 'o')
            x = np.linspace(self.volumes[0], self.volumes[-1], 50)
            plt.plot(x, eos.fit0(x**-(1.0 / 3)), '-r')
            plt.show()
            
        bulk = self.atoms.copy()
        bulk.set_cell(x * bulk.cell, scale_atoms=True)
        self.write_optimized(bulk, e)

        return e, v, B, a

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
 def eos(self, atoms, name):
     opts = self.opts
     
     traj = PickleTrajectory(self.get_filename(name, 'traj'), 'w', atoms)
     eps = 0.01
     strains = np.linspace(1 - eps, 1 + eps, 5)
     v1 = atoms.get_volume()
     volumes = strains**3 * v1
     energies = []
     cell1 = atoms.cell
     for s in strains:
         atoms.set_cell(cell1 * s, scale_atoms=True)
         energies.append(atoms.get_potential_energy())
         traj.write(atoms)
     traj.close()
     eos = EquationOfState(volumes, energies, opts.eos_type)
     v0, e0, B = eos.fit()
     atoms.set_cell(cell1 * (v0 / v1)**(1 / 3), scale_atoms=True)
     data = {'volumes': volumes,
             'energies': energies,
             'fitted_energy': e0,
             'fitted_volume': v0,
             'bulk_modulus': B,
             'eos_type': opts.eos_type}
     return data

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
    def analyse(self):
        for name, data in self.data.items():
            if 'strains' in data:
                atoms = self.create_system(name)
                # use relaxed volume if present
                if 'relaxed volume' in data:
                    volume = data['relaxed volume']
                else:
                    volume = atoms.get_volume()
                volumes = data['strains']**3 * volume
                energies = data['energies']
                # allow selection of eos type independent of data
                if self.eos is not None:
                    eos = EquationOfState(volumes, energies, self.eos)
                else:
                    eos = EquationOfState(volumes, energies)
                try:
                    v, e, B = eos.fit()
                except (RuntimeError, ValueError):
                    pass
                else:
                    data['fitted energy'] = e
                    data['volume'] = v
                    data['B'] = B

                    if abs(v) < min(volumes) or abs(v) > max(volumes):
                        raise ValueError(name + ': fit outside of range! ' + \
                                         str(abs(v)) + ' not in ' + \
                                         str(volumes))

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
 def f(width, k, g):
     filename = 'Fe-FD-%.2f-%02d-%2d.traj' % (width, k, g)
     configs = read(filename + '@::2')
     # Extract volumes and energies:
     volumes = [a.get_volume() for a in configs]
     energies = [a.get_potential_energy() for a in configs]
     eos = EquationOfState(volumes, energies)
     v0, e0, B = eos.fit()
     return v0, e0, B

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
def relax(input_atoms, ref_db):
    atoms_string = input_atoms.get_chemical_symbols()

    # Open connection to the database with reference data
    db = connect(ref_db)

    # Load our model structure which is just FCC
    atoms = FaceCenteredCubic('X', latticeconstant=1.)
    atoms.set_chemical_symbols(atoms_string)

    # Compute the average lattice constant of the metals in this individual
    # and the sum of energies of the constituent metals in the fcc lattice
    # we will need this for calculating the heat of formation
    a = 0
    ei = 0
    for m in set(atoms_string):
        dct = db.get(metal=m)
        count = atoms_string.count(m)
        a += count * dct.latticeconstant
        ei += count * dct.energy_per_atom
    a /= len(atoms_string)
    atoms.set_cell([a, a, a], scale_atoms=True)

    # Since calculations are extremely fast with EMT we can also do a volume
    # relaxation
    atoms.set_calculator(EMT())
    eps = 0.05
    volumes = (a * np.linspace(1 - eps, 1 + eps, 9))**3
    energies = []
    for v in volumes:
        atoms.set_cell([v**(1. / 3)] * 3, scale_atoms=True)
        energies.append(atoms.get_potential_energy())

    eos = EquationOfState(volumes, energies)
    v1, ef, B = eos.fit()
    latticeconstant = v1**(1. / 3)

    # Calculate the heat of formation by subtracting ef with ei
    hof = (ef - ei) / len(atoms)

    # Place the calculated parameters in the info dictionary of the
    # input_atoms object
    input_atoms.info['key_value_pairs']['hof'] = hof
    # Raw score must always be set
    input_atoms.info['key_value_pairs']['raw_score'] = -hof
    input_atoms.info['key_value_pairs']['latticeconstant'] = latticeconstant

    # Setting the atoms_string directly for easier analysis
    atoms_string = ''.join(input_atoms.get_chemical_symbols())
    input_atoms.info['key_value_pairs']['atoms_string'] = atoms_string

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
 def analyse(self):
     for name, data in self.data.items():
         if 'strains' in data:
             atoms = self.create_system(name)
             volumes = data['strains']**3 * atoms.get_volume()
             energies = data['energies']
             eos = EquationOfState(volumes, energies)
             try:
                 v, e, B = eos.fit()
             except ValueError:
                 pass
             else:
                 data['fitted energy'] = e
                 data['volume'] = v
                 data['B'] = B

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
 def analyse(self):
     OptimizeTask.analyse(self)
     for name, data in self.data.items():
         if 'strains' in data:
             atoms = self.create_system(name)
             volumes = data['strains']**3 * atoms.get_volume()
             energies = data['energies']
             eos = EquationOfState(volumes, energies)
             try:
                 v, e, B = eos.fit()
             except ValueError:
                 self.results[name].extend([None, None])
             else:
                 self.results[name][1:] = [energies[2] - e, v,
                                           B * 1e24 / units.kJ]
         else:
             self.results[name].extend([None, None])

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
def test_calculator():
    """
    Take ASE structure, PySCF object,
    and run through ASE calculator interface. 
    
    This allows other ASE methods to be used with PySCF;
    here we try to compute an equation of state.
    """
    ase_atom=Diamond(symbol='C', latticeconstant=3.5668)

    # Set up a cell; everything except atom; the ASE calculator will
    # set the atom variable
    cell = pbcgto.Cell()
    cell.h=ase_atom.cell
    cell.basis = 'gth-szv'
    cell.pseudo = 'gth-pade'
    cell.gs=np.array([8,8,8])
    cell.verbose = 0

    # Set up the kind of calculation to be done
    # Additional variables for mf_class are passed through mf_dict
    mf_class=pbcdft.RKS
    mf_dict = { 'xc' : 'lda,vwn' }

    # Once this is setup, ASE is used for everything from this point on
    ase_atom.set_calculator(pyscf_ase.PySCF(molcell=cell, mf_class=mf_class, mf_dict=mf_dict))

    print "ASE energy", ase_atom.get_potential_energy()
    print "ASE energy (should avoid re-evaluation)", ase_atom.get_potential_energy()
    # Compute equation of state
    ase_cell=ase_atom.cell
    volumes = []
    energies = []
    for x in np.linspace(0.95, 1.2, 5):
        ase_atom.set_cell(ase_cell * x, scale_atoms = True)
        print "[x: %f, E: %f]" % (x, ase_atom.get_potential_energy())
        volumes.append(ase_atom.get_volume())
        energies.append(ase_atom.get_potential_energy())

    eos = EquationOfState(volumes, energies)
    v0, e0, B = eos.fit()
    print(B / kJ * 1.0e24, 'GPa')
    eos.plot('eos.png')

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
    def analyse(self):
        for name, data in self.data.items():
            if 'strains' in data:
                atoms = self.create_system(name)
                volumes = data['strains']**3 * atoms.get_volume()
                energies = data['energies']
                eos = EquationOfState(volumes, energies)
                try:
                    v, e, B = eos.fit()
                except ValueError:
                    pass
                else:
                    data['fitted energy'] = e
                    data['volume'] = v
                    data['B'] = B

                    if abs(v) < min(volumes) or abs(v) > max(volumes):
                        raise ValueError(name + ': fit outside of range! ' + \
                                         str(abs(v)) + ' not in ' + \
                                         str(volumes))

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
def eosanal(volumes,energies,CI):
    eos=EquationOfState(volumes,energies) #Performs the EOS calcs
    v0,e0,B0=eos.fit() #Gives us the values at minimum energy
    def func(x): #sets up the function to be solved
        return CI-erf(.701707*x)
    #Proposes the function dictating the normal random variable
    z=fsolve(func,0) #solves for the normal random variable
    v=np.power(volumes,-.3333333) #sets up the variable as the one in the EOS
    fit3=np.polyder(((np.poly1d(np.polyfit(v,energies,3)))),2)
    #Creates an equation that solves for the bulk modulus in the same manner as the EOS
    BM=fit3(v)/9*v**5
    #solves for the bulk modulus corresponding to each volume value
    var=[np.std(volumes),np.std(energies),np.std(BM)]
    #solves for the standard deviation of each set of data of interest
    SS=[np.count_nonzero(volumes),np.count_nonzero(energies),np.count_nonzero(BM)]
    #determines the sample size of each set of data of interest.
    R=[z*var[0]*SS[0]**-.5,z*var[1]*SS[1]**-.5,z*var[2]*SS[2]**-.5]
    #Determines the radius of the confidence interval
    Limits=[v0-R[0],v0+R[0],e0-R[1],e0+R[1],B0-R[2],B0+R[2]]
    #sets up each confidence interval
    print 'The {0} confidence interval around the volume at minimum energy of the selected structure is between {1} and {2} A^3.The {0} confidence interval around the minimum energy of the selected structure is between {3} and {4} eV.The {0} confidence interval around the bulk modulus of the selected structure is between {5} and {6} eV/A^3  '.format(CI,Limits[0],Limits[1],Limits[2],Limits[3],Limits[4],Limits[5]) #displays the results
    return;

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
from ase.lattice.cubic import FaceCenteredCubic
from ase.calculators.emt import EMT
from ase.utils.eos import EquationOfState
from ase.db import connect

db = connect('refs.db')

metals = ['Al', 'Au', 'Cu', 'Ag', 'Pd', 'Pt', 'Ni']
for m in metals:
    atoms = FaceCenteredCubic(m)
    atoms.set_calculator(EMT())
    e0 = atoms.get_potential_energy()
    a = atoms.cell[0][0]

    eps = 0.05
    volumes = (a * np.linspace(1 - eps, 1 + eps, 9))**3
    energies = []
    for v in volumes:
        atoms.set_cell([v**(1. / 3)] * 3, scale_atoms=True)
        energies.append(atoms.get_potential_energy())

    eos = EquationOfState(volumes, energies)
    v1, e1, B = eos.fit()

    atoms.set_cell([v1**(1. / 3)] * 3, scale_atoms=True)
    ef = atoms.get_potential_energy()

    db.write(atoms, metal=m,
             latticeconstant=v1**(1. / 3),
             energy_per_atom=ef / len(atoms))

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
def get_eos(self, static=False):
    '''calculate the equation of state for the attached atoms.

    Returns a dictionary of data for each step. You do not need to
    specify any relaxation parameters, only the base parameters for the
    calculations. Writes to eos.org with a report of output.

    if static is True, then run a final static calculation at high
    precision, with ismear=-5.
    '''

    # this returns if the data exists.
    if os.path.exists('eos.json'):
        with open('eos.json') as f:
            return json.loads(f.read())

    # we need an initial calculation to get us going.
    self.calculate()

    cwd = os.getcwd()
    data = {'cwd': os.getcwd()}  # dictionary to store results in

    org = []  # list of strings to make the org-file report
    org += ['#+STARTUP: showeverything']
    org += ['* Initial guess']
    org += [str(self)]
    org += ['',
            '[[shell:jaspsum -p {0}][view initial guess]]'.format(cwd)]

    with open('eos.org', 'w') as f:
        f.write('\n'.join(org))

    atoms = self.get_atoms()
    original_atoms = atoms.copy()  # save for comparison later.
    v_init = atoms.get_volume()

    # ############################################################
    # ## Step 1
    # ############################################################
    org += ['* step 1 - relax ions and shape']
    volumes1, energies1 = [], []
    ready = True
    factors = [-0.15, -0.07, 0.0, 0.07, 0.15]
    for i, f in enumerate(factors):
        wd = cwd + '/step-1/f-{0}'.format(i)
        self.clone(wd)

        with jasp(wd,
                  isif=4,
                  ibrion=2,
                  ediffg=-0.05, ediff=1e-6,
                  nsw=50,
                  atoms=atoms) as calc:
            try:
                # add org-link to calculation
                org += ['[[shell:jaspsum {0}][{0}]]'.format(wd)]

                atoms.set_volume(v_init * (1 + f))
                volumes1.append(atoms.get_volume())
                energies1.append(atoms.get_potential_energy())
                calc.strip()

            except (VaspSubmitted, VaspQueued):
                ready = False

    if not ready:
        log.info('Step 1 is still running')
        raise VaspRunning

    data['step1'] = {}
    data['step1']['volumes'] = volumes1
    data['step1']['energies'] = energies1
    with open('eos.json', 'w') as f:
        f.write(json.dumps(data))

    # create an org-table of the data.
    org += ['',
            '#+tblname: step1',
            '| volume (A^3) | Energy (eV) |',
            '|-']
    for v, e in zip(volumes1, energies1):
        org += ['|{0}|{1}|'.format(v, e)]
    org += ['']

    with open('eos.org', 'w') as f:
        f.write('\n'.join(org))

    eos1 = EquationOfState(volumes1, energies1)

    try:
        v1, e1, B1 = eos1.fit()
    except:
        with open('error', 'w') as f:
            f.write('Error fitting the equation of state')

    data['step1']['eos'] = (v1, e1, B1)
    with open('eos.json', 'w') as f:
        f.write(json.dumps(data))

    # create a plot
#.........这里部分代码省略.........

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
# reoptimize/check volume
#
volumes = []
energies = []
for x in np.linspace(0.98, 1.02, 5):
    atoms.set_cell(init_cell * x, scale_atoms=True)
    volumes.append(atoms.get_volume() / atoms.get_number_of_atoms())
    energies.append(atoms.get_potential_energy() / atoms.get_number_of_atoms())
print "per atom volumes:", volumes
print "per atom energies:", energies

# fit EOS
#
from ase.utils.eos import EquationOfState
eos = EquationOfState(volumes, energies)
vpaf, epaf, B1 = eos.fit()
print "vpaf:", vpaf, "A^3"
print "epaf:", epaf, "eV"
print "B1:", B1 / kJ * 1.0e24, "GPa"

# get optimal lattice parameter from optimal volume
#
volrc = abs(np.linalg.det(refcell))
optlp = pow(vpaf * atoms.get_number_of_atoms() / volrc, 1. / 3.)
print "optlp:", optlp, "\n"

# get actual energy at optimal volume
#
opt_cell = optlp * refcell
atoms.set_cell(opt_cell, scale_atoms=True)
epao = atoms.get_potential_energy() / atoms.get_number_of_atoms()

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import fit [as 别名]
def get_e_v(fname):
    data = np.loadtxt(fname, usecols=(1, 3, 5, 6, 7))
    volumes = data[:, 1]
    energies = data[:, 4]
    return volumes, energies


def custom_plot(volumes, energies, eos):
    plot.plot(volumes, energies, 'ro')
    x = np.linspace(min(eos.v), max(eos.v), 100)
    y = eval(eos.eos_string)(x, eos.eos_parameters[0], eos.eos_parameters[1],
                             eos.eos_parameters[2], eos.eos_parameters[3])
    plot.plot(x, y, label='fit')
    plot.xlabel('Volume ($\AA^3$)')
    plot.ylabel('Energy (eV)')
    plot.legend(loc='best')
    plot.savefig('eos.png')
    plot.show()


if __name__ == '__main__':
    volumes, energies = get_e_v('test_bulk.txt')
    # eos = 'sjeos', 'murnaghan', 'birch', 'taylor', 'vinet' etc.
    eos = EquationOfState(volumes, energies, eos='murnaghan')
    v0, e0, B = eos.fit()
    # the ASE units for the bulk modulus is eV/Angstrom^3
    print('optimum volume, energy and bulk moduls', v0, e0, B)
    # plot
    # eos.plot(filename= "eos_fit")
    custom_plot(volumes, energies, eos)