import logging
import math
from collections import namedtuple
from datetime import datetime
from functools import wraps, reduce
from multiprocessing import Pool as ProcessPool
from operator import add
from pathlib import Path
from typing import List, Union
from xml.dom import minidom
from xml.dom.minidom import DocumentType, parseString
import numpy as np
from lxml import etree
from pandas import DataFrame
from gvasp.common.base import Atoms, Lattice
from gvasp.common.constant import COLUMNS_32, COLUMNS_8, ORBITALS, RED, RESET, HIGH_SYM
from gvasp.common.descriptor import ValueDescriptor
from gvasp.common.error import StructureNotEqualError, GridNotEqualError, AnimationError, FrequencyError, \
AttributeNotRegisteredError, ParameterError, PotDirNotExistError
from gvasp.common.parameter import Parameter
from gvasp.common.setting import RootDir
from gvasp.common.structure import Structure
from gvasp.common.utils import remove_mapping, is_subset_recommend_pot, str_list
from gvasp.lib import dos_cython, file_bind
POTENTIAL = ['PAW_LDA', 'PAW_PBE', 'PAW_PW91', 'USPP_LDA', 'USPP_PW91']
logger = logging.getLogger(__name__)
[docs]
class StructInfoFile(MetaFile):
def __new__(cls, *args, **kwargs):
if cls is StructInfoFile:
raise TypeError(f"<{cls.__name__} class> may not be instantiated")
return super().__new__(cls)
@property
def structure(self):
return Structure.from_POSCAR(self.name)
[docs]
class CellFile(StructInfoFile):
@property
def structure(self): # overwrite <structure method>
return Structure.from_cell(self.name)
[docs]
def to_POSCAR(self):
self.structure.write_POSCAR(name="POSCAR")
[docs]
class ARCFile(MetaFile):
[docs]
@staticmethod
def write(name: str, structure: List[Structure], lattice: Lattice):
a, b, c = lattice.length
alpha, beta, gamma = lattice.angle
with open(name, "w") as f:
f.write("!BIOSYM archive 3\n")
f.write("PBC=ON\n")
for frame in range(len(structure)):
structure[frame].atoms.set_coord(lattice) # first set frac_coord
structure[frame].atoms.cart_coord = structure[frame].atoms.count * [None] # then clear cart_coord
atoms = structure[frame].atoms.set_coord(Lattice.arc_lattice(lattice)) # set coord from arc_lattice
f.write("Auto Generated CAR File\n")
f.write(f'!DATE {datetime.now().strftime("%a %b %d %H:%M:%S %Y")}\n')
f.write(f"PBC {a:.5f} {b:.5f} {c:.5f} {alpha:.5f} {beta:.5f} {gamma:.5f} (P1)\n")
for atom in atoms:
formula, order = atom.formula, atom.order
element = formula + str(order + 1)
x, y, z = atom.cart_coord
f.write(f"{element:5s} {x:14.10f} {y:14.10f} {z:14.10f} XXXX 1 xx {formula:2s} 0.0000\n")
f.write("end\n")
f.write("end\n")
[docs]
class SubmitFile(MetaFile):
def __init__(self, name):
super(SubmitFile, self).__init__(name=name)
self.title = "AutoGenerated"
self.task, self._task = None, None
self.incar = None
self.kpoints = None
self.constrain = [None, None]
self.head_lines, self.env_lines, self.run_line, self.finish_line = None, None, None, ""
self.pre_process_lines, self.post_process_lines = None, None
self.vasp_std_line, self.vasp_gam_line = None, None
self.submit2write = []
@property
def build(self):
_head, _env = [], []
_vasp_line = ""
_pre_process, _post_process = [], []
for index, line in enumerate(self.strings):
if line.startswith("#SBATCH") or line.startswith("#!"):
if line.startswith("#SBATCH -J"):
_head.append(f"#SBATCH -J {self.title} \n")
else:
_head.append(line)
elif "source" in line or "module" in line:
_env.append(line)
elif "EXEC=" in line:
_vasp_line = line
elif "mpirun" in line:
self.run_line = line
elif "finish" in line:
self.finish_line += line
elif "Pre-Process" in line:
_pre_process.append(index)
elif "Post-Process" in line:
_post_process.append(index)
assert len(_pre_process) == 0 or (len(_pre_process) == 2 and _pre_process[1] > _pre_process[0])
assert len(_post_process) == 0 or (len(_post_process) == 2 and _post_process[1] > _post_process[0])
if "vasp_std" in _vasp_line:
self.vasp_std_line, self.vasp_gam_line = _vasp_line, _vasp_line.replace("vasp_std", "vasp_gam")
elif "vasp_gam" in _vasp_line:
self.vasp_std_line, self.vasp_gam_line = _vasp_line.replace("vasp_gam", "vasp_std"), _vasp_line
self.head_lines = "".join(_head)
self.env_lines = "".join(_env)
if len(_pre_process):
self.pre_process_lines = "".join(self.strings[_pre_process[0]:_pre_process[1] + 1])
else:
self.pre_process_lines = "# User-defined Pre-Process \n# End Pre-Process\n"
if len(_post_process):
self.post_process_lines = "".join(self.strings[_post_process[0]:_post_process[1] + 1])
else:
self.post_process_lines = "# User-defined Post-Process \n# End Post-Process\n"
return self
@property
def check_success_lines(self):
_lines = ["success=`grep accuracy OUTCAR | wc -l` \n",
"if [ $success -ne 1 ];then \n",
f" echo '{self.task} Task Failed!' \n",
" exit 1 \n",
"fi \n"]
return "".join(_lines)
@property
def backup_lines(self):
_lines = ["cp POSCAR POSCAR_300 \n",
"cp CONTCAR CONTCAR_300 \n",
"cp OUTCAR OUTCAR_300 \n",
"mv CONTCAR POSCAR \n"]
if self.task == "ConTS":
_lines.insert(2, "cp fort.188 fort.188_300 \n")
return "".join(_lines)
@property
def modify_lines(self):
_lines = [f"sed -i 's/ENCUT = 300.0/ENCUT = {self.incar._ENCUT}/' INCAR\n",
f"sed -i 's/PREC = Low/PREC = {self.incar._PREC}/' INCAR\n",
f"sed -i 's/1 1 1/{str_list(self.kpoints.number)}/' KPOINTS\n"]
if self.task == "ConTS":
_lines.insert(0, "dis=`grep 'distance after opt' OUTCAR | tail -1 | awk '{print $NF}'` \n")
_lines.insert(-1, f'sed -i "6c\{self.constrain[0]} {self.constrain[1]} $dis" fort.188\n')
return "".join(_lines)
@property
def bader_lines(self):
_lines = ["gvasp sum || chgsum.pl AECCAR0 AECCAR2 || return 1 \n",
"bader CHGCAR -ref CHGCAR_sum \n"]
return "".join(_lines)
@property
def spin_lines(self):
_lines = ["gvasp split || chgsplit.pl CHGCAR || return 1 \n",
"gvasp grd -d -1 || return 1 \n"]
return "".join(_lines)
[docs]
def pipe(self, attrs: list) -> str:
_cfg = []
for attr in attrs:
try:
_cfg.append(getattr(self, attr))
except AttributeError:
_cfg.append(attr)
return "".join(_cfg)
[docs]
class Fort188File(MetaFile):
@property
def constrain(self):
return self.strings[5].split()
@constrain.setter
def constrain(self, info: list):
self.strings[5] = " ".join(info) + "\n"
[docs]
def write(self, name="fort.188"):
with open(name, "w") as f:
f.writelines(self.strings)
[docs]
class INCAR(MetaFile, Parameter):
def __init__(self, name):
super(INCAR, self).__init__(name=name)
self._init_attr()
def __getattr__(self, item):
if item in self.__dict__.keys():
return self.__dict__[item]
elif item in self.__class__.__dict__.keys():
return self.__class__.__dict__[item]
def _init_attr(self):
"""
Initialize the attributes, if attribute not found, raise AttributeNotRegisteredError
"""
for line in self.strings:
if not line.strip(): # catch blanklines
continue
if line.split()[0][0] != "#":
real_line = line.split("#")[0]
attr_name = real_line.split("=")[0].strip()
attr_value = real_line.split("=")[1].strip()
for param_type in INCAR._type_trans.keys():
if attr_name in INCAR._type_trans[param_type]['name']:
attr_value = INCAR._type_trans[param_type]['func'](attr_value)
break
else:
raise AttributeNotRegisteredError(f"{attr_name} is not registered")
setattr(self, attr_name, attr_value)
[docs]
def write(self, name):
"""
Write interface, callback every _write* func
"""
self._write_base(name)
self._write_scf(name)
self._write_hse(name)
self._write_vdw(name)
self._write_dipole(name)
self._write_sol(name)
self._write_opt(name)
self._write_md(name)
self._write_charge(name)
self._write_wf(name)
self._write_density(name)
self._write_freq(name)
self._write_stm(name)
self._write_vtst(name)
self._write_neb(name)
self._write_dimer(name)
if self.LHFCALC is not True:
self._write_plusU(name)
@formatter(Parameter._baseParam)
def _write_base(self, name):
with open(name, "w") as f:
f.write(f"#----------/Base Parameter/----------# \n")
@formatter(Parameter._scfParam)
def _write_scf(self, name):
with open(name, "a+") as f:
f.write(f"#----------/SCF Parameter/----------# \n")
@formatter(Parameter._hseParam)
def _write_hse(self, name):
with open(name, "a+") as f:
f.write(f"#----------/HSE06 Parameter/----------# \n")
@formatter(Parameter._vdwParam)
def _write_vdw(self, name):
with open(name, "a+") as f:
f.write(f"#----------/VDW Parameter/----------# \n")
@formatter(Parameter._dipoleParam)
def _write_dipole(self, name):
with open(name, "a+") as f:
f.write(f"#----------/Dipole Parameter/----------# \n")
@formatter(Parameter._solParam)
def _write_sol(self, name):
with open(name, "a+") as f:
f.write(f"#----------/Solvation Parameter/----------# \n")
@formatter(Parameter._optParam)
def _write_opt(self, name):
with open(name, "a+") as f:
f.write(f"#----------/Optimize Parameter/----------# \n")
@formatter(Parameter._mdParam)
def _write_md(self, name):
with open(name, "a+") as f:
f.write(f"#----------/MD Parameter/----------# \n")
@formatter(Parameter._chgParam)
def _write_charge(self, name):
with open(name, "a+") as f:
f.write(f"#----------/Charge Parameter/----------# \n")
@formatter(Parameter._wfParam)
def _write_wf(self, name):
with open(name, "a+") as f:
f.write(f"#----------/WorkFunc Parameter/----------# \n")
@formatter(Parameter._dosParam)
def _write_density(self, name):
with open(name, "a+") as f:
f.write(f"#----------/DOS Parameter/----------# \n")
@formatter(Parameter._freqParam)
def _write_freq(self, name):
with open(name, "a+") as f:
f.write(f"#----------/Frequency Parameter/----------# \n")
@formatter(Parameter._stmParam)
def _write_stm(self, name):
with open(name, "a+") as f:
f.write(f"#----------/STM Parameter/----------# \n")
@formatter(Parameter._vtstParam)
def _write_vtst(self, name):
with open(name, "a+") as f:
f.write(f"#----------/VTST optimizer Parameter/----------# \n")
@formatter(Parameter._nebParam)
def _write_neb(self, name):
with open(name, "a+") as f:
f.write(f"#----------/NEB Parameter/----------# \n")
@formatter(Parameter._dimerParam)
def _write_dimer(self, name):
with open(name, "a+") as f:
f.write(f"#----------/Dimer Parameter/----------# \n")
@formatter(Parameter._plusUParam)
def _write_plusU(self, name):
with open(name, "a+") as f:
f.write(f"#----------/+U Parameter/----------# \n")
[docs]
class KPOINTS(MetaFile):
def __init__(self, name, task='opt', parse=True):
super(KPOINTS, self).__init__(name=name)
self.task = task
self.title, self.strategy, self.center, self.number, self.weight = None, None, None, None, None
self._parse() if parse else None
def _parse(self):
if self.task in ['opt']:
self.title = self.strings[0].strip()
self.strategy = self.strings[1].strip()
self.center = self.strings[2].strip()
self.number = list(map(int, self.strings[3].split()))
self.weight = list(map(float, self.strings[4].split()))
[docs]
def write(self, name):
with open(name, 'w') as f:
f.write(self.title + " \n")
f.write(self.strategy + " \n")
f.write(self.center + " \n")
f.write(" ".join(map(str, self.number)) + " \n")
f.write(" ".join(map(str, self.weight)) + " \n")
[docs]
@staticmethod
def min_number(structure: Structure, length=20.0):
lattice = structure.lattice
kpoints = np.ceil(length / lattice.length).astype(int)
structure.atoms.set_coord(lattice)
z_max = np.max(structure.atoms.cart_coord[:, 2])
z_min = np.min(structure.atoms.cart_coord[:, 2])
vacuum = lattice.length[2] - (z_max - z_min)
if (vacuum) >= 5:
kpoints[2] = 1
return kpoints
[docs]
@staticmethod
def from_strings(strings):
_kpoints = KPOINTS(name="KPOINTS", parse=False)
_kpoints._strings = strings
_kpoints._parse()
return _kpoints
[docs]
class POTCAR(MetaFile):
def __init__(self, name):
super(POTCAR, self).__init__(name=name)
self.potential, self.element, self.valence = None, None, None
if Path(self.name).exists():
self._parse()
def __add__(self, other):
potential = self.potential if isinstance(self.potential, list) else [self.potential]
element = self.element if isinstance(self.element, list) else [self.element]
valence = self.valence if isinstance(self.valence, list) else [self.valence]
strings = self._strings
potential += other.potential
element += other.element
valence += other.valence
strings += other._strings
potcar = POTCAR("POTCAR")
potcar.potential = potential
potcar.element = element
potcar.valence = valence
potcar._strings = strings
return potcar
def _parse(self):
self.potential = [line.split()[2] for line in self.strings if line.find("TITEL") != -1]
self.element = [line.split()[3] for line in self.strings if line.find("TITEL") != -1]
self.valence = [float(line.split()[5]) for line in self.strings if line.find("ZVAL") != -1]
[docs]
@staticmethod
def cat(potentials, elements: List[str], potdir=f"{RootDir}/pot"):
if potdir is None or not Path(potdir).exists():
raise PotDirNotExistError(f"potdir = {potdir} is not exist (should named as `pot`), please check")
if (isinstance(potentials, str) and potentials not in POTENTIAL) or \
(isinstance(potentials, list) and len(set(potentials).difference(set(POTENTIAL)))):
raise TypeError(f"potentials should be {POTENTIAL}")
if isinstance(potentials, str):
potentials = [potentials] * len(elements)
elif isinstance(potentials, list) and len(potentials) == len(elements):
pass
else:
raise ParameterError(f"{potentials} and {elements} is not match")
potcar_objects = []
for potential, element in zip(potentials, elements):
marker = is_subset_recommend_pot(element)
if potential == "PAW_PBE" and marker:
potcar_objects.append(POTCAR(name=(Path(potdir) / potential / f"{Path(element)}_{marker}" / "POTCAR")))
logger.warning(
f"VASP recommend `{element} ({potential})` use *_{marker} potential "
f"[https://www.vasp.at/wiki/index.php/Available_PAW_potentials#Recommended_potentials_for_DFT_calculationsss]")
else:
potcar_objects.append(POTCAR(name=(Path(potdir) / potential / f"{Path(element)}" / "POTCAR")))
return reduce(add, potcar_objects)
[docs]
def write(self, name):
with open(name, "w") as f:
f.writelines(self.strings)
[docs]
class XSDFile(MetaFile):
def __init__(self, name):
super(XSDFile, self).__init__(name=name)
self._xml = None
self.element, self.spin, self.frac_coord, self.selective_dynamics, self.lattice = None, None, None, None, None
self.constrain = None
self._parse()
@property
def strings(self):
if self._strings is None:
with open(self.name, "r") as f:
self._strings = f.read()
return self._strings
def _parse(self):
self._xml = etree.XML(self.strings.encode("utf-8"))
try:
SpaceGroupName = self._xml.xpath("//SpaceGroup//@Name")[0] # Cu2O xsd file (why generate no reason)
except IndexError:
SpaceGroupName = self._xml.xpath("//SpaceGroup//@GroupName")[0]
if SpaceGroupName == "P1":
Atom3d = self._xml.xpath("//Atom3d[@Components]")
IdentifyMapping = self._xml.xpath("//MappingFamily//IdentityMapping")[0]
identify_coord = [float(value) for index, value in
enumerate(IdentifyMapping.attrib['Constraint'].split(","))
if (index + 1) % 4 == 0]
Components = [item.split(",")[0] for item in self._xml.xpath("//Atom3d//@Components")]
Name = [atom.attrib.get('Name', atom.attrib['Components']) for atom in Atom3d]
FormalSpin = [int(atom.attrib.get('FormalSpin', '0')) for atom in Atom3d]
XYZ = []
for atom in Atom3d:
item = atom.attrib.get('XYZ', None)
if item is not None:
XYZ.append(list(map(float, item.split(","))))
else:
XYZ.append(identify_coord)
RestrictedProperties = [atom.attrib.get('RestrictedProperties', 'T T T') for atom in Atom3d]
TF = [item.replace("FractionalXYZ", "F F F").split() for item in RestrictedProperties]
else:
Atom3d = self._xml.xpath("//MappingFamily//Atom3d")
IdentifyMapping = self._xml.xpath("//MappingFamily//Atom3d//..")
atoms = []
for atom, identify in zip(Atom3d, IdentifyMapping):
if atom.attrib.get('Components', None) is not None:
formula = atom.attrib['Components']
formal_spin = int(atom.attrib.get('FormalSpin', '0'))
name = atom.attrib.get('Name', atom.attrib['Components'])
restricted_property = atom.attrib.get('RestrictedProperties', 'T T T')
tf = restricted_property.replace("FractionalXYZ", "F F F").split()
identify_coord = [float(value) for index, value in enumerate(identify.attrib['Constraint'].split(","))
if (index + 1) % 4 == 0]
atoms.append((formula, identify_coord, formal_spin, name, tf))
identity_atoms = remove_mapping(atoms)
Components = [atom[0] for atom in identity_atoms]
XYZ = [atom[1] for atom in identity_atoms]
FormalSpin = [atom[2] for atom in identity_atoms]
Name = [atom[3] for atom in identity_atoms]
TF = [atom[4] for atom in identity_atoms]
assert len(XYZ) == len(Components) == len(Name) == len(FormalSpin) == len(TF)
constrain = [index for index, name in enumerate(Name) if name.endswith("_c")]
SpaceGroup = self._xml.xpath("//SpaceGroup")[0]
Vector = [list(map(float, SpaceGroup.attrib[key].split(","))) for key in SpaceGroup.keys() if "Vector" in key]
self.element = Components
self.constrain = constrain
self.spin = FormalSpin
self.frac_coord = XYZ
self.selective_dynamics = TF
self.lattice = Vector
@property
def structure(self):
results = sorted(zip(range(len(self.element)), self.element, self.frac_coord), key=lambda x: (x[1], x[2][2]))
sorted_order, sorted_element, sorted_frac_coord = list(zip(*results))
sorted_constrain = list(map(lambda x: True if x in self.constrain else False, sorted_order))
sorted_selective_dynamics = np.array(self.selective_dynamics)[list(sorted_order)]
sorted_spin = np.array(self.spin)[list(sorted_order)]
atoms = Atoms(formula=sorted_element, frac_coord=sorted_frac_coord, selective_matrix=sorted_selective_dynamics,
constrain=sorted_constrain, spin=sorted_spin)
lattice = Lattice(np.array(self.lattice))
return Structure(atoms=atoms, lattice=lattice)
[docs]
@staticmethod
def write(contcar: Union[str, Path], outcar: Union[str, Path], name='output.xsd'):
structure = CONTCAR(contcar).structure
outcar = OUTCAR(outcar)
doc = minidom.Document()
doctype = DocumentType(qualifiedName="XSD []")
XSD = doc.createElement('XSD')
XSD.setAttribute('Version', "5.0")
XSD.setAttribute('WrittenBy', "Materials Studio 5.0")
AtomisticTreeRoot = doc.createElement('AtomisticTreeRoot')
AtomisticTreeRoot.setAttribute('ID', "1")
AtomisticTreeRoot.setAttribute('NumProperties', "58")
AtomisticTreeRoot.setAttribute('NumChildren', "1")
PropertyString = '<root>' \
'<Property Name="AngleAxisType" DefinedOn="AngleBetweenPlanesBender" Type="Enumerated"/>' \
'<Property Name="AngleEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="BeadDocumentID" DefinedOn="MesoMoleculeSet" Type="String"/>' \
'<Property Name="BendBendEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="BendTorsionBendEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="BondEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="EFGAsymmetry" DefinedOn="Atom" Type="Double"/>' \
'<Property Name="EFGQuadrupolarCoupling" DefinedOn="Atom" Type="Double"/>' \
'<Property Name="ElectrostaticEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="FaceMillerIndex" DefinedOn="GrowthFace" Type="MillerIndex"/>' \
'<Property Name="FacetTransparency" DefinedOn="GrowthFace" Type="Float"/>' \
'<Property Name="FermiLevel" DefinedOn="ScalarFieldBase" Type="Double"/>' \
'<Property Name="Force" DefinedOn="Matter" Type="CoDirection"/>' \
'<Property Name="FrameFilter" DefinedOn="Trajectory" Type="String"/>' \
'<Property Name="HarmonicForceConstant" DefinedOn="HarmonicRestraint" Type="Double"/>' \
'<Property Name="HarmonicMinimum" DefinedOn="HarmonicRestraint" Type="Double"/>' \
'<Property Name="HydrogenBondEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="ImportOrder" DefinedOn="Bondable" Type="UnsignedInteger"/>' \
'<Property Name="InversionEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="IsBackboneAtom" DefinedOn="Atom" Type="Boolean"/>' \
'<Property Name="IsChiralCenter" DefinedOn="Atom" Type="Boolean"/>' \
'<Property Name="IsOutOfPlane" DefinedOn="Atom" Type="Boolean"/>' \
'<Property Name="IsRepeatArrowVisible" DefinedOn="ElectrodeWire" Type="Boolean"/>' \
'<Property Name="KineticEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="LineExtentPadding" DefinedOn="BestFitLineMonitor" Type="Double"/>' \
'<Property Name="LinkageGroupName" DefinedOn="Linkage" Type="String"/>' \
'<Property Name="ListIdentifier" DefinedOn="PropertyList" Type="String"/>' \
'<Property Name="NMRShielding" DefinedOn="Atom" Type="Double"/>' \
'<Property Name="NonBondEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="NormalMode" DefinedOn="Bondable" Type="Direction"/>' \
'<Property Name="NormalModeFrequency" DefinedOn="Bondable" Type="Double"/>' \
'<Property Name="NumScanSteps" DefinedOn="LinearScan" Type="UnsignedInteger"/>' \
'<Property Name="OrbitalCutoffRadius" DefinedOn="Bondable" Type="Double"/>' \
'<Property Name="PlaneExtentPadding" DefinedOn="BestFitPlaneMonitor" Type="Double"/>' \
'<Property Name="PotentialEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="QuantizationValue" DefinedOn="ScalarFieldBase" Type="Double"/>' \
'<Property Name="RelativeVelocity" DefinedOn="Matter" Type="Direction"/>' \
'<Property Name="RepeatArrowScale" DefinedOn="ElectrodeWire" Type="Float"/>' \
'<Property Name="RestraintEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="ScanEnd" DefinedOn="LinearScan" Type="Double"/>' \
'<Property Name="ScanStart" DefinedOn="LinearScan" Type="Double"/>' \
'<Property Name="SeparatedStretchStretchEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="SimulationStep" DefinedOn="Trajectory" Type="Integer"/>' \
'<Property Name="StretchBendStretchEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="StretchStretchEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="StretchTorsionStretchEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="Temperature" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="ThreeBodyNonBondEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="TorsionBendBendEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="TorsionEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="TorsionStretchEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="TotalEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="Units" DefinedOn="ScalarFieldBase" Type="String"/>' \
'<Property Name="ValenceCrossTermEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="ValenceDiagonalEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="VanDerWaalsEnergy" DefinedOn="ClassicalEnergyHolder" Type="Double"/>' \
'<Property Name="_Stress" DefinedOn="MatterSymmetrySystem" Type="Matrix"/>' \
'<Property Name="_TrajectoryStress" DefinedOn="MatterSymmetrySystem" Type="Matrix"/>' \
'</root>'
Propertys = parseString(PropertyString).getElementsByTagName("Property")
for Property in Propertys:
AtomisticTreeRoot.appendChild(Property)
atoms_num = structure.atoms.count
energy = outcar.last_energy
force = outcar.last_force if outcar.last_force is not None else 0
mag = outcar.last_mag if outcar.last_force is not None else 0
SymmetrySystem = doc.createElement("SymmetrySystem")
SymmetrySystem.setAttribute("ID", "2")
SymmetrySystem.setAttribute("Mapping", "3")
SymmetrySystem.setAttribute("Children", f"ci({atoms_num + 2}):{atoms_num + 2}+4")
SymmetrySystem.setAttribute("Normalized", "1")
SymmetrySystem.setAttribute("Name", f"E:{energy:.3f} F:{force:.3f} M:{mag:.2f}")
SymmetrySystem.setAttribute("XYZ", "0,0,0")
SymmetrySystem.setAttribute("OverspecificationTolerance", "0.05")
SymmetrySystem.setAttribute("PeriodicDisplayType", "Original")
SymmetrySystem.setAttribute("HasSymmetryHistory", "1")
MappingSet = doc.createElement("MappingSet")
MappingSet.setAttribute("ID", f"{atoms_num + 6}")
MappingSet.setAttribute("SymmetryDefinition", f"{atoms_num + 4}")
MappingSet.setAttribute("ActiveSystem", "2")
MappingSet.setAttribute("NumFamilies", "1")
MappingSet.setAttribute("OwnsTotalConstraintMapping", "1")
MappingSet.setAttribute("TotalConstraintMapping", "3")
MappingFamily = doc.createElement("MappingFamily")
MappingFamily.setAttribute("ID", f"{atoms_num + 7}")
MappingFamily.setAttribute("NumImageMappings", "0")
IdentityMapping = doc.createElement("IdentityMapping")
IdentityMapping.setAttribute("ID", f"{atoms_num + 8}")
IdentityMapping.setAttribute("Element", "1,0,0,0,0,1,0,0,0,0,1,0")
IdentityMapping.setAttribute("Constraint", "1,0,0,0,0,1,0,0,0,0,1,0")
IdentityMapping.setAttribute("MappedObjects", f"ci({atoms_num + 1}):{atoms_num}+4,{atoms_num + 5}")
IdentityMapping.setAttribute("DefectObjects", f"{atoms_num + 4},{atoms_num + 9}")
IdentityMapping.setAttribute("NumImages", f"{atoms_num + 1}")
IdentityMapping.setAttribute("NumDefects", "2")
for atom in structure.atoms:
Atom3d = doc.createElement("Atom3d")
Atom3d.setAttribute("ID", f"{atom.order + 4}")
Atom3d.setAttribute("Mapping", f"{atoms_num + 8}")
Atom3d.setAttribute("Parent", "2")
if (atom.selective_matrix == ['F', 'F', 'F']).all():
Atom3d.setAttribute("RestrictedBy", f"{atoms_num + 5}")
Atom3d.setAttribute("RestrictedProperties", "FractionalXYZ")
Atom3d.setAttribute("Parent", "2")
Atom3d.setAttribute("Name", f"{atom.formula}{atom.order + 1}")
Atom3d.setAttribute("XYZ", f"{atom.frac_coord[0]:.12f},{atom.frac_coord[1]:.12f},{atom.frac_coord[2]:.12f}")
Atom3d.setAttribute("Charge", "0")
Atom3d.setAttribute("Components", f"{atom.formula}")
Atom3d.setAttribute("FormalSpin", "0")
IdentityMapping.appendChild(Atom3d)
CompleteRestriction = doc.createElement("CompleteRestriction")
CompleteRestriction.setAttribute("ID", f"{atoms_num + 5}")
CompleteRestriction.setAttribute("Mapping", f"{atoms_num + 8}")
CompleteRestriction.setAttribute("Parent", "2")
fix_atoms = [atom.order + 4 for atom in structure.atoms if (atom.selective_matrix == ['F', 'F', 'F']).all()]
if len(fix_atoms):
CompleteRestriction.setAttribute("RestrictsObjects",
f"ci({len(fix_atoms)}):{','.join([str(item) for item in fix_atoms])}")
CompleteRestriction.setAttribute("RestrictsProperties",
f"{','.join(['FractionalXYZ'] * len(fix_atoms))}")
else:
CompleteRestriction.setAttribute("Visible", "0")
matrix = structure.lattice.matrix
SpaceGroup = doc.createElement("SpaceGroup")
SpaceGroup.setAttribute("ID", f"{atoms_num + 4}")
SpaceGroup.setAttribute("Parent", "2")
SpaceGroup.setAttribute("Children", f"{atoms_num + 9}")
SpaceGroup.setAttribute("Name", "P1")
SpaceGroup.setAttribute("AVector", f"{matrix[0, 0]:.10f},{matrix[0, 1]:.10f},{matrix[0, 2]:.10f}")
SpaceGroup.setAttribute("BVector", f"{matrix[1, 0]:.10f},{matrix[1, 1]:.10f},{matrix[1, 2]:.10f}")
SpaceGroup.setAttribute("CVector", f"{matrix[2, 0]:.10f},{matrix[2, 1]:.10f},{matrix[2, 2]:.10f}")
SpaceGroup.setAttribute("Color", "136,0,204,255")
SpaceGroup.setAttribute("OrientationBase", "C along Z, A in XZ plane")
SpaceGroup.setAttribute("Centering", "3D Primitive-Centered")
SpaceGroup.setAttribute("Lattice", "3D Triclinic")
SpaceGroup.setAttribute("GroupName", "P1")
SpaceGroup.setAttribute("Operators", "1,0,0,0,0,1,0,0,0,0,1,0")
SpaceGroup.setAttribute("DisplayRange", "0,1,0,1,0,1")
SpaceGroup.setAttribute("CylinderRadius", "0.2")
SpaceGroup.setAttribute("LabelAxes", "1")
SpaceGroup.setAttribute("ActiveSystem", "2")
SpaceGroup.setAttribute("ITNumber", "1")
SpaceGroup.setAttribute("LongName", "P 1")
SpaceGroup.setAttribute("Qualifier", "Origin-1")
SpaceGroup.setAttribute("SchoenfliesName", "C1-1")
SpaceGroup.setAttribute("System", "Triclinic")
SpaceGroup.setAttribute("Class", "1")
SpaceGroup.setAttribute("DisplayStyle", "Dashed")
SpaceGroup.setAttribute("LineThickness", "4")
ReciprocalLattice3D = doc.createElement("ReciprocalLattice3D")
ReciprocalLattice3D.setAttribute("ID", f"{atoms_num + 9}")
ReciprocalLattice3D.setAttribute("Parent", f"{atoms_num + 4}")
ReciprocalLattice3D.setAttribute("Visible", "0")
IdentityMapping.appendChild(CompleteRestriction)
IdentityMapping.appendChild(SpaceGroup)
IdentityMapping.appendChild(ReciprocalLattice3D)
MappingRepairs = doc.createElement("MappingRepairs")
MappingRepairs.setAttribute("NumRepairs", "0")
MappingFamily.appendChild(IdentityMapping)
MappingFamily.appendChild(MappingRepairs)
InfiniteMapping = doc.createElement("InfiniteMapping")
InfiniteMapping.setAttribute("ID", "3")
InfiniteMapping.setAttribute("Element", "1,0,0,0,0,1,0,0,0,0,1,0")
InfiniteMapping.setAttribute("MappedObjects", "2")
MappingSet.appendChild(MappingFamily)
MappingSet.appendChild(InfiniteMapping)
OriginalObjects = doc.createElement("OriginalObjects")
OriginalObjects.setAttribute("ID", f"{atoms_num + 10}")
SetCollection = doc.createElement("SetCollection")
SetCollection.setAttribute("Objects", f"ci({atoms_num}):{atoms_num}+4")
OriginalObjects.appendChild(SetCollection)
SymmetrySystem.appendChild(MappingSet)
SymmetrySystem.appendChild(OriginalObjects)
AtomisticTreeRoot.appendChild(SymmetrySystem)
XSD.appendChild(AtomisticTreeRoot)
doc.appendChild(doctype)
doc.appendChild(XSD)
stem = Path(name).stem
stem = stem + "-y" if outcar.finish else stem + "-n"
with open(f"{stem}.xsd", "w") as f:
doc.writexml(f, indent='\t', addindent='\t', newl='\n', encoding="latin1")
logger.info(f"{RED}transform to {stem}.xsd file{RESET}")
[docs]
class POSCAR(StructInfoFile):
[docs]
@staticmethod
def dist(pos1: str, pos2: str):
"""
Calculate the distance of two POSCAR, distance = sqrt(sum((i-j)**2))
@param:
pos1: first POSCAR file name
pos2: second POSCAR file name
"""
structure1 = POSCAR(name=pos1).structure
structure2 = POSCAR(name=pos2).structure
return Structure.dist(structure1, structure2)
[docs]
@staticmethod
def align(pos1: str, pos2: str):
"""
Tailor the atoms' order to make the distance of pos1 and pos2 minimum
@param:
pos1: first POSCAR file name
pos2: second POSCAR file name
"""
structure1 = POSCAR(name=pos1).structure
structure2 = POSCAR(name=pos2).structure
logger.info(f"Align before: dist = {Structure.dist(structure1, structure2)}")
structure1_new, structure2_new = Structure.align(structure1, structure2)
logger.info(f"Align before: dist = {Structure.dist(structure1_new, structure2_new)}")
structure1_new.write_POSCAR(f"{pos1}_sort")
structure2_new.write_POSCAR(f"{pos2}_sort")
logger.info(f"New structure have been written to *_sort files")
[docs]
class CONTCAR(POSCAR):
pass
[docs]
class XDATCAR(StructInfoFile):
def __init__(self, name):
super().__init__(name)
self.lattice = Lattice.from_string(self.strings[2:5])
element_name = self.strings[5].split()
element_count = [int(item) for item in self.strings[6].split()]
self.element = sum([[name] * count for name, count in zip(element_name, element_count)], [])
self.frames = [i for i in range(len(self.strings)) if self.strings[i].find("Direct") != -1]
self._structure = []
@property
def structure(self): # overwrite <structure method>
if len(self._structure) == 0:
for frame in self.frames:
frac_coord = np.array([[float(item) for item in line.split()] for line in
self.strings[frame + 1:frame + 1 + len(self.element)]])
atoms = Atoms(formula=self.element, frac_coord=frac_coord)
self._structure.append(Structure(atoms=atoms, lattice=self.lattice))
return self._structure
[docs]
def movie(self, name):
"""Transform the XDATCAR to arc file"""
ARCFile.write(name=name, structure=self.structure, lattice=self.lattice)
[docs]
class DOSCAR(MetaFile):
ISPIN = ValueDescriptor("ISPIN", [1, 2])
LORBIT = ValueDescriptor("LORBIT", [0, 1, 2, 5, 10, 11, 12, 13, 14])
def __init__(self, name, ISPIN=2, LORBIT=12):
super(DOSCAR, self).__init__(name=name)
self.ISPIN = ISPIN
self.LORBIT = LORBIT
self.NAtom = int(self.strings[0].split()[0])
self.Emax, self.Emin, self.NDOS, self.fermi = tuple(map(float, self.strings[5].split()[:4]))
self.NDOS = int(self.NDOS)
if self.LORBIT not in [10, 12]:
logger.error(f"LORBIT = {self.LORBIT} is not supported in this version!")
exit(1)
[docs]
def load(self):
def merge_dos(energy_list, Total_up, Total_down, atom_list, length):
"""TODO: need optimize, deprecate DataFrame"""
atom_data = [energy_list]
columns, orbitals = [], []
Total_Dos = DataFrame(index=energy_list, columns=['tot_up', 'tot_down'], dtype='object')
Total_Dos['tot_up'] = Total_up
Total_Dos['tot_down'] = Total_down if self.ISPIN == 2 else [0.0] * len(Total_Dos['tot_up'])
if self.LORBIT == 10:
columns = COLUMNS_8[:length] if self.ISPIN == 2 else COLUMNS_8[:length * 2]
orbitals = []
elif self.LORBIT == 12:
columns = COLUMNS_32[:length] if self.ISPIN == 2 else COLUMNS_32[:length * 2]
orbitals = ORBITALS[1:int(math.sqrt(length / 2))] if self.ISPIN == 2 else ORBITALS[
1:int(math.sqrt(length))]
for data in atom_list:
data_ispin = np.zeros(shape=(len(energy_list), len(columns)))
if self.ISPIN == 2:
data_ispin = np.array(data)
else:
data_ispin[:, ::2] += np.array(data)
DATA = DataFrame(data_ispin, index=energy_list, columns=columns)
DATA['up'] = 0.0
DATA['down'] = 0.0
# LORBIT = 12, sum projected-orbitals, e.g., px+py+pz = p
for orbital in orbitals:
DATA[orbital + '_up'] = 0.0
DATA[orbital + '_down'] = 0.0
orbital_p_up = [item for item in DATA.columns.values if
item.startswith(orbital) and item.endswith('up') and item != f'{orbital}_up'
and item != 'up']
orbital_p_down = [item for item in DATA.columns.values if
item.startswith(orbital) and item.endswith(
'down') and item != f'{orbital}_down'
and item != 'down']
for item in orbital_p_up:
DATA[f'{orbital}_up'] += DATA[item]
for item in orbital_p_down:
DATA[f'{orbital}_down'] += DATA[item]
DATA['up'] += DATA[f'{orbital}_up']
DATA['down'] += DATA[f'{orbital}_down']
# LORBIT = 10
if self.LORBIT == 10:
for item in columns:
if item.endswith("up") and not item.startswith("s"):
DATA['up'] += DATA[item]
elif item.endswith("down") and not item.startswith("s"):
DATA['down'] += DATA[item]
# for both LORBIT=10/12
DATA['up'] += DATA['s_up']
DATA['down'] += DATA['s_down']
atom_data.append(DATA)
self.TDOS = Total_Dos # DataFrame(NDOS, 2) / DataFrame(NDOS, 1)
self.LDOS = atom_data # energy + List(NAtom, NDOS, NOrbital+8)
return self
return merge_dos(*dos_cython.load(self.strings, self.NDOS, self.ISPIN, self.fermi))
[docs]
class EIGENVAL(MetaFile):
def __init__(self, name):
super(EIGENVAL, self).__init__(name=name)
self.NKPoint, self.NBand = tuple(map(int, self.strings[5].split()[1:]))
self.KPoint_coord = None
self.KPoint_dist = None
self.KPoint_label = None
self.energy = None
self._parse()
def _parse(self):
"""
load Eigenval obtain the band-energy
@return:
self.energy: shape=(NKPoint, NBand, 2)
"""
self.KPoint_coord = []
self.energy = []
energy_flag = False
sBand = []
for index, line in enumerate(self.strings):
if index % (self.NBand + 2) == 7:
self.KPoint_coord.append(tuple(map(float, line.split()[0:3])))
energy_flag = True
continue
elif energy_flag:
sBand.append(list(map(float, line.split()[1:3])))
if len(sBand) == self.NBand:
energy_flag = False
self.energy.append(sBand)
sBand = []
self.KPoint_coord = np.array(self.KPoint_coord)
self.KPoint_dist = []
self.KPoint_label = []
for index, k_point in enumerate(self.KPoint_coord):
if index == 0:
self.KPoint_dist.append(0.)
else:
k_before = self.KPoint_coord[index - 1]
dist = np.sum((k_point - k_before) ** 2) ** 0.5
self.KPoint_dist.append(self.KPoint_dist[-1] + dist)
for label, value in HIGH_SYM.items():
if np.sum((np.array(value) - k_point) ** 2) ** 0.5 <= 1E-02 and (
not len(self.KPoint_label) or self.KPoint_label[-1] != label):
self.KPoint_label.append(label)
break
else:
self.KPoint_label.append("")
self.energy = np.array(self.energy)
return self
[docs]
def write(self, directory='band_data'):
"""
Write band-data to file, each band corresponding to one file and named as band_{index}
@params:
dir: save directory, default: $CWD/band_data
"""
Path(directory).mkdir(exist_ok=True)
for index in range(self.NBand):
np.savetxt(f"{directory}/band_{index + 1}", self.energy[:, index])
logger.info(f"Band data has been saved to {directory} directory")
[docs]
class CHGBase(StructInfoFile):
"""
Subclass of StructInfoFile, inherit <structure property>
"""
def __new__(cls, *args, **kwargs):
if cls is CHGBase:
raise TypeError(f"<{cls.__name__} class> may not be instantiated")
return super().__new__(cls)
def __init__(self, name):
super(CHGBase, self).__init__(name=name)
self.NGX, self.NGY, self.NGZ = None, None, None
self.density = None
def __add__(self, other):
if self.__class__.__name__.startswith("AECCAR") and other.__class__.__name__.startswith("AECCAR"):
if self.density is None and other.density is None:
pool = ProcessPool(processes=2)
results = [pool.apply_async(self.load), pool.apply_async(other.load)] # new instance
self, other = [item.get() for item in results]
pool.close()
pool.join()
elif self.density is None:
self.load()
elif other.density is None:
other.load()
if self.structure != other.structure:
raise StructureNotEqualError(f"{self.name}.structure is not equal to {other.name}.structure")
if (self.NGX, self.NGY, self.NGZ) != (other.NGX, other.NGY, other.NGZ):
raise GridNotEqualError(f"{self.name}.NGrid is not equal to {other.name}.NGrid")
density_sum = self.density + other.density
return CHGCAR_sum.from_array("CHGCAR_sum", self.structure, (self.NGX, self.NGY, self.NGZ), density_sum)
else:
raise TypeError(
f"unsupported operand type(s) for +: {self.__class__.__name__} and {other.__class__.__name__}")
[docs]
def load(self):
"""
load Electronic-Density
@return:
self.density: shape=(NGX, NGY, NGZ)
"""
info = file_bind.load(self.name)
self.NGX, self.NGY, self.NGZ = info.NGX, info.NGY, info.NGZ
self.density = info.density
assert self.density.size == self.NGX * self.NGY * self.NGZ, "Load density failure, size is not consistent"
self.density = self.density.reshape((self.NGX, self.NGY, self.NGZ), order="F")
return self
[docs]
def write(self, title=None, factor=1.0):
"""
write CHGCAR_* file from array
@param:
system: specify the structure system
factor: coordination factor
"""
self.structure.write_POSCAR(name=self.__class__.__name__, title=title, factor=factor)
density_fortran = self.density.reshape(-1, order="F")
# solve residue problem
residue = density_fortran.size % 5
if residue:
density_fortran1 = density_fortran[:-residue].reshape(-1, 5)
density_fortran2 = density_fortran[-residue:].reshape(1, -1)
else:
density_fortran1 = density_fortran.reshape(-1, 5)
density_fortran2 = None
with open(self.__class__.__name__, "a+") as f:
f.write(f"{self.NGX:>5}{self.NGY:>5}{self.NGZ:>5}\n")
np.savetxt(f, density_fortran1, fmt="%18.11E")
if density_fortran2 is not None:
np.savetxt(f, density_fortran2, fmt="%18.11E")
[docs]
class AECCAR0(CHGBase):
pass
[docs]
class AECCAR2(CHGBase):
pass
[docs]
class CHGCAR_sum(CHGBase):
def __init__(self, name):
super(CHGCAR_sum, self).__init__(name=name)
self._structure = None
@property
def structure(self):
if self._structure is None:
self._structure = super(CHGCAR_sum, self).structure()
return self._structure
@structure.setter
def structure(self, _structure):
self._structure = _structure
[docs]
@staticmethod
def from_array(name: str, structure, NGrid: tuple[int, int, int], density):
instance = CHGCAR_sum(name=name)
instance.structure = structure
instance.NGX, instance.NGY, instance.NGZ = NGrid
instance.density = density
return instance
[docs]
class CHGCAR_tot(CHGBase):
pass
[docs]
class CHGCAR_mag(CHGBase):
[docs]
@staticmethod
def to_grd(name="vasp.grd", DenCut=-1):
"""
transform CHGCAR_mag to grd file
param:
name: specify the name of grd file
DenCut: density lower than DenCut will be set to zero (default: -1: disable the DenCut option)
"""
file_bind.to_grd(name, DenCut)
[docs]
class CHGCAR_diff(CHGBase):
[docs]
def line_potential(self, direction='z'):
"""
Calculate Charge Density Difference along one direction, default: z-axis
Args:
direction (str): which axis you want to Average the CCD
Returns:
line_potential (np.array[:]): CCD along one axis
"""
mapping = {'x': 0, 'y': 1, 'z': 2}
if getattr(self, 'density', None) is None:
self.load()
NGXYZ = (self.NGX, self.NGY, self.NGZ)
if mapping.get(direction, None) is None:
raise KeyError(f"{direction} is not supported, should be [x, y, z]")
# TODO: python: malloc.c:2617: sysmalloc: Assertion failed for linux system (win success).
# `(old_top == initial_top (av) && old_size == 0) || ((unsigned long) (old_size) >= MINSIZE &&
# prev_inuse (old_top) && ((unsigned long) old_end & (pagesize - 1)) == 0)'
# Guess the `file_bind.so` file error
potential_swap = np.swapaxes(self.density / self.structure.lattice.length[mapping[direction]], -1,
mapping[direction])
return np.linspace(start=0, stop=self.structure.lattice.length[mapping[direction]],
num=NGXYZ[mapping[direction]]), np.mean(potential_swap, axis=(0, 1))
[docs]
class CHGCAR(StructInfoFile):
def __init__(self, name):
super(CHGCAR, self).__init__(name=name)
self.NGX, self.NGY, self.NGZ, self.NGrid = None, None, None, None
self.density_tot, self.density_mag = None, None
self._head = None
self._density_tot_strings, self._density_mag_strings = None, None
[docs]
def load(self):
"""
load Electronic-Density
@return:
self.NGrid: NGX * NGY * NGZ
self.density: shape=(NGX, NGY, NGZ)
self._density_strings: density with strings format
"""
start = len(self.structure.atoms) + 9
self.NGX, self.NGY, self.NGZ = tuple(map(int, self.strings[start].split()))
self._head = self.strings[:start + 1]
index = np.where(np.array(self.strings) == self.strings[start])[0]
assert len(index) == 2, f"Search indicator failure"
self.NGrid = self.NGX * self.NGY * self.NGZ
count = self.NGrid / 5 if self.NGrid % 5 == 0 else self.NGrid / 5 + 1
count = int(count)
self._density_tot_strings = self.strings[index[0] + 1:index[0] + 1 + count]
self._density_mag_strings = self.strings[index[1] + 1:index[1] + 1 + count]
# solve residue problem
self._density_tot_strings_list = np.char.split(self._density_tot_strings).tolist()
residue = 5 - len(self._density_tot_strings_list[-1])
self._density_tot_strings_list[-1] += ['0.'] * residue
self._density_mag_strings_list = np.char.split(self._density_mag_strings).tolist()
self._density_mag_strings_list[-1] += ['0.'] * residue
self.density_tot = np.append([], self._density_tot_strings_list).astype(float)
self.density_tot = self.density_tot[:-residue].reshape((self.NGX, self.NGY, self.NGZ), order="F")
self.density_mag = np.append([], self._density_mag_strings_list).astype(float)
self.density_mag = self.density_mag[:-residue].reshape((self.NGX, self.NGY, self.NGZ), order="F")
[docs]
def split(self):
"""split CHGCAR to CHGCAR_tot && CHGCAR_mag"""
if getattr(self, "_head", None) is None:
self.load()
with open("CHGCAR_tot", "w") as tot, open("CHGCAR_mag", "w") as mag:
tot.writelines(self._head)
tot.writelines(self._density_tot_strings)
mag.writelines(self._head)
mag.writelines(self._density_mag_strings)
[docs]
class ACFFile(MetaFile):
def __init__(self, name):
super(ACFFile, self).__init__(name=name)
self.cart_coord = []
self.charge = []
self._parse()
def _parse(self):
for line in self.strings[2:-4]:
item = line.split()
self.cart_coord.append(list(map(float, item[1:4])))
self.charge.append(float(item[4]))
self.cart_coord = np.array(self.cart_coord)
self.charge = np.array(self.charge)
[docs]
class LOCPOT(StructInfoFile):
def __init__(self, name):
super(LOCPOT, self).__init__(name=name)
self.NGX, self.NGY, self.NGZ, self.NGrid = None, None, None, None
self.potential = None
self.lattice = None
self._head = None
self._potential_strings = None
[docs]
def load(self):
"""
load Electrostatic Potential
Returns:
self.NGrid (int): value = NGX * NGY * NGZ
self.potential (np.array[:, :, :]): record the electrostatic potential
self.lattice (Lattice): <Lattice class> instance
"""
start = len(self.structure.atoms) + 9
self.NGX, self.NGY, self.NGZ = tuple(map(int, self.strings[start].split()))
self._head = self.strings[:start + 1]
self.lattice = Lattice.from_string(self._head[2:5])
index = np.where(np.array(self.strings) == self.strings[start])[0]
self.NGrid = self.NGX * self.NGY * self.NGZ
count = self.NGrid / 5 if self.NGrid % 5 == 0 else self.NGrid / 5 + 1
count = int(count)
self._potential_strings = self.strings[index[0] + 1:index[0] + 1 + count]
# solve residue problem
self._potential_strings = np.char.split(self._potential_strings).tolist()
residue = 5 - len(self._potential_strings[-1])
self._potential_strings[-1] += ['0.'] * residue
self.potential = np.append([], self._potential_strings).astype(float)
if residue:
self.potential = self.potential[:-residue].reshape((self.NGX, self.NGY, self.NGZ), order="F")
else:
self.potential = self.potential.reshape((self.NGX, self.NGY, self.NGZ), order="F")
[docs]
def line_potential(self, direction='z'):
"""
Calculate the electrostatic potential along one direction, default: z-axis
Args:
direction (str): which axis you want to calculate the electrostatic potential
Returns:
line_potential (np.array[:]): electrostatic potential along one axis
"""
mapping = {'x': 0, 'y': 1, 'z': 2}
if getattr(self, 'potential', None) is None:
self.load()
NGXYZ = (self.NGX, self.NGY, self.NGZ)
if mapping.get(direction, None) is None:
raise KeyError(f"{direction} is not supported, should be [x, y, z]")
potential_swap = np.swapaxes(self.potential, -1, mapping[direction])
return np.linspace(start=0, stop=self.lattice.length[mapping[direction]],
num=NGXYZ[mapping[direction]]), np.mean(potential_swap, axis=(0, 1))
[docs]
class OUTCAR(MetaFile):
def __init__(self, name):
super(OUTCAR, self).__init__(name=name)
element_name = [line.split()[3].split("_")[0] for line in self.strings if "TITEL" in line]
element_count = [list(map(int, line.split()[4:])) for line in self.strings if "ions per" in line][0]
self.element = sum([[name] * count for name, count in zip(element_name, element_count)], [])
self.lattice = {Lattice.from_string(self.strings[index + 1:index + 4]) for index, line in
enumerate(self.strings) if "direct lattice vectors" in line}
self.lattice = list(self.lattice)[0] if len(self.lattice) == 1 else self.lattice
self._frequency = [i for i in range(len(self.strings)) if "Hz" in self.strings[i]]
self._neb = [line for line in self.strings if "NEB:" in line]
self._fort = [line for line in self.strings if "fort.1881" in line]
self.spin, self.bands, self.kpoints, self.fermi, self.steps = None, None, None, None, None
self.energy, self.force, self.mag = None, None, None
self.last_energy, self.last_force, self.last_mag = None, None, None
self.finish = "reached" in self.strings[-30]
self._parse_base()
self.frequency = None
if len(self._frequency):
self._parse_freq()
self.kpoint_info, self.band_info = None, None
if self.finish:
self._parse_band()
self.tangent, self.last_tangent = 0., 0.
if len(self._neb):
self._parse_neb()
self.condist, self.last_condist = 0., 0.
if len(self._fort):
self._parse_fort()
def _parse_base(self):
self.spin = [int(line.split()[2]) for line in self.strings if "ISPIN" in line][0]
self.bands, self.kpoints = \
[(int(line.split()[-1]), int(line.split()[3])) for line in self.strings if
"NBANDS" in line and "the" not in line][0]
steps = [(index, int(line.split()[2].split("(")[0]), int(line.split()[3].split(")")[0]))
for index, line in enumerate(self.strings) if "Iteration" in line]
self.steps = namedtuple("Steps", ("index", "ionic", "electronic"))(*list(map(tuple, zip(*steps))))
self.fermi = [float(line.split()[2]) for line in self.strings if "E-fermi" in line][-1]
self.energy = [float(line.split()[3]) for line in self.strings if "energy without" in line]
self.force = [float(line.split()[5]) for line in self.strings if "FORCES: max atom" in line]
if self.spin == 2:
self.mag = [float(line.split()[5]) for line in self.strings if "number of electron " in line]
else:
self.mag = []
self.last_energy = self.energy[-1] if len(self.energy) else None
self.last_force = self.force[-1] if len(self.force) else None
self.last_mag = self.mag[-1] if len(self.mag) else None
def _parse_freq(self):
"""
Parse frequency information from OUTCAR
@return:
register self.frequency attr (type: namedtuple)
"""
image, wave_number, vib_energy, coord, vibration = [], [], [], [], []
for index in self._frequency:
item = self.strings[index].split()
image.append(False) if item[1] == "f" else image.append(True)
wave_number.append(float(item[-4]))
vib_energy.append(float(item[-2]))
item = list(map(lambda x: [float(i) for i in x],
np.char.split(self.strings[index + 2:index + 2 + len(self.element)])))
coord.append(np.array(item)[:, :3])
vibration.append(np.array(item)[:, 3:])
self.frequency = namedtuple("Frequency",
("image", "wave_number", "vib_energy", "coord", "vibration"))(image,
np.array(wave_number),
np.array(vib_energy),
np.array(coord),
np.array(vibration))
def _parse_band(self):
"""
Parse band information from OUTCAR
@return:
register self.kpoint_info (type: tuple(List[KPoint(coord, value)], List[KPoint(coord, value)]))
"""
def spin_obtain(kpoint_group):
"""
According to the index-list of kpoint in OUTCAR, parse the band_info
@param:
kpoint_group: list, record the index of 'k-point' field occurrence in OUTCAR
@return:
spin: List[KPoint(coord, value)], record the each k-point info
"""
spin = []
trans_func = lambda x: [int(x[0]), float(x[1]) - self.fermi, float(x[2])] # format transform
# parse band_info
for index in kpoint_group:
coord = list(map(float, band_info[index].split()[3:]))
value = list(map(trans_func, np.char.split(band_info[index + 2:index + 2 + self.bands])))
spin.append(KPoint(coord, np.array(value)))
return spin
def transform_band(spin):
band = []
for kpoint in spin:
band.append(kpoint.value[:, 1])
band = np.array(band)
band = band.transpose((1, 0))
return band
content = self.strings[self.steps.index[-1]:] # calculate bandgap from last step
start_index = [index for index, line in enumerate(content) if "E-fermi" in line][0]
end_index = [index for index, line in enumerate(content) if "-----" in line and index > start_index][0]
band_info = content[start_index:end_index] # band_info content in last step
KPoint = namedtuple("KPoint", ("coord", "value")) # including `kpoint coord` and `each band energy`
if self.spin == 2:
kpoint_index = [index for index, line in enumerate(band_info) if "k-point" in line]
half_len = int(len(kpoint_index) / 2)
kpoint_up, kpoint_down = kpoint_index[:half_len], kpoint_index[half_len:] # 'k-point' occurrence index
spin_up = spin_obtain(kpoint_up)
spin_down = spin_obtain(kpoint_down)
self.kpoint_info = namedtuple("KPoint_info", ("up", "down"))(spin_up, spin_down)
self.band_info = namedtuple("Band_info", ("up", "down"))(transform_band(spin_up), transform_band(spin_down))
else:
logger.warning("Non-spin polarized calculation may have error, please check!!")
def _parse_neb(self):
self.tangent = [float(line.split()[-1]) for line in self.strings if line.find("tangent") != -1]
self.last_tangent = self.tangent[-1]
def _parse_fort(self):
self.condist = [float(line.split()[-1]) for line in self.strings if line.find("distance after opt") != -1]
self.last_condist = self.condist[-1]
[docs]
def bandgap(self, cutoff=0.01):
"""
Calculated the bandgap from OUTCAR file
@param:
cutoff: any occupy lower than cutoff will be treated as the empty state
@return:
type: type of bandgap, ['direct', 'indirect']
bandgap: value of bandgap
"""
def MO_obtain(spin):
"""
According to the kpoint_info, return the homo and lumo
@param:
spin: List[KPoint(coord, value)], record the each k-point info
@return:
homo_spin: KPoint(coord, value), record the homo information
lumo_spin: KPoint(coord, value), record the lumo information
"""
# calculate homo and lumo
homo_spin, lumo_spin = None, None
for item in spin:
lumo_index = np.where(item.value[:, 2] <= cutoff)[0][0]
homo_index = lumo_index - 1
homo, lumo = item.value[homo_index], item.value[lumo_index]
homo_spin = KPoint(item.coord, homo) if homo_spin is None or homo_spin.value[1] < homo[1] else homo_spin
lumo_spin = KPoint(item.coord, lumo) if lumo_spin is None or lumo_spin.value[1] > lumo[1] else lumo_spin
return homo_spin, lumo_spin
KPoint = namedtuple("KPoint", ("coord", "value")) # including `kpoint coord` and `each band energy`
if self.spin == 2:
homo_up, lumo_up = MO_obtain(self.kpoint_info[0])
homo_down, lumo_down = MO_obtain(self.kpoint_info[1])
homo_real = homo_up if homo_up.value[1] > homo_down.value[1] else homo_down
lumo_real = lumo_up if lumo_up.value[1] < lumo_down.value[1] else lumo_down
if homo_real.coord == lumo_real.coord:
return "direct", lumo_real.value[1] - homo_real.value[1]
else:
return "indirect", lumo_real.value[1] - homo_real.value[1]
else:
raise NotImplementedError("Non-spin polarized calculation has not been implemented")
[docs]
def animation_freq(self, freq: [str, int] = "image", frames: int = 30, scale: float = 0.6):
"""
Generate freq.arc file from OUTCAR
@param:
freq: specify which freq you want to animate, accept [int, str] arguments
frames: specify how many points you want to interpolate along one direction, default = 30
scale: determine the vibration scale, default = 0.6
"""
if self.frequency is None:
raise AnimationError(f"{self.name} don't include frequency information")
if isinstance(freq, int) and freq not in range(len(self._frequency)):
raise FrequencyError(f"freq{freq} is not in {self.name}, should be {list(range(len(self._frequency)))}")
if isinstance(freq, str) and freq != "image":
raise FrequencyError(f"`{freq}` is not supported, should be `image`")
freq = [index for index, item in enumerate(self.frequency.image) if item] if isinstance(freq, str) else [freq]
if len(freq) == 0:
raise FrequencyError(f"`{freq}` frequency is not found")
# generate the directions for vibration
direction_001 = np.linspace(start=0, stop=scale, num=frames) # 0-1 direction
direction_010 = np.linspace(start=scale, stop=0, num=frames) # 1-0 direction
direction_101 = np.linspace(start=0, stop=-scale, num=frames) # 0-^1 direction
direction_110 = np.linspace(start=-scale, stop=0, num=frames) # ^1-0 direction
direction_all = np.concatenate([direction_001, direction_010, direction_101, direction_110])
direction_all = direction_all[:, np.newaxis, np.newaxis]
# generate multi-frames new coordinates
coord_freq = []
for freq_index in freq:
logger.info(f"Processing freq{freq_index + 1} file ...")
coord = np.repeat(self.frequency.coord[freq_index][np.newaxis, :], direction_all.shape[0], axis=0)
vibration = np.repeat(self.frequency.vibration[freq_index][np.newaxis, :], direction_all.shape[0], axis=0)
coord_freq.append(coord + vibration * direction_all)
# generate the *.arc file
if not isinstance(self.lattice, Lattice):
raise NotImplementedError("we here only considered the lattice unchangeable for the whole calculation")
for freq_index, coord_index in zip(freq, coord_freq):
structure = [Structure(atoms=Atoms(formula=self.element, cart_coord=cart_coord), lattice=self.lattice)
for cart_coord in coord_index]
ARCFile.write(name=f"freq{freq_index + 1}.arc", structure=structure, lattice=self.lattice)
logger.info(f"All freq transform to corresponding *.arc files")
[docs]
class MODECAR(MetaFile):
[docs]
@staticmethod
def write_from_freq(freq: int, scale: float, outcar="OUTCAR"):
"""
Generate MODECAR file from the `image-freq`
@param:
freq: which freq you want to generate MODECAR
scale: scale factor, may lower than 1.0
outcar: name of the OUTCAR
"""
outcar = OUTCAR(name=outcar)
frequency = outcar.frequency.vibration[freq]
np.savetxt("MODECAR", frequency * scale, fmt="%10.5f")
[docs]
@staticmethod
def write_from_POSCAR(pos1: str, pos2: str):
"""
Generate MODECAR file from the two `POSCAR` file
@param:
pos1: name of first POSCAR file
pos2: name of second POSCAR file
"""
structure1 = POSCAR(pos1).structure
structure2 = POSCAR(pos2).structure
assert structure1 == structure2, f"{pos1} and {pos2} are not structure match"
diff = structure1 - structure2
dist = Structure.dist(structure1, structure2)
diff_norm = diff / dist
np.savetxt("MODECAR", diff_norm, fmt="%20.10E")