Python EquationOfState.plot方法代码示例

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [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 plot [as 别名]
            0.02, 0.04, 0.06]
 calculators = [] # reset list
 for f in factors:
     newatoms = atoms.copy()
     newatoms.set_volume(v1*(1 + f))
     label = 'bulk/cu-mp/step2-{0}'.format(COUNTER)
     COUNTER += 1
     calc = jasp(label,
                 xc='PBE',
                 encut=350,
                 kpts=(6,6,6),
                 isym=2,
                 atoms=newatoms)
     calculators.append(calc)
 pool = multiprocessing.Pool(processes=NCORES)
 out = pool.map(do_calculation, calculators)
 pool.close()
 pool.join() # wait here for calculations to finish
 # proceed with analysis
 V += [x[0] for x in out]
 E += [x[1] for x in out]
 V = np.array(V)
 E = np.array(E)
 f = np.array(V)/v1
 # only take points within +- 10% of the minimum
 ind = (f >=0.90) & (f <= 1.1)
 eos = EquationOfState(V[ind], E[ind])
 v2, e2, B = eos.fit()
 print('step2: v2 = {v2}'.format(**locals()))
 eos.plot('images/cu-mp-eos.png')

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]

#.........这里部分代码省略.........
        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
    f = eos1.plot(show=False)
    f.subplots_adjust(left=0.18, right=0.9, top=0.9, bottom=0.15)
    plt.xlabel(u'volume ($\AA^3$)')
    plt.ylabel(u'energy (eV)')
    plt.title(u'E: %.3f eV, V: %.3f $\AA^3$, B: %.3f GPa' %
              (e1, v1, B1 / GPa))

    plt.text(eos1.v0, max(eos1.e), 'EOS: %s' % eos1.eos_string)
    f.savefig('eos-step1.png')

    org += ['[[./eos-step1.png]]',
            '']

    min_energy_index = np.argmin(energies1)

    if min_energy_index in [0, -1]:
        log.warn('Your minimum energy is at an endpoint.'
                 'This indicates something is wrong.')

    with open('eos.org', 'w') as f:
        f.write('\n'.join(org))
    # ########################################################
    # #  STEP 2
    # ########################################################
    # step 2 - isif=4, ibrion=1. now we allow the shape of each cell to
    # change, and we use the best guess from step 1 for minimum volume.
    ready = True
    volumes2, energies2 = [], []
    factors = [-0.09, -0.06, -0.03, 0.0, 0.03, 0.06, 0.09]

    org += ['* step 2 - relax ions and shape with improved minimum estimate']

    for i, f in enumerate(factors):

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]
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)

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]
        v = atoms.get_volume()
        volumes.append(v)
        f.write('{0:1.3f} {1:1.5f} {2:1.5f}'.format(a, v, e))
    except:
        print('aims failed when a = {0:1.3f}'.format(a))
        ready = False
if not ready:
    import sys; sys.exit()


print '#+tblname: pt-bcc-latt'
print r'| lattice constant ($\AA$) | Total Energy (eV) |'
for lc, e in zip(LC, energies):
    print '| {0} | {1} |'.format(lc, e)

import matplotlib.pyplot as plt
plt.plot(LC, energies)
plt.xlabel('Lattice constant ($\AA$)')
plt.ylabel('Total energy (eV)')
plt.savefig('images/pt-bcc-latt.png')

eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print '''
v0 = {0} A^3
E0 = {1} eV
B  = {2} eV/A^3'''.format(v0, e0, B)
f.write('v0 = {0} A^3 \n E0 = {1} eV \n B  = {2} eV/A^3'.format(v0, e0, B))
eos.plot('images/pt-bcc-eos.png')
f.close()

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]
from vasp import Vasp
from ase.utils.eos import EquationOfState
LC = [3.5, 3.55, 3.6, 3.65, 3.7, 3.75]
energies = []
volumes = []
for a in LC:
    calc = Vasp('bulk/Cu-{0}'.format(a))
    atoms = calc.get_atoms()
    volumes.append(atoms.get_volume())
    energies.append(atoms.get_potential_energy())
calc.stop_if(None in energies)
eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print '''
v0 = {0} A^3
E0 = {1} eV
B  = {2} eV/A^3'''.format(v0, e0, B)
eos.plot('images/Cu-fcc-eos.png')

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]
 calculators = [] # reset list
 for f in factors:
     newatoms = atoms.copy()
     newatoms.set_volume(v1*(1 + f))
     label = 'bulk/cu-mp2/step2-{0}'.format(COUNTER)
     COUNTER += 1
     calc = jasp(label,
                 xc='PBE',
                 encut=350,
                 kpts=(6,6,6),
                 isym=2,
                 debug=logging.DEBUG,
                 atoms=newatoms)
     calculators.append(calc)
 pool = multiprocessing.Pool(processes=3)
 out = pool.map(do_calculation, calculators)
 pool.close()
 pool.join() # wait here for calculations to finish
 # proceed with analysis
 V += [x[0] for x in out]
 E += [x[1] for x in out]
 V = np.array(V)
 E = np.array(E)
 f = np.array(V)/v1
 # only take points within +- 10% of the minimum
 ind = (f >=0.90) & (f <= 1.1)
 eos = EquationOfState(V[ind], E[ind])
 v2, e2, B = eos.fit()
 print('step2: v2 = {v2}'.format(**locals()))
 eos.plot('images/cu-mp2-eos.png',show=True)

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]
from ase.io import read 
from ase.units import kJ
from ase.utils.eos import EquationOfState
configs = read('[email protected]:5')  # read 5 configurations
# Extract volumes and energies:
volumes = [ag.get_volume() for ag in configs]
energies = [ag.get_potential_energy() for ag in configs]
eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print B / kJ * 1.0e24, 'GPa'
eos.plot('Ag-eos.png')

# 需要导入模块: from ase.utils.eos import EquationOfState [as 别名]
# 或者: from ase.utils.eos.EquationOfState import plot [as 别名]
cell.verbose = 0

# Set up the kind of calculation to be done
# Additional variables for mf_class are passed through mf_dict
# E.g. gamma-point SCF calculation can be set to
mf_class = pbcdft.RKS
# SCF with k-point sampling can be set to
mf_class = lambda cell: pbcdft.KRKS(cell, kpts=cell.make_kpts([2,2,2]))

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 plot [as 别名]
              spin='none',
              relativistic = 'atomic_zora scalar',
              kpts=(12, 12, 12),  # specifies the Monkhorst-Pack grid
              sc_accuracy_etot=1e-5,
              sc_accuracy_eev=1e-2,
              sc_accuracy_rho=1e-4,
              sc_accuracy_forces=1e-3)
    atoms.set_calculator(calc)
    try:
        e = atoms.get_potential_energy()
        energies.append(e)
        volumes.append(atoms.get_volume())
    except:
        pass
if len(energies) != len(LC):
    import sys; sys.exit()
import matplotlib.pyplot as plt
plt.plot(LC, energies)
plt.xlabel('Lattice constant ($\AA$)')
plt.ylabel('Total energy (eV)')
plt.savefig('images/Pt-fcc.png')

eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print '''
v0 = {0} A^3
E0 = {1} eV
B  = {2} eV/A^3'''.format(v0, e0, B)

eos.plot('images/Pt-fcc-eos.png')