Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
task_path = re.sub('confs', global_task_name, conf_path)
task_path = os.path.join(task_path, vasp_str)
os.makedirs(task_path, exist_ok=True)
cwd = os.getcwd()
os.chdir(task_path)
if os.path.isfile('POSCAR') :
os.remove('POSCAR')
os.symlink(os.path.relpath(equi_contcar), 'POSCAR')
os.chdir(cwd)
task_poscar = os.path.join(task_path, 'POSCAR')
# stress
equi_outcar = os.path.join(equi_path, 'OUTCAR')
stress = vasp.get_stress(equi_outcar)
np.savetxt(os.path.join(task_path, 'equi.stress.out'), stress)
# gen strcture
ss = Structure.from_file(task_poscar)
# gen defomations
norm_strains = [-norm_def, -0.5*norm_def, 0.5*norm_def, norm_def]
shear_strains = [-shear_def, -0.5*shear_def, 0.5*shear_def, shear_def]
dfm_ss = DeformedStructureSet(ss,
symmetry = False,
norm_strains = norm_strains,
shear_strains = shear_strains)
n_dfm = len(dfm_ss)
# gen incar
if 'relax_incar' in jdata.keys():
relax_incar_path = jdata['relax_incar']
assert(os.path.exists(relax_incar_path))
relax_incar_path = os.path.abspath(relax_incar_path)
fc = open(relax_incar_path).read()
kspacing =float(re.findall((r"KSPACING(.+?)\n"),fc)[0].replace('=',''))
kgamma =('T' in re.findall((r"KGAMMA(.+?)\n"),fc)[0])
equi_path = os.path.join(equi_path, vasp_str)
equi_contcar = os.path.join(equi_path, 'CONTCAR')
assert os.path.exists(equi_contcar),"Please compute the equilibrium state using vasp first"
task_path = re.sub('confs', global_task_name, conf_path)
task_path = os.path.join(task_path, vasp_str)
os.makedirs(task_path, exist_ok=True)
cwd = os.getcwd()
os.chdir(task_path)
if os.path.isfile('POSCAR') :
os.remove('POSCAR')
os.symlink(os.path.relpath(equi_contcar), 'POSCAR')
os.chdir(cwd)
task_poscar = os.path.join(task_path, 'POSCAR')
# gen strcture
print("task poscar: ", task_poscar)
ss = Structure.from_file(task_poscar)
# gen defects
vds = InterstitialGenerator(ss, insert_ele)
dss = []
for jj in vds :
dss.append(jj.generate_defect_structure(supercell))
# gen incar
if 'relax_incar' in jdata.keys():
relax_incar_path = jdata['relax_incar']
assert(os.path.exists(relax_incar_path))
relax_incar_path = os.path.abspath(relax_incar_path)
fc = open(relax_incar_path).read()
else :
fp_params = jdata['vasp_params']
ecut = fp_params['ecut']
ediff = fp_params['ediff']
npar = fp_params['npar']
def get_spacing(filename='POSCAR', cutoff=0.95):
"""
Returns the interlayer spacing for a 2D material.
"""
structure = Structure.from_file('POSCAR')
lines = open(filename).readlines()
c_axis = lines[4].split()
lattice_parameter = lines[1].split()
split_coords = [line.split() for line in lines[8:8+structure.num_sites]]
z_coords = list()
for coord in split_coords:
z_coord = float(coord[2])
if z_coord > cutoff:
z_coord -= 1
z_coords.append(z_coord)
max_height = max([z_height for z_height in z_coords])
min_height = min([z_height for z_height in z_coords])
spacing = ((1.0 + min_height) - max_height) * float(c_axis[2])\
* float(lattice_parameter[0])
def relax(submit=True, force_overwrite=False):
"""
Should be run before pretty much anything else, in order to get the
right energy of the 2D material.
"""
if force_overwrite or not utl.is_converged(os.getcwd()):
directory = os.getcwd().split('/')[-1]
# Ensure 20A interlayer vacuum
utl.add_vacuum(20 - utl.get_spacing(), 0.9)
# vdw_kernel.bindat file required for VDW calculations.
os.system('cp {} .'.format(KERNEL_PATH))
# KPOINTS
Kpoints.automatic_density(Structure.from_file('POSCAR'),
1000).write_file('KPOINTS')
kpts_lines = open('KPOINTS').readlines()
with open('KPOINTS', 'w') as kpts:
for line in kpts_lines[:3]:
kpts.write(line)
kpts.write(kpts_lines[3].split()[0] + ' '
+ kpts_lines[3].split()[1] + ' 1')
# INCAR
INCAR_DICT.update({'MAGMOM': utl.get_magmom_string()})
Incar.from_dict(INCAR_DICT).write_file('INCAR')
# POTCAR
utl.write_potcar()
# Submission script
if HIPERGATOR == 1:
utl.write_pbs_runjob(directory, 1, 16, '800mb', '6:00:00',
VASP_2D)
def plot_pourbaix_diagram(metastability=0.0, ion_concentration=1e-6, fmt='pdf'):
"""
Creates a Pourbaix diagram for the material in the cwd.
Args:
metastability (float): desired metastable tolerance energy
(meV/atom). <~50 is generally a sensible range to use.
ion_concentration (float): in mol/kg. Sensible values are
generally between 1e-8 and 1.
fmt (str): matplotlib format style. Check the matplotlib
docs for options.
"""
# Create a ComputedEntry object for the 2D material.
composition = Structure.from_file('POSCAR').composition
energy = Vasprun('vasprun.xml').final_energy
cmpd = ComputedEntry(composition, energy)
# Define the chemsys that describes the 2D compound.
chemsys = ['O', 'H'] + [elt.symbol for elt in composition.elements
if elt.symbol not in ['O', 'H']]
# Experimental ionic energies
# See ions.yaml for ion formation energies and references.
exp_dict = ION_DATA['ExpFormEnergy']
ion_correction = ION_DATA['IonCorrection']
# Pick out the ions pertaining to the 2D compound.
ion_dict = dict()
for elt in chemsys:
# bulksDir = os.listdir('inputs/bulks')
# else:
# bulks = []
twod = []
competitors = []
monos = {}
bulks = {}
results = {}
# making direcotries of types of chem sub directories
results['onRatio'] = []
results['hull'] = []
# read in monolayer and bulk motifs
for mon in monosDir:
monos[mon] = Structure.from_file(monolayer_path+'/'+mon)
results[mon] = []
for bulk in bulksDir:
bulks[bulk] = Structure.from_file(bulk_path+'/'+mon)
results[bulk+'_bulk'] = []
# perform the chem sub based on the input X and M lists
for base_name, mono in monos.items():
# os.mkdir(monolayer)
# os.chdir(monolayer)
for M in M_elements:
for X in X_elements:
coreElements = [M,X]
# naming the files .. can be set as poscar.comment
sub_spec = np.unique(mono.species)
if len(sub_spec) != len(coreElements):
print (base_name+' DOES NOT HAVE '+len(coreElements)+' ELEMENTS. ENDING LOOP')
break
# framework and store them in a dictionary ({formula: entry}).
if os.path.isdir('all_competitors'):
os.chdir('all_competitors')
for comp_dir in [
dir for dir in os.listdir(os.getcwd()) if os.path.isdir(dir) and
is_converged(dir)
]:
vasprun = Vasprun('{}/vasprun.xml'.format(comp_dir))
composition = vasprun.final_structure.composition
energy = vasprun.final_energy
finished_competitors[comp_dir] = ComputedEntry(composition, energy)
os.chdir('../')
for directory in directories:
os.chdir(directory)
composition = Structure.from_file('POSCAR').composition
try:
energy = Vasprun('vasprun.xml').final_energy
except:
energy = 100
my_entry = ComputedEntry(composition, energy) # 2D material
entries = MPR.get_entries_in_chemsys(
[elt.symbol for elt in composition]
)
# If the energies of competing phases have been calculated in
# the current framework, put them in the phase diagram instead
# of the MP energies.
for i in range(len(entries)):
formula = entries[i].composition.reduced_formula
if formula in finished_competitors:
entries[i] = finished_competitors[formula]
"""
Convenience method to perform quick loading of data from a filename. The
type of object returned depends the file type.
Args:
fname (string): A filename.
Returns:
Note that fname is matched using unix-style, i.e., fnmatch.
(Structure) if *POSCAR*/*CONTCAR*/*.cif
(Vasprun) *vasprun*
(obj) if *json* (passthrough to monty.serialization.loadfn)
"""
if (fnmatch(fname, "*POSCAR*") or fnmatch(fname, "*CONTCAR*") or
".cif" in fname.lower()) or fnmatch(fname, "*.vasp"):
return Structure.from_file(fname)
elif fnmatch(fname, "*vasprun*"):
from pymatgen.io.vasp import Vasprun
return Vasprun(fname)
elif fnmatch(fname, "*.json*"):
from monty.serialization import loadfn
return loadfn(fname)
def run_task(self, fw_spec):
mode = self["mode"]
disp = self["displacement"]
structure = Structure.from_file("POSCAR")
nm_eigenvecs = np.array(fw_spec["normalmodes"]["eigenvecs"])
nm_norms = np.array(fw_spec["normalmodes"]["norms"])
# displace the sites along the given normal mode
nm_displacement = nm_eigenvecs[mode, :, :] * disp / nm_norms[mode, :, np.newaxis]
for i, vec in enumerate(nm_displacement):
structure.translate_sites(i, vec, frac_coords=False)
# write the modified structure to poscar
structure.to(fmt="poscar", filename="POSCAR")
def plot_local_potential(axis=2, fmt='pdf'):
"""
Plot data from the LOCPOT file along any of the 3 primary axes.
Useful for determining surface dipole moments and electric
potentials on the interior of the material.
"""
ax = plt.figure(figsize=(16, 10)).gca()
locpot = Locpot.from_file('LOCPOT')
structure = Structure.from_file('CONTCAR')
vd = VolumetricData(structure, locpot.data)
abs_potentials = vd.get_average_along_axis(axis)
vacuum_level = max(abs_potentials)
vasprun = Vasprun('vasprun.xml')
cbm = vasprun.get_band_structure().get_cbm()['energy'] - vacuum_level
vbm = vasprun.get_band_structure().get_vbm()['energy'] - vacuum_level
potentials = [potential - vacuum_level for potential in abs_potentials]
axis_length = structure.lattice._lengths[axis]
positions = np.arange(0, axis_length, axis_length / len(potentials))
ax.plot(positions, potentials, linewidth=2, color='k')
ax.set_xlim(0, axis_length)
ax.set_ylim(-20, 0)