Python eos.EquationOfState类代码示例

 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)

 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

    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

    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))

 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

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

 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

 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])

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')

    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))

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;

    def analyse(self, atomsfile=None):
        try:
            BulkTask.analyse(self)
        except ValueError: # allow fit outside of range
            pass

        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 ValueError:
                    pass
                else:
                    data['fitted energy'] = e
                    data['volume'] = v
                    data['B'] = B
                # with respect tot the reference volume
                data['volume error [%]'] = (data['volume'] / atoms.get_volume() - 1) * 100
                if self.collection.B:
                    i = self.collection.labels.index(self.collection.xc) - 1
                    B = self.collection.B[name][i] * units.kJ * 1e-24
                    data['B error [%]'] = (data['B'] / B - 1) * 100
                else:
                    data['B error [%]'] = None
                data['strukturbericht'] = self.collection.data[name][0]
                data['crystal structure'] = strukturbericht[data['strukturbericht']]
                # calculate lattice constant from volume
                cs = data['crystal structure']
                if cs == 'bcc':
                    a0 = (volume*2)**(1/3.)
                    a = (data['volume']*2)**(1/3.)
                elif cs == 'cesiumchloride':
                    a0 = (volume)**(1/3.)
                    a = (data['volume'])**(1/3.)
                elif cs in ['fcc',
                            'diamond',
                            'zincblende',
                            'rocksalt',
                            'fluorite']:
                    a0 = (volume*4)**(1/3.)
                    a = (data['volume']*4)**(1/3.)
                i = self.collection.labels.index(self.collection.xc) - 1
                a0_ref = self.collection.data[name][i]
                if 'relaxed volume' not in data:
                    # no volume relaxation performed - volume equals the reference one
                    assert abs(a0 - a0_ref) < 1.e-4
                data['lattice constant'] = a
                data['lattice constant error [%]'] = (a - a0_ref) / a0_ref * 100

        if atomsfile:
            # MDTMP: TODO
            atomdata = read_json(atomsfile)
            for name, data in self.data.items():
                atoms = self.create_system(name)
                e = -data['energy']
                for atom in atoms:
                    e += atomdata[atom.symbol]['energy']
                e /= len(atoms)
                data['cohesive energy'] = e
                if self.collection.xc == 'PBE':
                    eref = self.collection.data[name][7]
                else:
                    eref = self.collection.data[name][9]
                data['cohesive energy error [%]'] = (e / eref - 1) * 100

            self.summary_keys += ['cohesive energy',
                                  'cohesive energy error [%]']

from ase.units import kJ
from ase.utils.eos import EquationOfState

lestim = lat0
volumes = []
energies = []
elstr=el1+el2+'-'+str
for x in np.linspace(0.98, 1.02, 5):
    mys.set_cell(lestim * x, scale_atoms=True)
    volumes.append(mys.get_volume()/mys.get_number_of_atoms())
    energies.append(mys.get_potential_energy()/mys.get_number_of_atoms())
print "Volumespa:", volumes
print "Energiespa:", energies

eos = EquationOfState(volumes, energies)
vpa1, epa1, B1 = eos.fit()
lpopt1 = pow(vpa1*(n1+n2)/vpc, 1/3.)

print 'vpa1:', vpa1, 'A^3'
print 'epa1:', epa1, 'eV'
print 'B1:', B1 / kJ * 1.0e24, 'GPa'
print 'lpopt1:', lpopt1
#eos.plot(elstr+'-eos.pdf')#,show=True)

if binary == 1:
    hof = epa1+esub1*n1/(n1+n2)+esub2*n2/(n1+n2)
else:
    hof = epa1/n1-esub1
print "hof1:", hof*1000, "meV/atom\n"

        try:
                e=atoms.get_potential_energy()
                energies.append(e)
                ones.append(1)
        except (VaspSubmitted,VaspQueued):
                ready=False
if not ready:
    import sys; sys.exit()
for m in LC:
    with jasp('bulk/Cl-{0}'.format(m)) as calc:
        atoms=calc.get_atoms()
        volumes.append(atoms.get_volume())

# all of the above sets up the problem for the code below. It simply creates an FCC Cl. Changing the above can easily make a different structure.
print volumes
print energies
eos=EquationOfState(volumes,energies) #Performs the EOS calcs
v0,e0,B0=eos.fit() #Gives us the values at minimum energy
CI=.95 #Determines the level of Confidence Interval
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#solves for the standard deviation of each set of data of interest

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')
    # show()


if __name__ == '__main__':
    # load from file
    # volumes = np.loadtxt('filename')[:,0]
    # energies = np.loadtxt('filename')[:,1]
    # volumes = np.array([13.72, 14.83, 16.0, 17.23, 18.52])
    # energies = np.array([-56.29, -56.41, -56.46, -56.46, -56.42])
    volumes, energies = get_e_v('VOLUME')
    # 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)