diff --git a/doc/DoxygenLayout.xml b/doc/DoxygenLayout.xml index fbe3c41e1..ac353bf86 100644 --- a/doc/DoxygenLayout.xml +++ b/doc/DoxygenLayout.xml @@ -1,5 +1,5 @@ - + @@ -24,11 +24,10 @@ - + - @@ -63,6 +62,7 @@ + @@ -83,8 +83,7 @@ - - + @@ -95,6 +94,7 @@ + @@ -107,12 +107,11 @@ - + - @@ -124,6 +123,7 @@ + @@ -137,9 +137,8 @@ - + - @@ -161,6 +160,7 @@ + @@ -183,12 +183,12 @@ - + - + diff --git a/python/casm/README.md b/python/casm/README.md index 53f44668c..4c6219514 100644 --- a/python/casm/README.md +++ b/python/casm/README.md @@ -112,6 +112,9 @@ Individual Package dependencies include: - **prisms_jobs** ([https://prisms-center.github.io/prisms_jobs_docs](https://prisms-center.github.io/prisms_jobs_docs)) +- **atomate** ([https://github.com/hackingmaterials/atomate])(https://github.com/hackingmaterials/atomate) + +- ** pymatgen** ([https://github.com/materialsproject/pymatgen])(https://github.com/materialsproject/pymatgen) Generating html documentation ----------------------------- diff --git a/python/casm/casm/aims/__init__.py b/python/casm/casm/aims/__init__.py new file mode 100644 index 000000000..0a6dc6c21 --- /dev/null +++ b/python/casm/casm/aims/__init__.py @@ -0,0 +1,3 @@ +"""A module for interacting with FHI-aims""" + +__all__ = dir() diff --git a/python/casm/casm/aims/aims.py b/python/casm/casm/aims/aims.py new file mode 100644 index 000000000..25b7d007a --- /dev/null +++ b/python/casm/casm/aims/aims.py @@ -0,0 +1,294 @@ +import os +import shutil +import re +import subprocess +import sys +import time +import gzip +import warnings + +import casm.project +import casm.aimswrapper +import casm.aims +import casm.aims.io.geometry + +from casm.aims.io.geometry import Geometry +from casm.aims.io.basis import basis_settings +from casm.aims.io.io import DEFAULT_AIMS_GZIP_LIST, DEFAULT_AIMS_COPY_LIST, \ + DEFAULT_AIMS_REMOVE_LIST, DEFAULT_AIMS_MOVE_LIST + + +class AimsError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class AimsWarning(Warning): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +def continue_job(jobdir, contdir, settings): + """Use the files in job directory 'jobdir', to setup a job in directory 'contdir'. + + Args: + jobdir: path to current job directory + contdir: path to new directory to continue the job + settings: + copy: Files copied from 'jobdir' to 'contdir' + It also copies either geometry.in.next_step or if that does not exist geometry.in + move: Files moved from 'jobdir' to 'contdir' + keep: Files (along with those in 'copy') to keep in 'jobdir'. The rest are removed. + """ + + configdir = os.getcwd() + _res = os.path.split(configdir) + cfgname = os.path.split(_res[0])[1] + "/" + _res[1] + casm_dirs = casm.project.DirectoryStructure(configdir) + casm_sets = casm.project.ProjectSettings(configdir) + clex = casm_sets.default_clex + aimsfiles = casm.aimswrapper.aims_input_file_names(casm_dirs, cfgname, clex) + controlfile, prim_posfile, super_posfile, basisfile = aimsfiles + + print("Continue FHI-aims job:\n Original: " + jobdir + "\n Continuation: " + contdir) + sys.stdout.flush() + + # remove duplicates + move = list(set(settings['move'] + DEFAULT_AIMS_MOVE_LIST)) + copy = list(set(settings['copy'] + DEFAULT_AIMS_COPY_LIST)) + remove = list(set(settings['remove'] + DEFAULT_AIMS_REMOVE_LIST)) + compress = list(set(settings['compress'])) + + # Check that the user isn't being contradictory or silly + for f in move: + if f == "geometry.in" or f == "geometry.in.next_step": + raise AimsError("Error in casm.vasp.general.continue_job(). " + "Do not include geo.in and geo.in.next_step 'move'; use 'backup' if you want a backup") + if f in remove: + if f in DEFAULT_AIMS_MOVE_LIST: + raise AimsError("Error in casm.aims.general.continue_job(). " + "%s cannot be removed, FHI-aims will not run!!!" % f) + else: + warnings.warn("Warning: %s found in both 'move' and 'remove'. " + "The file will not be removed." % f, AimsWarning) + remove = list(set(remove) - set(f)) + if f in copy: + if f in DEFAULT_AIMS_MOVE_LIST: + warnings.warn("Warning: %s found in both 'move' and 'copy'. " + "The fill will be moved only." % f, AimsWarning) + copy = list(set(copy) - set(f)) + else: + warnings.warn("Warning: %s found in both 'move' and 'copy'. " + "The file will be copied only." % f, AimsWarning) + move = list(set(move) - set(f)) + if f in compress: + if f in DEFAULT_AIMS_GZIP_LIST: + raise AimsError("Error in casm.aims.general.continue_job(). " + "%s cannot be compressed, FHI-aims will not run!!!" % f) + else: + raise AimsError("Error in casm.aims.general.continue_job(). " + "%s found in both 'move' and 'compress', but these options contradict. " + "Did you mean 'backup'?" % f) + + # make the new contdir + if not os.path.isdir(contdir): + os.mkdir(contdir) + + # copy/move files + if os.path.isfile(os.path.join(jobdir, "geometry.in.next_step")) and \ + os.path.getsize(os.path.join(jobdir, "geometry.in.next_step")) > 0: + shutil.move(os.path.join(jobdir, "control.in"), os.path.join(contdir, "control.in")) + my_basis = basis_settings(basisfile) + newgeo = Geometry(os.path.join(jobdir, "geometry.in.next_step")) + newgeo.write(os.path.join(contdir, "geometry.in"), my_basis) + print(" cp geometry.in.next_step -> geometry.in (and added default moments)") + else: + shutil.copyfile(os.path.join(jobdir, "geometry.in"), os.path.join(contdir, "geometry.in")) + print(" no geometry.in.next_step, took old geometry.in (this means no relaxation step was performed)") + shutil.move(os.path.join(jobdir, "control.in"), os.path.join(contdir, "control.in")) + + +def complete_job(jobdir, settings): + """Clean up and compress output + + Args: + jobdir: path to current job directory + settings: casm settings + """ + print("Completing FHI-aims job: " + jobdir) + sys.stdout.flush() + + # compress files + print(" gzip: ", end='') + for file in settings["compress"]: + if os.path.isfile(os.path.join(jobdir, file)): + print(file, end='') + # Open target file, target file.gz + f_in = open(os.path.join(jobdir, file), 'rb') + f_out = gzip.open(os.path.join(jobdir, file)+'.gz', 'wb') + # Compress, close files + f_out.writelines(f_in) + f_out.close() + f_in.close() + # Remove original target file + os.remove(os.path.join(jobdir, file)) + print(" gzipping DONE") + print("\n") + sys.stdout.flush() + + +class FreezeError(object): + """VASP appears frozen""" + def __init__(self): + self.pattern = None + + def __str__(self): + return "FHI-aims appears to be frozen" + + @staticmethod + def error(jobdir=None): + """ Check if aims is frozen + + Args: + jobdir: job directory + Returns: + True if: + 1) no file has been modified for 5 minutes + 2) 'LOOP+' exists in OUTCAR and no output file has been modified + in 5x the time for the slowest loop + """ + + # Check if any files modified in last 300 s + most_recent = None + most_recent_file = None + for f in os.listdir(jobdir): + t = time.time() - os.path.getmtime(os.path.join(jobdir, f)) + if most_recent is None: + most_recent = t + most_recent_file = f + elif t < most_recent: + most_recent = t + most_recent_file = f + + print("Most recent file output (" + most_recent_file + "):", most_recent, " seconds ago.") + sys.stdout.flush() + if most_recent < 300: + return False + + @staticmethod + def fix(err_jobdir, new_jobdir, settings): + """ Fix by killing the job and resubmitting.""" + print(" Kill job and continue...") + continue_job(err_jobdir, new_jobdir, settings) + + +def error_check(jobdir, stdoutfile): + """ Check vasp stdout for errors """ + err = dict() + + # Error to check line by line, only look for first of each type + sout = open(stdoutfile, 'r') + + # Error to check for once + possible = [FreezeError()] + for p in possible: + if p.error(jobdir=jobdir): + err[p.__class__.__name__] = p + + sout.close() + if len(err) == 0: + return None + else: + return err + + +def run(jobdir=None, stdout="std.out", stderr="std.err", command=None, ncpus=None, + poll_check_time=5.0, err_check_time=60.0): + """ Run FHI-aims using subprocess. + + The 'command' is executed in the directory 'jobdir'. + + Args: + jobdir: directory to run. If jobdir is None, the current directory is used. + stdout: filename to write to. If stdout is None, "std.out" is used. + stderr: filename to write to. If stderr is None, "std.err" is used. + command: (str or None) FHI-aims execution command + If command != None: then 'command' is run in a subprocess + Else, if ncpus == 1, then command = "aims" + Else, command = "mpirun -np {NCPUS} aims" + ncpus: (int) if '{NCPUS}' is in 'command' string, then 'ncpus' is substituted in the command. + if ncpus==None, $PBS_NP is used if it exists, else 1 + poll_check_time: how frequently to check if the vasp job is completed + err_check_time: how frequently to parse vasp output to check for errors + """ + print("Begin FHI-aims run:") + sys.stdout.flush() + + if jobdir is None: + jobdir = os.getcwd() + + currdir = os.getcwd() + os.chdir(jobdir) + + if ncpus is None: + if "PBS_NP" in os.environ: + ncpus = os.environ["PBS_NP"] + else: + ncpus = 1 + + if command is None: + if ncpus == 1: + command = "aims" + else: + command = "mpirun -np {NCPUS} aims" + + if re.search("NCPUS", command): + command = command.format(NCPUS=str(ncpus)) + + print(" jobdir:", jobdir) + print(" exec:", command) + sys.stdout.flush() + + sout = open(os.path.join(jobdir, stdout), 'w') + serr = open(os.path.join(jobdir, stderr), 'w') + err = None + p = subprocess.Popen(command.split(), stdout=sout, stderr=serr) + + # wait for process to end, and periodically check for errors + last_check = time.time() + while p.poll() is None: + time.sleep(poll_check_time) + if time.time() - last_check > err_check_time: + last_check = time.time() + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + # FreezeErrors are fatal and usually not helped with abort_scf + if "FreezeError" in err.keys(): + print(" FHI-aims seems frozen, killing job") + sys.stdout.flush() + p.kill() + + # close output files + sout.close() + serr.close() + + os.chdir(currdir) + + print("Run complete") + sys.stdout.flush() + + # check finished job for errors + if err is None: + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + print(" Found errors:", end='') + for e in err: + print(e, end='') + print("\n") + return err diff --git a/python/casm/casm/aims/io/__init__.py b/python/casm/casm/aims/io/__init__.py new file mode 100755 index 000000000..da4e0102a --- /dev/null +++ b/python/casm/casm/aims/io/__init__.py @@ -0,0 +1,3 @@ +"""Tools for FHI-aims input and output""" + +__all__ = dir() diff --git a/python/casm/casm/aims/io/aimsrun.py b/python/casm/casm/aims/io/aimsrun.py new file mode 100644 index 000000000..12f53ff8e --- /dev/null +++ b/python/casm/casm/aims/io/aimsrun.py @@ -0,0 +1,133 @@ +import os +import gzip + + +class AimsRunError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class AimsRun: + """ An object + + The AimsRun class contains: + self.total_energy: final total energy + self.forces: final forces on atoms (2d list of double) + self.atom_type: list of atom types (list of str) + self.atoms_per_type: list of number of atoms for each atom type (list of int) + self.lattice: final lattice (2d list of double) + self.rec_lat: final reciprocal lattice (2d list of double) + self.basis: final basis (2d list of double) + self.coord_mode: coordinate mode for basis (always 'direct') + self.is_complete: did the calculation run to completion? (bool) + self.eigenvalues: eigenvalues and energies for doing band structure plots (list of 2d list of double) + self.all_e_0: energy (e_0_energy) for each electronic step of each ionic step + self.nelm: NELM (max. number of electronic steps) + """ + def __init__(self, filename): + """ Create a run object """ + self.filename = filename + self.lattice = [] + self.total_energy = None + self.forces = [] + self.atom_type = [] + self.atoms_per_type = [] + self.lattice = [] + self.rec_lat = [] + self.basis = [] + self.coord_mode = None + self.is_complete = False + self.efermi = None + self.all_e_0 = [] + + # Parse geometry first + geofile = self.filename.replace("std.out", "geometry.in.next_step") + + if not os.path.isfile(geofile): + err_str = 'Initial geometry is already relaxed, this is not what you want\n' + err_str += 'Quitting, please perturb input geometry to force relaxation.' + raise AimsRunError(err_str) + + with open(geofile, 'rt') as f: + atom = [] + name = [] + for line in f: + line = line.split() + if len(line) > 1: + if line[0] == 'lattice_vector': + self.lattice.append([float(x) for x in line[1:4]]) + if line[0] == 'atom_frac': + self.coord_mode = 'direct' + atom.append([float(x) for x in line[1:4]]) + name.append(line[4]) + if line[0] == 'atom': + self.coord_mode = 'cartesian' + atom.append([float(x) for x in line[1:4]]) + name.append(line[4]) + + symbols = [name[0]] + for i in range(len(name)): + tmp = [] + for j in range(len(symbols)): + tmp.append(symbols[j]) + if tmp.count(name[i]) == 0: + symbols.append(name[i]) + + nums = {} + for i in range(len(symbols)): + k = 0 + for j in range(len(atom)): + if symbols[i] == name[j]: + k += 1 + nums[symbols[i]] = k + + for i in range(len(symbols)): + self.atom_type.append(symbols[i]) + self.atoms_per_type.append(nums[symbols[i]]) + for j in range(len(name)): + if name[j] == symbols[i]: + self.basis.append(atom[j]) + + # Parse output now + if os.path.isfile(self.filename): + if self.filename.split('.')[-1].lower() == 'gz': + f = gzip.open(self.filename, 'rt') + else: + f = open(self.filename) + elif os.path.isfile(self.filename + '.gz'): + f = gzip.open(self.filename + '.gz', 'rt') + else: + raise AimsRunError('file not found: ' + self.filename) + + print('read file', self.filename) + + k = 0 + read_forces = False + for line in f: + if 'Total atomic forces' in line: + print(line) + read_forces = True + self.forces = [] + k = 0 + if read_forces: + line = line.split() + if k == 0: # skip line with text trigger + k += 1 + continue + self.forces.append([float(x) for x in line[2:5]]) + k += 1 + if k >= len(atom): + read_forces = False + + if '| Total energy :' in line: + self.all_e_0.append(float(line.split()[6])) + + if 'Total energy of the DFT / Hartree-Fock s.c.f. calculation :' in line: + self.total_energy = float(line.split()[11]) + + def is_complete(self): + """ Return True if FHI-aims ran to completion """ + return self.is_complete diff --git a/python/casm/casm/aims/io/basis.py b/python/casm/casm/aims/io/basis.py new file mode 100644 index 000000000..d2d0838c6 --- /dev/null +++ b/python/casm/casm/aims/io/basis.py @@ -0,0 +1,119 @@ +import os +import re +import glob +import casm +import casm.project +import casm.aimswrapper + +""" +Sample basis File: +BASIS_DIR_PATH = /absolute/path/to/species_defaults +SPECIES ALIAS init_mom +Mn3 Mn 3 +Mn4 Mn 4 +""" + + +class BasisError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class IndividualBasis: + """ + The IndividualBasis class contains: + self.name: the name as listed in the POS file + self.alias: the species file lists the name, for convenience only + self.tags: (dict) the tags that need to be modified in the geometry.in for this specie + (i.e. adding initial_moment) + All values are stored as strings. + self.basisdir_base: common directory for all basis inputs + self.basis_location: location of basis directory relative to self.basisdir_base + self.basisdir: directory containing particular basis (self.basisdir_base joined with self.basis_location) + self.basis_file: the filename of the basis set for the specie + """ + + def __init__(self, values, tags, basisdir_base): + """ Construct an IndividualBasis. + + Args: + values: (str list) entries in basis file row + tags: (str list) column names 4+ in basis file, basis tags that need to be modified + basisdir_base: (str) common directory for all basis inputs + """ + + configdir = os.getcwd() + _res = os.path.split(configdir) + cfgname = os.path.split(_res[0])[1] + "/" + _res[1] + casm_dirs = casm.project.DirectoryStructure(configdir) + casm_sets = casm.project.ProjectSettings(configdir) + clex = casm_sets.default_clex + setfile = casm_dirs.settings_path_crawl("relax.json", cfgname, clex) + settings = casm.aimswrapper.read_settings(setfile) + + if len(values) != (len(tags)): + raise BasisError("Length of values != length of tags.\nvalues = " + str(values) + "\ntags = " + str(tags)) + self.name = values[0] + self.alias = values[1] + self.init_mom = values[2] + + try: + self.write_basis = not (float(values[2]) == 0) + except ValueError: + raise BasisError("Could not read basis: " + str(values)) + self.basisdir_base = os.path.join(basisdir_base, settings["basis"]) + + file_pattern = str(self.basisdir_base)+"/*_"+str(values[1]) + "_*" + for basis_file in glob.glob(file_pattern): + self.filename = basis_file + + self.basisdir = self.basisdir_base + self.tags = dict() + for i, key in enumerate(tags): + self.tags[key] = values[i] + + +def basis_settings(filename): + """ Returns a dict of IndividualBasis objects, with keys equal to their names. """ + try: + file = open(filename) + except IOError: + raise BasisError("Could not open: '" + filename + "'") + + # Read BASIS_DIR_PATH from first line + line = file.readline() + m = re.match("BASIS_DIR_PATH\s*=\s*(.*)", line) + if not m: + err_str = 'Could not read BASIS_DIR_PATH.\n' + err_str += 'Expected: BASIS_DIR_PATH = /path/to/location/of/species_defaults/\n' + err_str += 'Found: "' + line + '"' + raise BasisError(err_str) + basis_dir_path = m.group(1) + + # Parsing the header + header = file.readline().strip() + column_names = header.split() + if len(column_names) < 3: + raise BasisError("Insufficient number of columns in basis file") + tags = column_names[:3] + basis_settings_loc = dict() + for line in file: + if line.strip(): + values = line.strip().split() + basis_settings_loc[values[1]] = IndividualBasis(values, tags, basis_dir_path) + file.close() + + return basis_settings_loc + + +def write_basis(filename, basis): + file = open(filename, 'a') + for s in basis: + with open(basis[s].filename, 'r') as bf: + species_data = bf.readlines() + for line in species_data: + file.write(line) + file.close() diff --git a/python/casm/casm/aims/io/control.py b/python/casm/casm/aims/io/control.py new file mode 100644 index 000000000..d780909f9 --- /dev/null +++ b/python/casm/casm/aims/io/control.py @@ -0,0 +1,39 @@ +class ControlError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Control: + """ + The Control class contains: + tags: a dict of all settings + """ + + def __init__(self, filename): + """ Construct a Control object from 'filename'""" + + self.ctrl_content = [] + + try: + f = open(filename, 'r') + except IOError: + raise ControlError("Could not open file: '" + filename + "'") + + # read control.in as is + for line in f: + self.ctrl_content.append(line.strip()) + f.close() + + def write(self, filename): + try: + outfile = open(filename, 'w') + except IOError as e: + raise e + for line in self.ctrl_content: + line += '\n' + outfile.write(line) + + outfile.close() diff --git a/python/casm/casm/aims/io/geometry.py b/python/casm/casm/aims/io/geometry.py new file mode 100755 index 000000000..c7479101c --- /dev/null +++ b/python/casm/casm/aims/io/geometry.py @@ -0,0 +1,302 @@ +import os +import numpy as np + +from casm.vasp.io.poscar import Poscar, Site +from casm.project import DirectoryStructure, ProjectSettings +import casm.aimswrapper + + +class GeometryError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Geometry: + """ + The Geometry class contains: + self._lattice: the lattice; lattice vectors are stored as rows in a numpy array + self._reciprocal_lattice: the reciprocal lattice; lattice vectors are stored as rows in a numpy array + self.sel_dyn: True or False, Selective Dynamics Flag + self.type_atoms: lists the specie names as in the POS (ex. [Mn3 Mn4]) + self.type_atoms_alias: lists the POTCAR names for the species + self.num_atoms: lists the atoms as in the POS (ex. [1 1]) + self.coord_mode: Contains the coordinate mode text from POSCAR, with whitespace stripped from beginning and end + self.basis: a list of Site objects + + """ + def __init__(self, filename): + """ Construct a Geometry object from 'filename' + + Args: + filename = file to read + """ + self.basis = [] + self._lattice = np.zeros((3, 3)) + self._reciprocal_lattice = np.zeros((3, 3)) + self.coord_mode = '' + self.sel_dyn = False + self.SD_FLAG = False + self.num_atoms = [] + self.type_atoms = [] + self.type_atoms_alias = list(self.type_atoms) + + try: + file = open(filename, 'rb') + except IOError: + raise GeometryError("Could not read file: " + filename) + + if filename.split('/')[-1] == 'POS': + self.read_pos(filename) + else: + self._read_geo(file) + + file.close() + + @staticmethod + def read_casm_settings(): + configdir = os.getcwd() + _res = os.path.split(configdir) + cfgname = os.path.split(_res[0])[1] + "/" + _res[1] + casm_dirs = DirectoryStructure(configdir) + casm_sets = ProjectSettings(configdir) + clex = casm_sets.default_clex + setfile = casm_dirs.settings_path_crawl("relax.json", cfgname, clex) + settings = casm.aimswrapper.read_settings(setfile) + return settings + + def write(self, filename, species_data): + """ Write geometry to filename """ + try: + file = open(filename, 'w') + except IOError: + raise GeometryError("Could not write: " + filename) + + settings = self.read_casm_settings() + + file.write('#auto-generated geometry.in by CASM\n') + for i in range(3): + file.write("lattice_vector %.8f %.8f %.8f\n" % + (self._lattice[i, 0], self._lattice[i, 1], self._lattice[i, 2])) + + if self.coord_mode[0].lower() == 'd': + for s in self.basis: + file.write('atom_frac %8.8f %8.8f %8.8f %s\n' % + (s.position[0], s.position[1], s.position[2], s.occupant)) + file.write(' initial_moment %3.3f\n' % float(species_data[s.occupant].init_mom)) + if settings['is_slab']: + if s.position[2] < float(settings['fix_pos']) / self._lattice[2, 2]: + file.write(' constrain_relaxation .true.\n') + else: + for s in self.basis: + file.write('atom_frac %8.8f %8.8f %8.8f %s\n' % + (s.position[0], s.position[1], s.position[2], s.occupant)) + file.write(' initial_moment %3.3f\n' % float(species_data[s.occupant].init_mom)) + if settings['is_slab']: + if s.position[2] < float(settings['fix_pos']): + file.write(' constrain_relaxation .true.\n') + + def lattice(self, index=None): + """ Returns the lattice, or lattice vector 'index', as numpy array""" + if index is not None: + return self._lattice[index, :] + return self._lattice + + def reciprocal_lattice(self, index=None): + """ Returns the reciprocal lattice vector 'index', as numpy array""" + if index is not None: + return self._reciprocal_lattice[index, :] + return self._reciprocal_lattice + + def volume(self, lattice=None): + """ Returns scalar triple product of lattice vectors """ + if lattice is None: + lattice = self._lattice + return np.inner(lattice[0, :], np.cross(lattice[1, :], lattice[2, :])) + + def reciprocal_volume(self, reciprocal_lattice=None): + """ Returns scalar triple product of reciprocal lattice vector """ + if reciprocal_lattice is None: + reciprocal_lattice = self._reciprocal_lattice + return self.volume(reciprocal_lattice) + + def basis_dict(self): + """ Return a dictionary where keys are unique specie 'alias' and values are lists of atoms.""" + basis_dict = dict() + for i, atom in enumerate(self.type_atoms_alias): + start = sum(self.num_atoms[0:i]) + end = start + self.num_atoms[i] + if atom not in basis_dict.keys(): + basis_dict[atom] = self.basis[start:end] + else: + basis_dict[atom] += self.basis[start:end] + return basis_dict + + def unsort_dict(self): + """ Return a dict to unsort atom order. + + Returns 'unsort_dict', for which: unsorted_dict[orig_index] == sorted_index; + + unsorted_dict[sorted_index] == orig_index + + For example: + 'unsort_dict[0]' returns the index into the unsorted POSCAR of the first atom in the sorted POSCAR + """ + # create basis_dict, but with initial position in POSCAR instead of coordinate + basis_dict = dict() + for i, atom in enumerate(self.type_atoms_alias): + start = sum(self.num_atoms[0:i]) + end = start + self.num_atoms[i] + if atom not in basis_dict.keys(): + basis_dict[atom] = range(start, end) + else: + basis_dict[atom] += range(start, end) + + orig_pos = [] + for atom in sorted(basis_dict.keys()): + orig_pos += basis_dict[atom] + + new_pos = range(0, len(self.basis)) + + return dict(zip(new_pos, orig_pos)) + + def _read_geo(self, file): + """ Called by __init__ to read the lattice into self._lattice + Args: + file: an open geometry.skel being read from + + self._lattice contains lattice vectors stored as rows in a numpy array (easy inversion) + """ + atom_read = False + atom_name = None + cart = None + lat = [] + cont = file.readlines() + for line in cont: + pos = np.empty(3) + if b'lattice_vector' in line: + lat.append([float(x) for x in line.split()[1:4]]) + + if b'atom ' in line: + self.coord_mode = 'cartesian' + cart = True + atom_read = True + word = line.split() + try: + pos[0] = float(word[1]) + pos[1] = float(word[2]) + pos[2] = float(word[3]) + atom_name = word[4].decode('UTF-8') + except ValueError: + raise GeometryError("Error reading basis coordinate: '" + line + "'") + + if b'atom_frac ' in line: + self.coord_mode = 'direct' + cart = False + atom_read = True + word = line.split() + try: + pos[0] = float(word[1]) + pos[1] = float(word[2]) + pos[2] = float(word[3]) + atom_name = word[4].decode('UTF-8') + except ValueError: + raise GeometryError("Error reading basis coordinate: '" + line + "'") + + if atom_read: + sd_flags = 'T T T' + self.basis.append(Site(cart=cart, position=np.array(pos), + SD_FLAG=sd_flags, occupant=atom_name, occ_alias=atom_name)) + atom_read = False + + # done reading, analyze now + all_names = [] + for a in self.basis: + all_names.append(a.occupant) + + tmp = [] + symbols = [all_names[0]] + for i in range(len(all_names)): + tmp = [] + for j in range(len(symbols)): + tmp.append(symbols[j]) + if tmp.count(all_names[i]) == 0: + symbols.append(all_names[i]) + + nums = {} + for i in range(len(symbols)): + k = 0 + for j in range(len(all_names)): + if symbols[i] == all_names[j]: + k += 1 + nums[symbols[i]] = k + + self.type_atoms = symbols + if len(self.type_atoms) == 0: + raise GeometryError("No atom names found") + try: + self.num_atoms = [nums[sym] for sym in symbols] + except ValueError: + raise GeometryError("Could not read number of each atom type") + if len(self.num_atoms) != len(self.type_atoms): + raise GeometryError("Atom type and number lists are not the same length") + self.type_atoms_alias = list(self.type_atoms) + + self._lattice = np.array(lat) + if self._lattice.shape != (3, 3): + raise GeometryError("Lattice shape error: " + np.array_str(self._lattice)) + self._reciprocal_lattice = 2 * np.pi * np.linalg.inv(np.transpose(self._lattice)) + self.scaling = 1.0 + + del tmp, symbols, nums + + def check_constraints(self, cont, ln, total_lines): + flags = ['T', 'T', 'T'] + + if ln >= total_lines: + return flags[0] + ' ' + flags[1] + ' ' + flags[2] + + next_line = cont[ln] + k = 0 + + if b'initial_moment' in next_line: + raise ValueError('Initial Moments will be added automatically, please remove from geometry.skel.') + + if b'constrain_relaxation ' in next_line: + limit = np.min([3, total_lines - ln]) + self.sel_dyn = True + check_lines = cont[ln:ln+limit] + check = check_lines[k].strip() + while b'constrain_relaxation ' in check: + try: + word = check.split() + if word[1] == b'.true.': + flags = ['F', 'F', 'F'] + if word[1].lower() == b'x': + flags[0] = 'F' + if word[1].lower() == b'y': + flags[1] = 'F' + if word[1].lower() == b'z': + flags[2] = 'F' + except ValueError: + raise GeometryError("Error reading constrints from line '" + check + "'") + + k += 1 + + if k == len(check_lines): + break + check = check_lines[k].strip() + + return flags[0] + ' ' + flags[1] + ' ' + flags[2] + + def read_pos(self, file): + self.basis = Poscar(file).basis + self._lattice = Poscar(file).lattice() + self._reciprocal_lattice = Poscar(file).reciprocal_lattice() + self.coord_mode = Poscar(file).coord_mode + self.SD_FLAG = Poscar(file).SD_FLAG + self.num_atoms = Poscar(file).num_atoms + self.type_atoms = Poscar(file).type_atoms + self.type_atoms_alias = Poscar(file).type_atoms_alias diff --git a/python/casm/casm/aims/io/io.py b/python/casm/casm/aims/io/io.py new file mode 100755 index 000000000..69f0f228e --- /dev/null +++ b/python/casm/casm/aims/io/io.py @@ -0,0 +1,108 @@ +import os +import sys + +from casm.aims.io.steps import Steps +from casm.aims.io.basis import basis_settings, write_basis +from casm.aims.io.parser import Parser +from casm.aims.io.control import Control +from casm.aims.io.kpoints import Kpoints +from casm.aims.io.geometry import Geometry + +AIMS_INPUT_FILE_LIST = ["control.in", "geometry.in"] + +DEFAULT_AIMS_COPY_LIST = [] +DEFAULT_AIMS_MOVE_LIST = [] +DEFAULT_AIMS_GZIP_LIST = ["std.out"] +DEFAULT_AIMS_REMOVE_LIST = [] + + +class AimsIOError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +def job_complete(jobdir=None): + """Return True if aims job at path 'jobdir' is complete""" + if jobdir is None: + jobdir = os.getcwd() + runfile = os.path.join(jobdir, "std.out") + if not os.path.isfile(runfile) and not os.path.isfile(runfile + ".gz"): + return False + if Parser(runfile).complete: + return True + return False + + +def ionic_steps(jobdir=None): + """Find the number of ionic steps completed in 'jobdir'""" + try: + steps = Steps(os.path.join(jobdir, "std.out")) + return len(steps.E) + except IOError: + raise AimsIOError("Could not read number of ionic steps from " + os.path.join(jobdir, "std.out")) + + +def write_aims_input(dirpath, controlfile, prim_posfile, super_posfile, basisfile, strict_kpoints): + """ Update FHI-aims input files in directory 'dirpath' """ + print("Setting up FHI-aims input files:", dirpath) + + print(" Reading basis:", basisfile) + basis_set = basis_settings(basisfile) + + print(" Reading supercell POS:", super_posfile) + super_cell = Geometry(super_posfile) + prim = Geometry(prim_posfile) + + print(" Reading control.in:", controlfile) + super_ctrl = Control(controlfile) + + print(" Parsing k_grid:", controlfile) + prim_kpoints = Kpoints(controlfile) + + print(" Generating supercell KPOINTS") + if strict_kpoints: + super_kpts = prim_kpoints + else: + super_kpts = prim_kpoints.super_kpoints(prim, super_cell) + + # write main input files + print(" Writing supercell geometry.in:", os.path.join(dirpath, 'geometry.in')) + super_cell.write(os.path.join(dirpath, 'geometry.in'), basis_set) + + print(" Writing control.in:", os.path.join(dirpath, 'control.in')) + super_ctrl.write(os.path.join(dirpath, 'control.in')) + + print(" Parsing supercell k_grid and k_offset into control.in:", os.path.join(dirpath, 'control.in')) + super_kpts.parse(os.path.join(dirpath, 'control.in')) + + print(" Adding basis set data to control.in:", os.path.join(dirpath, 'control.in')) + write_basis(os.path.join(dirpath, 'control.in'), basis_set) + + print(" DONE\n") + sys.stdout.flush() + + +def write_abort_scf(mode='e', jobdir=None): + """ Write abort_scf file with two modes: + mode = 'e' for stop at the next electronic step + mode = 'i' for stop at the next ionic step + """ + if jobdir is None: + jobdir = os.getcwd() + if mode.lower()[0] == 'e': + filename = os.path.join(jobdir, 'abort_scf') + stop_string = " " + elif mode.lower()[0] == 'i': + filename = os.path.join(jobdir, 'abort_opt') + stop_string = " " + else: + raise AimsIOError("Invalid abort mode specified: " + str(mode)) + + try: + with open(filename, 'w') as f: + f.write(stop_string) + except IOError as e: + raise e diff --git a/python/casm/casm/aims/io/kpoints.py b/python/casm/casm/aims/io/kpoints.py new file mode 100755 index 000000000..7360638b8 --- /dev/null +++ b/python/casm/casm/aims/io/kpoints.py @@ -0,0 +1,151 @@ +import os +import copy + +import numpy as np + +from casm.aims.io.geometry import Geometry +import casm.aimswrapper +from casm.project import DirectoryStructure, ProjectSettings + + +class KpointsError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Kpoints: + """ + The Kpoints class contains: + self.header: (str) the first line from the KPOINTS file being read + self.num_points: (int) contains the value in the second line (0=>automatic) + self.subdivisions: (list of int) the number of kpoints along each of the vectors in reciprocal space + or the kpoint total "length" if the mode is Automatic + self.automode: (str) Gamma/Monkhorst-Pack/Automatic + self.shift: (list of float) the shifts that are added to the automatically generated points + """ + def __init__(self, filename): + """ Constructs a Kpoints object""" + try: + file = open(filename) + except IOError: + raise KpointsError('IOError' + filename) + + for line in file: + try: + if "k_grid" in line: + self.subdivisions = [float(x) for x in line.split()[1:4]] + else: + pass + except IOError: + raise KpointsError("Error reading k_grid from line: '" + line + "'\nIn file: '" + filename + "'") + + try: + if "k_offset" in line: + self.shift = [float(x) for x in line.split()[1:4]] + else: + pass + except IOError: + raise KpointsError("Error reading k_offset from line: '" + line + "'\nIn file: '" + filename + "'") + file.close() + + def super_kpoints(self, prim, super_cell): + """ Assuming 'self' is the kpoints associated with a PRIM, it uses a scaling method to calculate + the kpoint-mesh for a supercell, such that it has a equal or greater kpoint + density than the prim. + + Returns: + super_kpoints: a Kpoints object for the supercell + + Args: + prim: Poscar object for the prim + super_cell: a Kpoints object for the supercell + """ + configdir = os.getcwd() + _res = os.path.split(configdir) + cfgname = os.path.split(_res[0])[1] + "/" + _res[1] + casm_dirs = DirectoryStructure(configdir) + casm_sets = ProjectSettings(configdir) + clex = casm_sets.default_clex + setfile = casm_dirs.settings_path_crawl("relax.json", cfgname, clex) + settings = casm.aimswrapper.read_settings(setfile) + + super_kpoints = copy.deepcopy(self) + + if prim is None: + raise KpointsError("This error means there is no geometry.skel...") + + super_kpoints.subdivisions = [1, 1, 1] + + # calculate prim volumetric kpoint densities + prim_density = self.density(prim) + + # calculate recip lattice vector lengths + super_recip_vec_lengths = [np.linalg.norm(super_cell.reciprocal_lattice(x)) for x in range(3)] + + # while supercell kpoint density is less than prim kpoint density + while super_kpoints.density(super_cell) < prim_density: + # increase the number of subdivisions along the least dense super recip vector + linear_density = [super_kpoints.subdivisions[x] / super_recip_vec_lengths[x] for x in range(3)] + min_index = linear_density.index(min(linear_density)) + super_kpoints.subdivisions[min_index] += 1 + + # set all subdivisions to be at similar linear density + scale = super_kpoints.subdivisions[min_index] / super_recip_vec_lengths[min_index] + for i in range(2): + super_kpoints.subdivisions[i] = int(np.ceil(scale * super_recip_vec_lengths[i] - 0.1)) + + if settings["is_slab"]: # slab override for FHI-aims ONLY + ox = prim.lattice(0) + oy = prim.lattice(1) + + sx = super_cell.lattice(0) + sy = super_cell.lattice(1) + + norm_prim_x = np.linalg.norm(ox) + norm_prim_y = np.linalg.norm(oy) + + norm_supr_x = np.linalg.norm(sx) + norm_supr_y = np.linalg.norm(sy) + + super_kpoints.subdivisions[0] = int(np.ceil(self.subdivisions[0] * (norm_prim_x / norm_supr_x))) + super_kpoints.subdivisions[1] = int(np.ceil(self.subdivisions[1] * (norm_prim_y / norm_supr_y))) + super_kpoints.subdivisions[2] = 1 + + return super_kpoints + + def density(self, kpts_obj): + """ Return the kpoint density with respect to a Kpoints object. """ + return ((self.subdivisions[0] * self.subdivisions[1] * self.subdivisions[2]) / + Geometry.reciprocal_volume(kpts_obj)) + + def parse(self, filename): + """ Parse k_grid and k_offset """ + try: + file = open(filename, 'r') + except IOError: + raise KpointsError("Write failed") + + tmp = file.readlines() + file.close() + + ln = 0 + for line in tmp: + if "k_grid" in line: + tmp[ln] = "k_grid " + str(self.subdivisions[0]) + " " + str(self.subdivisions[1]) + \ + " " + str(self.subdivisions[2]) + "\n" + if "k_offset" in line: + tmp[ln] = "k_offset " + str(self.shift[0]) + " " + str(self.shift[1]) + \ + " " + str(self.shift[2]) + "\n" + ln += 1 + + file = open(filename, 'w') + for line in tmp: + file.write(line) + file.close() + + del tmp + + return diff --git a/python/casm/casm/aims/io/parser.py b/python/casm/casm/aims/io/parser.py new file mode 100755 index 000000000..fa12a8f9f --- /dev/null +++ b/python/casm/casm/aims/io/parser.py @@ -0,0 +1,44 @@ +import os +import gzip + + +class ParserError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Parser(object): + """Parse FHI-aims output files. + + Currently, just contains: + self.complete = True/False + self.kpoints = list of int, or none + """ + def __init__(self, filename): + self.complete = False + self.read(filename) + + def read(self, file): + if os.path.isfile(file): + if file.split(".")[-1].lower() == "gz": + f = gzip.open(file) + else: + f = open(file) + elif os.path.isfile(file + ".gz"): + f = gzip.open(file + ".gz") + else: + raise ParserError("file not found: " + file) + + for line in f: + try: + if "Have a nice day." in line: + self.complete = True + return True + except IOError: + raise ParserError("Error reading 'Have a nice day.' from line: '" + line + "'\n" + "NOT CONVERGED ERROR In file: '" + file + "'") + + f.close() diff --git a/python/casm/casm/aims/io/steps.py b/python/casm/casm/aims/io/steps.py new file mode 100755 index 000000000..88a9f012f --- /dev/null +++ b/python/casm/casm/aims/io/steps.py @@ -0,0 +1,38 @@ +import os +import gzip + + +class StepsError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Steps(object): + """Parse run file for ionic steps + + Currently, just contains: + self.E = list of ionic step Total energy values + """ + def __init__(self, filename): + self.filename = filename + self.E = [] + self.read() + + def read(self): + """Parse file. Currently just collects E0 as list 'self.E'""" + if os.path.isfile(self.filename): + if self.filename.split(".")[-1].lower() == "gz": + f = gzip.open(self.filename) + else: + f = gzip.open(self.filename + ".gz") + else: + raise StepsError("file not found: " + self.filename) + + self.E = [] + for line in f: + if b"| Total energy :" in line: + self.E.append(float(line.split()[6])) + f.close() diff --git a/python/casm/casm/aims/relax.py b/python/casm/casm/aims/relax.py new file mode 100755 index 000000000..36d6c4cb4 --- /dev/null +++ b/python/casm/casm/aims/relax.py @@ -0,0 +1,315 @@ +import os +import sys + +import casm +import casm.project +import casm.aimswrapper + +from casm.aims.aims import continue_job, run, complete_job +from casm.aims.io.basis import basis_settings +from casm.aims.io.geometry import Geometry +from casm.aims.io.io import AIMS_INPUT_FILE_LIST, job_complete +from casm.aims.io.parser import Parser + + +class AimsRelaxError: + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class AimsRelax(object): + """The Relax class contains functions for setting up, executing, and parsing a elaxation. + + The relaxation is initialized in a directory containing input files, called 'relaxdir'. + It then creates the following directory structure: + .../relaxdir/ + run.0/ + run.1/ + ... + run.final/ + + 'run.i' directories are only created when ready. + 'run.final' contains the final run + + Contains: + self.relaxdir (.../relax) + self.rundir (list of .../relax/run.i) + self.finaldir (.../relax/run.final) + """ + def __init__(self, relaxdir=None, settings=None): + """ + Construct a relaxation job object. + + Args: + relaxdir: path to relaxation directory + settings: dictionary-like object containing settings, or if None, it reads + the json file: .../relaxdir/relax.json + + possible settings keys are: + used by .run() function: + "ncpus": number of ncpus to run mpi on + "aims_cmd": (default "aims") shell command to execute FHI-aims + or None to use default mpirun + "strict_kpoint": force strict copying of kpoints in control.skel file, + otherwise kpoints are scaled based on supercell size + used by not_converging(): + "run_limit": (default 10) maximum number of runs to allow + before setting status to "not_converging" + """ + print("Constructing a FHI-aims Relax object") + sys.stdout.flush() + + # store path to .../relaxdir, and create if not existing + if relaxdir is None: + relaxdir = os.getcwd() + self.relaxdir = os.path.abspath(relaxdir) + + print(" Relax directory:", self.relaxdir) + sys.stdout.flush() + + # find existing .../relaxdir/run.run_index directories, store paths in self.rundir list + self.rundir = [] + self.errdir = [] + self.update_rundir() + self.update_errdir() + + if settings is None: + self.settings = dict() + else: + self.settings = settings + + # set default settings: + if "aims_cmd" not in self.settings: + self.settings["aims_cmd"] = None + if "compress" not in self.settings: + self.settings["compress"] = [] + + storedir = 'run.' + str(settings["basis"]) + self.finaldir = os.path.join(self.relaxdir, storedir) + + print("FHI-aims Relax object constructed\n") + sys.stdout.flush() + + def add_rundir(self): + """Make a new run.i directory""" + os.mkdir(os.path.join(self.relaxdir, "run." + str(len(self.rundir)))) + self.update_rundir() + self.update_errdir() + + def update_rundir(self): + """Find all .../config/vasp/relax/run.i directories, store paths in self.rundir list""" + self.rundir = [] + run_index = len(self.rundir) + while os.path.isdir(os.path.join(self.relaxdir, "run." + str(run_index))): + self.rundir.append(os.path.join(self.relaxdir, "run." + str(run_index))) + run_index += 1 + + def add_errdir(self): + """Move run.i to run.i_err.j directory""" + os.rename(self.rundir[-1], self.rundir[-1] + "_err." + str(len(self.errdir))) + self.update_errdir() + + def update_errdir(self): + """Find all .../config/vasp/relax/run.i_err.j directories, store paths in self.errdir list""" + self.errdir = [] + if len(self.rundir) == 0: + pass + else: + err_index = len(self.errdir) + while os.path.isdir(self.rundir[-1] + "_err." + str(err_index)): + self.errdir.append(self.rundir[-1] + "_err." + str(err_index)) + err_index += 1 + + def setup(self, initdir, settings): + """ mv all files and directories (besides initdir) into initdir """ + print("Moving files into initial run directory:", initdir) + + if settings["basis"] == 'light': + initdir = os.path.abspath(initdir) + for p in os.listdir(self.relaxdir): + if p in AIMS_INPUT_FILE_LIST and os.path.join(self.relaxdir, p) != initdir: + os.rename(os.path.join(self.relaxdir, p), os.path.join(initdir, p)) + print("\n") + + if settings["basis"] == 'tight': + configdir = os.getcwd() + _res = os.path.split(configdir) + cfgname = os.path.split(_res[0])[1] + "/" + _res[1] + casm_dirs = casm.project.DirectoryStructure(configdir) + casm_sets = casm.project.ProjectSettings(configdir) + clex = casm_sets.default_clex + aimsfiles = casm.aimswrapper.aims_input_file_names(casm_dirs, cfgname, clex) + bin_0, bin_1, bin_2, basisfile = aimsfiles + + print('Taking geometry.in.next_step from previous light run...') + + old_light_dir = os.path.join(self.relaxdir, 'run.light') + old_geo_file = os.path.join(old_light_dir, 'geometry.in.next_step') + + if not os.path.isdir(old_light_dir): + if not os.path.isfile(old_geo_file): + err_str = 'The directory ', old_light_dir, ' exists which points to a previously converged run.' + err_str += 'However, I cannot find ', old_geo_file, ' which usually means that the' \ + ' starting structure' + err_str += 'was already converged. This is very unusual and you need to check what went on.' + err_str += 'When you are sure that the structure converged properly, ' \ + 'copy the correct geometry file to ', old_geo_file, ' and restart...' + raise AimsRelaxError(err_str) + err_str = 'Previous converged light run does not exist, not running tight on unrelaxed structures...' + err_str += 'Stopping, please run light first' + raise AimsRelaxError(err_str) + + my_basis = basis_settings(basisfile) + newgeo = Geometry(os.path.join(old_light_dir, "geometry.in.next_step")) + newgeo.write(os.path.join(self.relaxdir, "geometry.in"), my_basis) + initdir = os.path.abspath(initdir) + for p in os.listdir(self.relaxdir): + if p in AIMS_INPUT_FILE_LIST and os.path.join(self.relaxdir, p) != initdir: + os.rename(os.path.join(self.relaxdir, p), os.path.join(initdir, p)) + + print('') + sys.stdout.flush() + + def complete(self): + """Check if the relaxation is complete. + + Completion criteria: "Have a nice day." in std.out + """ + outfile = os.path.join(self.finaldir, "std.out") + if not os.path.isfile(outfile): + return False + if not Parser(outfile).complete(): + return False + return True + + def converged(self): + """Check if configuration is relaxed. + + This is called when self.rundir[-1] is complete. + + Convergence criteria: FHI-aims returns "Have a nice day." + """ + outcarfile = os.path.join(self.rundir[-1], "std.out") + if not os.path.isfile(outcarfile): + return False + if not Parser(outcarfile).complete(): + return False + return True + + def not_converging(self): + """Check if configuration is not converging. + + This is called when self.rundir[-1] is complete and not a constant volume job and self.converged() == False. + + Not converging criteria: >= 10 runs without completion + """ + if len(self.rundir) >= int(self.settings["run_limit"]): + return True + return False + + def run(self): + """ Perform a series of FHI-aims jobs to relax a structure. """ + + print("Begin FHI-aims relaxation run") + sys.stdout.flush() + + # get current status of the relaxation: + status, task = self.status() + print("\n++ status:" + status + " next task:", task) + sys.stdout.flush() + + while status == "incomplete": + if task == "setup": + self.add_rundir() + self.setup(self.rundir[-1], self.settings) + elif task == "relax": + self.add_rundir() + continue_job(self.rundir[-2], self.rundir[-1], self.settings) + elif task == "constant": + self.add_rundir() + continue_job(self.rundir[-2], self.rundir[-1], self.settings) + else: + # probably hit walltime + print('Assuming walltime has been hit last run, continuing there...') + self.add_rundir() + continue_job(self.rundir[-2], self.rundir[-1], self.settings) + + while True: + # run FHI-aims + result = run(self.rundir[-1], ncpus=self.settings["ncpus"]) + # if no errors, continue + if result is None or result == self.not_converging: + break + # else, attempt to fix first error + self.add_errdir() + os.mkdir(self.rundir[-1]) + # self.add_rundir() +# err = result.itervalues().next() + +# print "\n++ status:", "error", " next task:", "fix_error" +# sys.stdout.flush() + +# print "Attempting to fix error:", str(err) +# err.fix(self.errdir[-1],self.rundir[-1], self.settings) +# print "" +# sys.stdout.flush() + +# if (self.settings["backup"] != None) and len(self.rundir) > 1: +# print "Restoring from backups:" +# for f in self.settings["backup"]: +# if os.path.isfile(os.path.join(self.rundir[-2], f + "_BACKUP.gz")): +# f_in = gzip.open(os.path.join(self.rundir[-2], f + "_BACKUP.gz", 'rb')) +# f_out = open(os.path.join(self.rundir[-1], f, 'wb')) +# f_out.write(f_in.read()) +# f_in.close() +# f_out.close() +# print f, " restored!" +# sys.stdout.flush() + + status, task = self.status() + print("\n++ status:" + status + " next task:", task) + sys.stdout.flush() + + if status == "complete": + if not os.path.isdir(self.finaldir): + # mv final results to relax.final + print("mv" + os.path.basename(self.rundir[-1]) + os.path.basename(self.finaldir)) + sys.stdout.flush() + os.rename(self.rundir[-1], self.finaldir) + self.rundir.pop() + complete_job(self.finaldir, self.settings) + + return status, task + + def status(self): + """ Determine the status of a vasp relaxation series of runs. Individual runs in the series + are in directories labeled "run.0", "run.1", etc. + + Returns a tuple: (status = "incomplete" or "complete" or "not_converging", + task = continuedir or "relax" or "constant" or None) + + The first value is the status of the entire relaxation. + + The second value is the current task, where 'continuedir' is the path to a + vasp job directory that is not yet completed, "relax" indicates another + volume relaxation job is required, and "constant" that a constant volume run is required. + """ + # if not yet started + if len(self.rundir) == 0: + return "incomplete", "setup" + + # check status of relaxation runs + self.update_rundir() + + # if the latest run is complete: + if job_complete(self.rundir[-1]): + return "complete", None + + # elif not converging, return 'not_converging' error + elif self.not_converging(): + return "not_converging", None + else: + return "incomplete", "relax" diff --git a/python/casm/casm/aimswrapper/__init__.py b/python/casm/casm/aimswrapper/__init__.py new file mode 100644 index 000000000..2963d91e2 --- /dev/null +++ b/python/casm/casm/aimswrapper/__init__.py @@ -0,0 +1,13 @@ +"""A wrapper for running FHI-aims through casm""" + +from casm.aimswrapper.aimswrapper import AimsWrapperError, read_settings, write_settings, \ + aims_input_file_names +from casm.aimswrapper.relax import Relax + +__all__ = [ + 'Relax', + 'AimsWrapperError', + 'read_settings', + 'write_settings', + 'aims_input_file_names' +] \ No newline at end of file diff --git a/python/casm/casm/aimswrapper/aimswrapper.py b/python/casm/casm/aimswrapper/aimswrapper.py new file mode 100755 index 000000000..dd6722fcc --- /dev/null +++ b/python/casm/casm/aimswrapper/aimswrapper.py @@ -0,0 +1,151 @@ +import json + +from casm.project.io import read_project_settings + +from casm.aims.io.io import DEFAULT_AIMS_COPY_LIST, DEFAULT_AIMS_MOVE_LIST + + +class AimsWrapperError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +def read_settings(filename): + """Returns a JSON object reading JSON files containing settings for FHI-aims PBS jobs. + + Returns: + settings = a JSON object containing the settings file contents + This can be accessed like a dict: settings["account"], etc. + ** All values are expected to be 'str' type. ** + + The required keys are: + "queue": queue to submit job in + "ppn": processors (cores) per node to request + "atom_per_proc": max number of atoms per processor (core) + "walltime": walltime to request (ex. "48:00:00") + + The optional keys are: + "account": account to submit job under (default None) + "pmem": string for requested memory (default None) + "priority": requested job priority (default "0") + "message": when to send messages about jobs (ex. "abe", default "a") + "email": where to send messages (ex. "me@fake.com", default None) + "qos": quality of service, 'qos' option (ex. "fluxoe") + "aims_cmd": FHI-aims execution command (default is "aims" (ncpus=1) or "mpirun -np {NCPUS} aims" (ncpus>1)) + "ncpus": number of cpus (cores) to run on (default $PBS_NP) + "run_limit": number of vasp runs until "not_converging" (default 10) + "err_types" : list of errors to check for, currently ony FrozenError + "prerun" : bash commands to run before aims.Relax.run (default None) + "postrun" : bash commands to run after aims.Relax.run completes (default None) + """ + try: + file = open(filename) + settings = read_project_settings(filename) + file.close() + except IOError as e: + print("Error reading settings file:", filename) + raise e + + required = ["is_slab", "basis"] + optional = ["fix_pos"] + + for key in required: + if key not in settings: + raise AimsWrapperError(key + "' missing from: '" + filename + "'") + + for key in optional: + if key not in settings: + if key.lower() in ["extra_input_files", "remove", "compress", "backup"]: + settings[key] = [] + elif key.lower() in ["move"]: + settings[key] = DEFAULT_AIMS_MOVE_LIST + elif key.lower() in ["copy"]: + settings[key] = DEFAULT_AIMS_COPY_LIST + else: + settings[key] = None + + if settings["basis"] != "light" and settings["basis"] != "tight": + raise AimsWrapperError("Basis setting for FHI-aims must be >light< or >tight<") + + if settings["is_slab"] == "True": + if 'fix_pos' not in settings: + err_str = "Error: When running a slab you need to define a 'fix_pos' key" + err_str += "which defines the position below which atoms in the slab will be fixed (all coordinates)." + err_str += "fix_pos: Units in absolute Angstrom coordinates (NOT fractional!)" + raise AimsWrapperError(err_str) + settings["is_slab"] = True + else: + settings["is_slab"] = False + + settings['strict_kpoints'] = False + + return settings + + +def write_settings(settings, filename): + """ Write 'settings' as json file, 'filename' """ + with open(filename, 'w') as f: + json.dump(settings, f, indent=4) + + +def aims_input_file_names(workdir, configname, clex): + """ + Collect casm.aimswrapper input files from the CASM project hierarchy + + Looks for: + + control.in: + The base input file used for calculations. Found via: + DirectoryStructure.settings_path_crawl + + POS: + The CASM-generated POS file giving the initial structure to be calculated. + + SPECIES: + The SPECIES file specifying FHI-aims basis settings for each species in the structure. + + Arguments + --------- + + workdir: casm.project.DirectoryStructure instance + CASM project directory hierarchy + + configname: str + The name of the configuration to be calculated + + clex: casm.project.ClexDescription instance + The cluster expansion being worked on. Used for the 'calctype' settings. + + + Returns + ------- + + filepaths: tuple(control.in, POS, SPECIES) + A tuple containing the paths to the aimswrapper input files + + + Raises + ------ + If any required file is not found. + + """ + # Find required input files in CASM project directory tree + controlfile = workdir.settings_path_crawl("control.skel", configname, clex) + prim_geometryfile = workdir.settings_path_crawl("geometry.skel", configname, clex) + super_poscarfile = workdir.POS(configname) + basisfile = workdir.settings_path_crawl("basis", configname, clex) + + # Verify that required input files exist + if controlfile is None: + raise AimsWrapperError("aims_input_file_names failed. No control.skel file found in CASM project.") + if prim_geometryfile is None: + raise AimsWrapperError("No reference geometry.skel found in CASM project.") + if super_poscarfile is None: + raise AimsWrapperError("aims_input_file_names failed. No POS file found for this configuration.") + if basisfile is None: + raise AimsWrapperError("aims_input_file_names failed. No SPECIES file found in CASM project.") + + return controlfile, prim_geometryfile, super_poscarfile, basisfile diff --git a/python/casm/casm/aimswrapper/relax.py b/python/casm/casm/aimswrapper/relax.py new file mode 100755 index 000000000..6dbfa2950 --- /dev/null +++ b/python/casm/casm/aimswrapper/relax.py @@ -0,0 +1,440 @@ +import os +import sys +import json + +from prisms_jobs import Job, JobDB, error_job, complete_job, JobsError, JobDBError, EligibilityError + +from casm.wrapper.misc import confname_as_jobname + +from casm.project import DirectoryStructure, ProjectSettings +from casm.misc.noindent import NoIndent, NoIndentEncoder + +from casm.aims.relax import AimsRelax +from casm.aims.io.io import write_aims_input +from casm.aimswrapper import AimsWrapperError, read_settings, write_settings, aims_input_file_names +from casm.aims.io.aimsrun import AimsRun +from casm.aims.io.geometry import Geometry + + +class Relax(object): + """The Relax class contains functions for setting up, executing, and parsing a relaxation. + + The relaxation creates the following directory structure: + config/ + calctype.name/ + run.0/ + run.1/ + ... + run.final/ + + 'run.i' directories are only created when ready. + 'run.final' is a final run + + This automatically looks for the settings files using casm.project.DirectoryStructure.settings_path_crawl + + Attributes + ---------- + + casm_settings: casm.project.ProjectSettings instance + CASM project settings + + casm_directories: casm.project.DirectoryStructure instance + CASM project directory hierarchy + + settings: dict + Settings for pbs and the relaxation, see casm.project.io.read_project_settings + + configdir: str + Directory where configuration results are stored. The result of: + casm.project.DirectoryStructure.configuration_dir(self.configname) + + configname: str + The name of the configuration to be calculated + + sort: boolean + True if sorting atoms in POSCAR by type + + clex: casm.project.ClexDescription instance + The cluster expansion being worked on. Used for the 'calctype' settings. + Currently, fixed to self.casm_settings.default_clex. + + """ + def __init__(self, configdir=None, auto=True, sort=True): + """ + Construct a relaxation job object. + + Arguments + ---------- + + configdir: str, optional, default=None + Path to configuration directory. If None, uses the current working directory + + sort: boolean, optional, default=True, + Use True to sort atoms in POSCAR by type + + """ + print("Construct a casm.aimswrapper.Relax instance:") + + if configdir is None: + configdir = os.getcwd() + print(" Input directory:" + configdir) + + # get the configname from the configdir path + _res = os.path.split(configdir) + self.configname = os.path.split(_res[0])[1] + "/" + _res[1] + print(" Configuration:" + self.configname) + + print(" Reading CASM settings") + self.casm_directories = DirectoryStructure(configdir) + self.casm_settings = ProjectSettings(configdir) + if self.casm_settings is None: + raise AimsWrapperError("Not in a CASM project. The '.casm' directory was not found.") + + if os.path.abspath(configdir) != self.configdir: + print("") + print("input configdir:" + configdir) + print("determined configname:" + self.configname) + print("expected configdir given configname:" + self.configdir) + raise AimsWrapperError("Mismatch between configname and configdir") + + # fixed to default_clex for now + self.clex = self.casm_settings.default_clex + + # store path to .../config/calctype.name, and create if not existing + self.calcdir = self.casm_directories.calctype_dir(self.configname, self.clex) + + print(" Calculations directory:" + self.calcdir) + if not os.path.isdir(self.calcdir): + os.mkdir(self.calcdir) + + # read the settings json file + print(" Reading relax.json...") + sys.stdout.flush() + setfile = self.casm_directories.settings_path_crawl("relax.json", self.configname, self.clex) + + if setfile is None: + raise AimsWrapperError("Could not find relax.json in settings directory") + else: + print(" Read settings from:" + setfile) + self.settings = read_settings(setfile) + + # set default settings if not present + if "aims_cmd" not in self.settings: + self.settings["aims_cmd"] = None + if "ppn" not in self.settings: + self.settings["ppn"] = None + if "nodes" not in self.settings: + self.settings["nodes"] = None + if "run_limit" not in self.settings: + self.settings["run_limit"] = None + if "prerun" not in self.settings: + self.settings["prerun"] = None + if "postrun" not in self.settings: + self.settings["postrun"] = None + + self.auto = auto + self.sort = sort + self.strict_kpoints = self.settings['strict_kpoints'] + + print(" DONE\n") + sys.stdout.flush() + + @property + def configdir(self): + return self.casm_directories.configuration_dir(self.configname) + + def setup(self): + """ Setup initial relaxation run + Uses the following files from the most local .../settings/calctype.name directory: + control.skel: input settings + geometry.skel: reference for KPOINTS object + basis: species info + Uses the following files from the .../config directory: + POS: structure of the configuration to be relaxed + """ + # Find required input files in CASM project directory tree + aimsfiles = aims_input_file_names(self.casm_directories, self.configname, self.clex) + controlfile, prim_posfile, super_posfile, basisfile = aimsfiles + + write_aims_input(self.calcdir, controlfile, prim_posfile, super_posfile, basisfile, self.strict_kpoints) + + def submit(self): + """Submit a job for this relaxation""" + print("Submitting configuration: " + self.configname) + # first, check if the job has already been submitted and is not completed + db = JobDB() + print("Calculation directory: ", self.calcdir) + sub_id = db.select_regex_id("rundir", self.calcdir) + + if sub_id is not []: + for j in sub_id: + job = db.select_job(j) + if job["jobstatus"] != "?": + print("JobID: " + job["jobid"], + " Jobstatus:" + job["jobstatus"] + " Not submitting.") + sys.stdout.flush() + return + + # second, only submit a job if relaxation status is "incomplete" + # construct the Relax object + relaxation = AimsRelax(self.calcdir, self.run_settings()) + # check the current status + status, task = relaxation.status() + + if status == "complete": + print("Status:", status, " Not submitting.") + sys.stdout.flush() + + # ensure job marked as complete in db + if self.auto: + for j in sub_id: + job = db.select_job(j) + if job["taskstatus"] == "Incomplete": + try: + complete_job(jobid=j) + except IOError as e: + raise AimsWrapperError(e) + + # ensure results report written + if not os.path.isfile(os.path.join(self.calcdir, "properties.calc.json")): + self.finalize() + + return + + elif status == "not_converging": + print("Status:" + status + " Not submitting.") + sys.stdout.flush() + return + + elif status != "incomplete": + raise AimsWrapperError("unexpected relaxation status: '" + + status + "' and task: '" + task + "'") + + print("Preparing to submit a FHI-aims relaxation job") + sys.stdout.flush() + + # cd to configdir, submit jobs from configdir, then cd back to currdir + currdir = os.getcwd() + os.chdir(self.calcdir) + + # determine the number of atoms in the configuration + print("Counting atoms in the Supercell") + sys.stdout.flush() +# pos = vasp.io.Poscar(os.path.join(self.configdir, "POS")) +# N = len(pos.basis) + + # construct command to be run + cmd = "" + if self.settings["prerun"] is not None: + cmd += self.settings["prerun"] + "\n" + cmd += "python -c \"import casm.aimswrapper; casm.aimswrapper.Relax('" + self.configdir + "').run()\"\n" + if self.settings["postrun"] is not None: + cmd += self.settings["postrun"] + "\n" + + print("Constructing the job") + sys.stdout.flush() + + # construct a pbs.Job + job = Job(name=confname_as_jobname(self.configdir), + account=self.settings["account"], + nodes=int(self.settings["nodes"]), + ppn=int(self.settings['ppn']), + walltime=self.settings["walltime"], + pmem=self.settings["pmem"], + qos=self.settings["qos"], + queue=self.settings["queue"], + message=self.settings["message"], + email=self.settings["email"], + priority=self.settings["priority"], + command=cmd, + auto=self.auto) + + print("Submitting") + sys.stdout.flush() + # submit the job + job.submit() + self.report_status("submitted") + + # return to current directory + os.chdir(currdir) + + print("CASM AimsWrapper relaxation job submission complete\n") + sys.stdout.flush() + + def run_settings(self): + """ Set default values based on runtime environment""" + settings = dict(self.settings) + + # set default values + if settings["ncpus"] is None: + settings["ncpus"] = int(settings["nodes"]) * int(settings["ppn"]) + + if settings["run_limit"] is None or settings["run_limit"] == "CASM_DEFAULT": + settings["run_limit"] = 10 + + return settings + + def run(self): + """ Setup input files, run a FHI-aims relaxation, and report results """ + # construct the Relax object + relaxation = AimsRelax(self.calcdir, self.run_settings()) + + # check the current status + status, task = relaxation.status() + + if status == "complete": + print("Status: " + status) + sys.stdout.flush() + # mark job as complete in db + if self.auto: + try: + complete_job() + except IOError as e: + raise AimsWrapperError(e) + + # write results to properties.calc.json + self.finalize() + return + + elif status == "not_converging": + print("Status:" + status) + self.report_status("failed", "run_limit") + print("Returning") + sys.stdout.flush() + return + + elif status == "incomplete": + if task == "setup": + self.setup() + self.report_status("started") + status, task = relaxation.run() + + else: + self.report_status("failed", "unknown") + raise AimsWrapperError("unexpected relaxation status: '" + status + "' and task: '" + task + "'") + + # once the run is done, update database records accordingly + if status == "not_converging": + # mark error + if self.auto: + try: + error_job("Not converging") + except (JobsError, JobDBError) as e: + raise AimsWrapperError(e) + + print("Not Converging!") + sys.stdout.flush() + self.report_status("failed", "run_limit") + + # print a local settings file, so that the run_limit can be extended if the + # convergence problems are fixed + + config_set_dir = self.casm_directories.configuration_calc_settings_dir(self.configname, self.clex) + + try: + os.makedirs(config_set_dir) + except IOError: + pass + settingsfile = os.path.join(config_set_dir, "relax.json") + write_settings(self.settings, settingsfile) + + print("Writing:" + settingsfile) + print("Edit the 'run_limit' property if you wish to continue.") + sys.stdout.flush() + return + + elif status == "complete": + # mark job as complete in db + try: + complete_job() + except (JobsError, JobDBError, EligibilityError) as e: + raise AimsWrapperError(e) + self.finalize() + + else: + self.report_status("failed", "unknown") + raise AimsWrapperError("Relaxation complete with unexpected status: '" + + status + "' and task: '" + task + "'") + + def report_status(self, status, failure_type=None): + """Report calculation status to status.json file in configuration directory. + + Args: + status: string describing calculation status. Currently used values are + not_submitted + submitted + complete + failed + failure_type: optional string describing reason for failure. Currently used values are + unknown + electronic_convergence + run_limit""" + + output = dict() + output["status"] = status +# output["basis_type"] = basis_type + if failure_type is not None: + output["failure_type"] = failure_type + + outputfile = os.path.join(self.calcdir, "status.json") + with open(outputfile, 'w') as file: + file.write(json.dumps(output, cls=NoIndentEncoder, indent=4, sort_keys=True)) + print("Wrote " + outputfile) + sys.stdout.flush() + + def finalize(self): + if self.is_converged(): + # write properties.calc.json + finaldir = "run." + str(self.settings["basis"]) + aimsdir = os.path.join(self.calcdir, finaldir) + super_posfile = os.path.join(self.configdir, "POS") + speciesfile = self.casm_directories.settings_path_crawl("basis", self.configname, self.clex) + output = self.properties(aimsdir, super_posfile, speciesfile) + outputfile = os.path.join(self.calcdir, "properties.calc.json") + with open(outputfile, 'w') as file: + file.write(json.dumps(output, cls=NoIndentEncoder, indent=4, sort_keys=True)) + print("Wrote " + outputfile) + sys.stdout.flush() + self.report_status(status='complete') + + def is_converged(self): + # Check for electronic convergence in completed calculations. Returns True or False. + relaxation = AimsRelax(self.calcdir, self.run_settings()) + rundirname = 'run.' + str(self.settings["basis"]) + dname = os.path.join(self.configdir, 'calctype.default', rundirname) + relaxation.rundir.append(dname) + for i in range(len(relaxation.rundir)): + try: + with open(os.path.join(self.calcdir, relaxation.rundir[-i-1], "std.out")) as runfile: + for line in runfile: + if "Have a nice day." in line: + return True + except IOError: + pass + + return False + + @staticmethod + def properties(aimsdir, super_posfile=None, basisfile=None): + """Report results to properties.calc.json file in configuration directory, + after checking for electronic convergence.""" + output = dict() + arun = AimsRun(os.path.join(aimsdir, "std.out")) + + if super_posfile is not None and basisfile is not None: + # basis_settings_loc = basis_settings(basisfile) + super_cell = Geometry(super_posfile) + # unsort_dict = super_cell.unsort_dict() + else: + # fake unsort_dict (unsort_dict[i] == i) + # unsort_dict = dict(zip(range(0, len(arun.basis)), range(0, len(arun.basis)))) + super_cell = Geometry(os.path.join(aimsdir, "geometry.in")) + + output["atom_type"] = super_cell.type_atoms + output["atoms_per_type"] = super_cell.num_atoms + output["coord_mode"] = arun.coord_mode + output["relaxed_forces"] = [NoIndent(v) for v in arun.forces] + output["relaxed_lattice"] = [NoIndent(v) for v in arun.lattice] + output["relaxed_basis"] = [NoIndent(v) for v in arun.basis] + output["relaxed_energy"] = arun.total_energy + return output diff --git a/python/casm/casm/feff/__init__.py b/python/casm/casm/feff/__init__.py new file mode 100644 index 000000000..8bb243a6b --- /dev/null +++ b/python/casm/casm/feff/__init__.py @@ -0,0 +1,3 @@ +from casm.feff.feff import Feff, check_consistent_settings +__all__ = ['Feff', + 'check_consistent_settings'] diff --git a/python/casm/casm/feff/feff.py b/python/casm/casm/feff/feff.py new file mode 100644 index 000000000..3a8668f01 --- /dev/null +++ b/python/casm/casm/feff/feff.py @@ -0,0 +1,560 @@ +import re +import os +import sys +import time +import json +import subprocess + +import numpy as np +import matplotlib.pyplot as plt +import scipy.interpolate as sci_int + +from casm.project import DirectoryStructure, ProjectSettings +from casm.project.io import read_project_settings, read_feff_settings + +from pymatgen.io.feff.sets import MPXANESSet +from pymatgen.core.structure import Structure +from pymatgen.symmetry.analyzer import SpacegroupAnalyzer + +from prisms_jobs import Job, JobDB + +import matplotlib +matplotlib.use('Agg') + + +class FeffError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Feff(object): + """ + Computes band structure with pymatgen standard settings + - k-point density is 1000 / atom + - high symmetry path is determined from cell geometry + """ + + def __init__(self, rundir=None): + if rundir is None: + raise FeffError('Can not create from nothing-directory') + + print("Constructing a FEFF XAS object") + sys.stdout.flush() + + print(" Reading CASM settings") + self.casm_directories = DirectoryStructure(rundir) + self.casm_settings = ProjectSettings(rundir) + if self.casm_settings is None: + raise FeffError("Not in a CASM project. The '.casm' directory was not found.") + + # fixed to default_clex for now + self.clex = self.casm_settings.default_clex + + _res = os.path.split(rundir) + self.configname = os.path.split(_res[0])[1] + "/" + _res[1] + + print(" Reading DFT and plot settings for configuration: " + self.configname) + sys.stdout.flush() + + setfile = self.casm_directories.settings_path_crawl("relax.json", self.configname, self.clex) + if setfile is None: + raise FeffError("Could not find relax.json in settings directory") + else: + print(" Read DFT settings from:" + setfile) + self.settings = read_project_settings(setfile) + + fefffile = self.casm_directories.settings_path_crawl("feff.json", self.configname, self.clex) + if fefffile is None: + raise FeffError("Could not find feff.json in settings directory") + else: + print(" Read FEFF settings from:" + fefffile) + self.feff_settings = read_feff_settings(fefffile) + + self.submit_dir = rundir + self.feff_dir = os.path.abspath(os.path.join(rundir, 'calctype.default', 'xanes')) + self.contcar_dir = os.path.abspath(os.path.join(rundir, 'calctype.default', 'run.final')) + self.new_incar = [] + self.status_file = os.path.abspath(os.path.join(rundir, 'calctype.default', 'status.json')) + + self.feff_program_sequence = ['rdinp', 'atomic', 'dmdw', 'pot', 'ldos', 'screen', 'opconsat', + 'xsph', 'fms', 'mkgtr', 'path', 'genfmt', 'ff2x', 'sfconv', + 'compton', 'eels'] + + print(" Run directory: %s" % self.feff_dir) + sys.stdout.flush() + + def setup(self): + print('Getting structure for FEFF computation in directory: ') + print(' %s' % self.feff_dir) + + st = Structure.from_file(os.path.join(self.contcar_dir, 'CONTCAR')) + + to_compute = [] + sp_compute = [] + + for s in SpacegroupAnalyzer(st).get_symmetry_dataset()['equivalent_atoms']: + if to_compute.count(str(s)) == 0: + to_compute.append(str(s)) + sp_compute.append(str(st.species[int(s)].symbol)) + + for ab in range(len(to_compute)): + edge = ['K'] + + for e in edge: + self.feff_settings['feff_user_tags']['EDGE'] = e + mp_xanes = MPXANESSet(absorbing_atom=int(to_compute[ab]), structure=st, + nkpts=self.feff_settings['feff_nkpts'], + radius=self.feff_settings['feff_radius'], + user_tag_settings=self.feff_settings['feff_user_tags']) + + mp_xanes.write_input(os.path.join(self.feff_dir, str(sp_compute[ab]) + '_' + + str(e) + '_' + str(to_compute[ab])), make_dir_if_not_present=True) + + def exec_feff(self, jobdir=None, poll_check_time=5.0, err_check_time=60.0): + """ Run FEFF program sequence + + The 'command' is executed in the directory 'jobdir'. + + Args: + jobdir: directory to run it + poll_check_time: how frequently to check if the job is completed + err_check_time: how frequently to check for freeze + """ + print("Begin FEFF run:") + sys.stdout.flush() + + if jobdir is None: + raise FeffError('Can NOT run FEFF in higher directories') + + currdir = os.getcwd() + os.chdir(jobdir) + + if "PBS_NP" in os.environ: + ncpus = int(os.environ["PBS_NP"]) + elif "SLURM_NPROCS" in os.environ: + ncpus = int(os.environ["SLURM_NPROCS"]) + else: + ncpus = 1 + + err = None + + for p in self.feff_program_sequence: + + if ncpus == 1: + command = os.path.join(self.feff_settings['feff_base_dir'], p) + else: + command = "mpirun -np {NCPUS} " + os.path.join(self.feff_settings['feff_base_mpi'], p) + + if re.search("NCPUS", command): + command = command.format(NCPUS=str(ncpus)) + + print(" jobdir:", jobdir) + print(" exec:", command) + sys.stdout.flush() + + stdout = p + '.log' + stderr = p + '.err' + + sout = open(os.path.join(jobdir, stdout), 'w') + serr = open(os.path.join(jobdir, stderr), 'w') + + p = subprocess.Popen(command.split(), stdout=sout, stderr=serr) + + # wait for process to end, and periodically check for errors + last_check = time.time() + while p.poll() is None: + time.sleep(poll_check_time) + if time.time() - last_check > err_check_time: + last_check = time.time() + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + # FreezeErrors are fatal and usually not helped with abort_scf + if "FreezeError" in err.keys(): + print(" FEFF run seems frozen, killing job") + sys.stdout.flush() + p.kill() + + # close output files + sout.close() + serr.close() + + os.chdir(currdir) + + print("Run ended") + sys.stdout.flush() + + return err + + def run(self): + st = Structure.from_file(os.path.join(self.contcar_dir, 'CONTCAR')) + + to_compute = [] + sp_compute = [] + + for s in SpacegroupAnalyzer(st).get_symmetry_dataset()['equivalent_atoms']: + if to_compute.count(str(s)) == 0: + to_compute.append(str(s)) + sp_compute.append(str(st.species[int(s)].symbol)) + + for ab in range(len(to_compute)): + edge = ['K'] + + for e in edge: + self.feff_settings['feff_user_tags']['EDGE'] = e + mp_xanes = MPXANESSet(absorbing_atom=int(to_compute[ab]), structure=st, + nkpts=self.feff_settings['feff_nkpts'], + radius=self.feff_settings['feff_radius'], + user_tag_settings=self.feff_settings['feff_user_tags']) + + mp_xanes.write_input(os.path.join(self.feff_dir, str(sp_compute[ab]) + '_' + + str(e) + '_' + str(to_compute[ab])), + make_dir_if_not_present=True) + self.exec_feff(jobdir=os.path.join(self.feff_dir, str(sp_compute[ab]) + '_' + + str(e) + '_' + str(to_compute[ab]))) + + self.plot_feff(plot_dir=self.feff_dir) + + def submit(self): + """Submit a job for this band structure computation""" + print("Submitting configuration: " + self.configname) + # first, check if the job has already been submitted and is not completed + db = JobDB() + print("Calculation directory: ", self.feff_dir) + sub_id = db.select_regex_id("rundir", self.feff_dir) + + if sub_id is not []: + for j in sub_id: + job = db.select_job(j) + if job["jobstatus"] != "?": + print("JobID: " + job["jobid"], + " Jobstatus:" + job["jobstatus"] + " Not submitting.") + sys.stdout.flush() + return + + with open(self.status_file, 'r') as f: + state = json.load(f) + + if 'feff' in state: + if state['feff'] == 'plotted': + print('XAS in ' + self.configname + ' is already plotted, skipping submit.') + return + status = state['status'] + + # check the current status + if status == "complete": + print("Preparing to submit the FEFF jobs") + sys.stdout.flush() + + # cd to configdir, submit jobs from configdir, then cd back to currdir + currdir = os.getcwd() + os.chdir(self.submit_dir) + + # construct command to be run + cmd = "" + if self.settings["prerun"] is not None: + cmd += self.settings["prerun"] + "\n" + cmd += "python -c \"from casm.feff.feff import Feff; Feff('" + self.submit_dir + "').run()\"\n" + if self.settings["postrun"] is not None: + cmd += self.settings["postrun"] + "\n" + + print("Constructing the job") + sys.stdout.flush() + + # construct the Job + job = Job(name=self.configname + '_FEFF', + account=self.settings["account"], + nodes=int(self.settings["nodes"]), + ppn=int(self.settings['ppn']), + walltime=self.settings["walltime"], + pmem=self.settings["pmem"], + queue=self.settings["queue"], + message=self.settings["message"], + email=self.settings["email"], + command=cmd) + + print("Submitting: " + self.submit_dir) + sys.stdout.flush() + # submit the job + job.submit() + + # return to current directory + os.chdir(currdir) + + print("CASM FEFF job submission complete\n") + sys.stdout.flush() + + return + else: + raise FeffError("Unexpected relaxation status: " + status + '\n' + 'To run FEFF the relaxation has to be \'complete\'') + + def plot_feff(self, plot_dir=None): + if plot_dir is None: + raise FeffError('Need to supply a diretory to plot in') + + currdir = os.getcwd() + os.chdir(plot_dir) + + st = Structure.from_file(os.path.join(self.contcar_dir, 'CONTCAR')) + + to_compute = [] + sp_compute = [] + for s in SpacegroupAnalyzer(st).get_symmetry_dataset()['equivalent_atoms']: + if to_compute.count(str(s)) == 0: + to_compute.append(str(s)) + sp_compute.append(st.species[int(s)].symbol) + + wt, ct = np.unique(SpacegroupAnalyzer(st).get_symmetry_dataset()['equivalent_atoms'], return_counts=True) + weights = dict(zip(wt, ct)) + + print('Plotting FEFF data in ' + plot_dir) + print(' Species are: ', sp_compute) + print(' Weights are: ', weights) + + num_grid = 2000 + data = dict() + + for ab in range(len(to_compute)): + edge = ['K'] + + for ed in edge: + emin = 1E31 + emax = -1E9 + + mufile = os.path.join(plot_dir, + str(sp_compute[ab]) + '_' + str(ed) + '_' + str(to_compute[ab]), 'xmu.dat') + + if not os.path.isfile(mufile): + print(' *** WARNING: File not found: %s' % mufile) + print(' You might want to rerun this for complete averages...') + continue + + w, e, k, mu0, chi, p = np.genfromtxt(mufile, unpack=True) + + for i in range(len(mu0)): + mu0 *= weights[int(to_compute[ab])] + + if self.feff_settings['feff_plot_use_omega']: + if np.amin(w) < emin: + emin = np.amin(w) + if np.amax(w) > emax: + emax = np.amax(w) + data[str(sp_compute[ab]) + str(ed) + str(to_compute[ab])] = {'data': [w, mu0], 'mins': [emin, emax]} + else: + if np.amin(e) < emin: + emin = np.amin(e) + if np.amax(e) > emax: + emax = np.amax(e) + data[str(sp_compute[ab]) + str(ed) + str(to_compute[ab])] = {'data': [e, mu0], 'mins': [emin, emax]} + + fname = os.path.join(plot_dir, str(sp_compute[ab]) + '_' + + str(ed) + '_' + str(to_compute[ab]), 'xanes_single.png') + self.plot_single(data, str(sp_compute[ab]) + str(ed) + str(to_compute[ab]), fname) + print(fname) + + # find the min / max x values to interpolate + absvals = dict() + for ab in range(len(to_compute)): + edge = ['K'] + for e in edge: + absmin = 1E31 + absmax = 0 + if not str(sp_compute[ab]) + str(e) + str(to_compute[ab]) in data: + print('Check: ', str(sp_compute[ab]) + str(e) + str(to_compute[ab])) + raise FeffError('Not all Spectra could be computed.\n If running in MPI mode' + 'try running in single core mode first before checking further.\n' + 'FEFF sometimes has convergence issues in MPI mode.') + if data[str(sp_compute[ab]) + str(e) + str(to_compute[ab])]['mins'][0] < absmin: + absmin = data[str(sp_compute[ab]) + str(e) + str(to_compute[ab])]['mins'][0] + if data[str(sp_compute[ab]) + str(e) + str(to_compute[ab])]['mins'][1] > absmax: + absmax = data[str(sp_compute[ab]) + str(e) + str(to_compute[ab])]['mins'][1] + absvals[str(sp_compute[ab]) + str(e)] = [absmin, absmax] + + species = [] + for ab in range(len(to_compute)): + if str(sp_compute[ab]) not in species: + species.append(str(sp_compute[ab])) + + # interpolate and average now + avg = np.empty((num_grid, 2, len(species))) + avg.fill(0) + + for si, sp in enumerate(species): + for ab in range(len(to_compute)): + if sp_compute[ab] != sp: + continue + edge = ['K'] + for e in edge: + if not str(sp_compute[ab]) + str(e) + str(to_compute[ab]) in data: + print('Check: ', str(sp_compute[ab]) + str(e) + str(to_compute[ab])) + continue + xi = np.linspace(absvals[str(sp_compute[ab]) + str(e)][0], + absvals[str(sp_compute[ab]) + str(e)][1], num_grid) + yi = sci_int.griddata(data[str(sp_compute[ab]) + str(e) + str(to_compute[ab])]['data'][0], + data[str(sp_compute[ab]) + str(e) + str(to_compute[ab])]['data'][1], + xi[:], method='linear', fill_value=0) + for i in range(num_grid): + avg[i, 0, si] = xi[i] + avg[i, 1, si] += yi[i] + + avg[:, 1, si] /= float(len(to_compute[ab])) + + for si, sp in enumerate(species): + edge = ['K'] + for e in edge: + fname = 'average_xanes_' + sp + '_' + str(e) + '.png' + print(fname) + lbl = r'Average ' + str(e) + '-Edge ' + sp + ' (' + str(si) + ' total)' + self.plot_avgs(avg, si, lbl, fname) + + os.chdir(currdir) + + with open(self.status_file, 'r') as f: + state = json.load(f) + + state['feff'] = 'plotted' + + with open(self.status_file, 'w') as f: + json.dump(state, f) + + def plot_single(self, in_data, loc_index, oname): + plt.rcParams['xtick.major.size'] = 8 + plt.rcParams['xtick.major.width'] = 3 + plt.rcParams['xtick.minor.size'] = 4 + plt.rcParams['xtick.minor.width'] = 3 + plt.rcParams['xtick.labelsize'] = 16 + plt.rcParams['ytick.major.size'] = 8 + plt.rcParams['ytick.major.width'] = 3 + plt.rcParams['ytick.minor.size'] = 4 + plt.rcParams['ytick.minor.width'] = 3 + plt.rcParams['ytick.labelsize'] = 16 + plt.rcParams['axes.linewidth'] = 3 + + g_y = self.gauss_broad(in_data[loc_index]['data'][0], in_data[loc_index]['data'][1], + self.fwhm2sigma(float(self.feff_settings['feff_plot_sigma']))) + scale = 1 / np.nanmax(g_y) + plt.plot(in_data[loc_index]['data'][0], g_y * scale, '-', + color='navy', label=r'O ($\sigma$ = 1 eV)', linewidth=1.0) + + plt.title(r'XANES', fontsize=16) + plt.xlabel(r'Energy [eV]', fontsize=16) + plt.ylabel(r'Absorption $\mu_{0}$', fontsize=16) + plt.legend(loc='upper right', ncol=1, fontsize=8) + plt.tight_layout() + plt.savefig(oname, dpi=170) + plt.clf() + + def plot_avgs(self, in_data, sp_num, lbls, oname): + plt.rcParams['xtick.major.size'] = 8 + plt.rcParams['xtick.major.width'] = 3 + plt.rcParams['xtick.minor.size'] = 4 + plt.rcParams['xtick.minor.width'] = 3 + plt.rcParams['xtick.labelsize'] = 16 + plt.rcParams['ytick.major.size'] = 8 + plt.rcParams['ytick.major.width'] = 3 + plt.rcParams['ytick.minor.size'] = 4 + plt.rcParams['ytick.minor.width'] = 3 + plt.rcParams['ytick.labelsize'] = 16 + plt.rcParams['axes.linewidth'] = 3 + + plt.title(r'Average XANES', fontsize=16) + plt.xlabel(r'Energy [eV]', fontsize=16) + plt.ylabel(r'Absorption $\mu_{0}$', fontsize=16) + + g_y = self.gauss_broad(in_data[:, 0, sp_num], in_data[:, 1, sp_num], + float(self.fwhm2sigma(self.feff_settings['feff_plot_sigma']))) + plt.plot(in_data[:, 0, sp_num], g_y / np.nanmax(g_y), '-', color='navy', label=lbls, linewidth=1.0) + + plt.tight_layout() + plt.savefig(oname, dpi=170) + plt.close() + + @staticmethod + def fwhm2sigma(fwhm): + return fwhm / np.sqrt(8 * np.log(2)) + + @staticmethod + def gauss_broad(x_data, y_data, sigma): + yvals = np.zeros(x_data.shape) + for ii in range(len(x_data)): + kernel = np.exp(-(x_data - x_data[ii]) ** 2 / (2 * sigma ** 2)) + kernel = kernel / sum(kernel) + yvals[ii] = sum(y_data * kernel) + return yvals + + +def check_consistent_settings(project_settings=None, feff_settings=None): + if project_settings is None and feff_settings is None: + raise FeffError('Need to specify both settings here.') + if project_settings['ppn'] != 1 or project_settings['nodes'] != 1: + if 'feff_base_mpi' not in feff_settings: + raise FeffError('You have to specify >feff_base_mpi< for FEFF ' + 'when running with more than 1 core.') + elif 'feff_base_dir' not in feff_settings: + raise FeffError('You have to specify >feff_base_dir< for FEFF.') + + +class FreezeError(object): + def __init__(self): + self.pattern = None + + def __str__(self): + return "FEFF appears to be frozen" + + @staticmethod + def error(jobdir=None): + """ Check if code is frozen + + Args: + jobdir: job directory + Returns: + True if: + 1) no file has been modified for 5 minutes + """ + + # Check if any files modified in last 300s + most_recent = None + most_recent_file = None + for f in os.listdir(jobdir): + t = time.time() - os.path.getmtime(os.path.join(jobdir, f)) + if most_recent is None: + most_recent = t + most_recent_file = f + elif t < most_recent: + most_recent = t + most_recent_file = f + + print(" -> Most recent file output (" + most_recent_file + "):", + np.round(most_recent, decimals=3), " seconds ago.") + sys.stdout.flush() + if most_recent < 1200: + return False + + @staticmethod + def fix(): + """ Fix by killing the job and resubmitting.""" + print(" Kill job and continue...") + raise FeffError('FEFF code was frozen [no output 20 Minutes], killed. FIXME if needed...') + + +def error_check(jobdir, stdoutfile): + """ Check vasp stdout for errors """ + err = dict() + + # Error to check line by line, only look for first of each type + sout = open(stdoutfile, 'r') + + # Error to check for once + possible = [FreezeError()] + for p in possible: + if p.error(jobdir=jobdir): + err[p.__class__.__name__] = p + + sout.close() + if len(err) == 0: + return None + else: + return err diff --git a/python/casm/casm/learn/__init__.py b/python/casm/casm/learn/__init__.py index 4bb4fe853..2747e972f 100644 --- a/python/casm/casm/learn/__init__.py +++ b/python/casm/casm/learn/__init__.py @@ -81,7 +81,7 @@ def create_halloffame(maxsize, rel_tol=1e-6): from casm.learn.feature_selection import fit_and_select from casm.learn.direct_selection import direct_fit -__all__ = __all__ = [ +__all__ = [ 'create_halloffame', 'EqualIndividual', 'empty_individual', diff --git a/python/casm/casm/learn/evolve.py b/python/casm/casm/learn/evolve.py index ccd1361d7..b4dc67ead 100644 --- a/python/casm/casm/learn/evolve.py +++ b/python/casm/casm/learn/evolve.py @@ -821,7 +821,6 @@ def _run(self): filename=self.evolve_params.halloffame_filename, n_halloffame=self.evolve_params.n_halloffame) - ## read or construct initial population self.toolbox.decorate("population", enforce_constraints(self.constraints)) self.pop = initialize_population(self.evolve_params.n_population, self.toolbox, @@ -1408,3 +1407,268 @@ def fit(self, X, y): + +class NeuralNetworkSelection(BaseEstimator, SelectorMixin): + + def __init__(self, + algorithm, + estimator, + hidden_layer_sizes, + cv=None, + evolve_params_kwargs=dict(), + constraints_kwargs=dict(), + alg_args=list(), + alg_kwargs=dict(), + stats=default_stats(), + verbose=True): + + self.algorithm = algorithm + self.alg_args = alg_args + self.alg_kwargs = alg_kwargs + + self.hidden_layer_sizes = hidden_layer_sizes + + self.estimator = estimator + self.cv = cv + + self.evolve_params = EvolutionaryParams(**evolve_params_kwargs) + self.constraints = Constraints(**constraints_kwargs) + + self.verbose = verbose + + self.stats = stats + + def _run(self): + """ + Run the specified evolutionary algorithm. + + Arguments + --------- + + algorithm: func, + The evolutionary algorithm to perform + + alg_args: List[func] + A list of functions used to generate positional arguments for 'algorithm' + by calling f(self). + + alg_kwargs: dict + A dict of key:function pairs used to generate keyword arguments for 'algorithm'. + + Example: + + alg_args = [ + operator.attrgetter("pop"), + operatot.attrgetter("toolbox")] + + alg_kwargs = { + "stats":operator.attrgetter("stats"), + "halloffame":operatot.attrgetter("halloffame") + } + + Then: + + in_args = [f(self) for f in alg_args] + + in_kwargs = dict() + for key, f in six.iteritems(alg_kwargs): + in_kwargs[key] = f(self) + + self.pop, self.logbook = algorithm(*in_args, **in_kwargs) + + + Returns + ------- + + self: returns an instance of self. + """ + ## read or construct hall of fame + self.halloffame = initialize_halloffame( + filename=self.evolve_params.halloffame_filename, + n_halloffame=self.evolve_params.n_halloffame) + + print('running...') + ## Run algorithm + + + quit() + + + ## Print end population + if self.verbose: + print("\nFinal population:") + casm.learn.print_population(self.pop_end) + + save_population(self.pop, filename=self.evolve_params.pop_end_filename, verbose=self.verbose) + + ## Print hall of fame + if self.verbose: + print("\nHall of Fame:") + casm.learn.print_population(self.halloffame) + + save_halloffame(self.halloffame, filename=self.evolve_params.halloffame_filename, verbose=self.verbose) + + return self + + def _get_support_mask(self): + """ + Return most fit inidividual found. + + Returns + ------- + + support: List[bool] of length n_features + This is a boolean list of shape [n_features], in which an element is True + iff its corresponding feature is selected for retention. + + """ + return self.halloffame[0] + + def get_halloffame(self): + """ + Returns the hall of fame of most fit individuals generated by 'fit'. + + Returns + ------- + + hall: deap.tools.HallOfFame + A Hall Of Fame containing the optimal solutions, as judged by CV score. + + """ + return self.halloffame + + + +class MLPRegressor(NeuralNetworkSelection): + """ + Population best first algorithm for regression by optimizing a CV score. + + Implements a best first search optimization for a population of individual + solutions. Each individual is associated with a 'status' that indicates + whether it has had children yet or not. At each step, the most fit individual + that hasn't had children has children and the population is updated to keep + only the 'n_population' most fit individuals. The algorithm stops when all + individuals in the population have had children. + + By default, children are generated by generating all the individual that differ + from the parent by +/- 1 selected feature. But any generating function may be + given that generates offspring from a single parent individual. + + + Attributes + ---------- + + estimator: estimator object implementing 'fit' + The estimator specified by the input settings. + + scoring: string + A string or a scorer callable object / function with signature + scorer(estimator, X, y). The parameter for sklearn.model_selection.cross_val_score, + default = None, uses estimator.score(). + + cv: cross-validation generator or an iterable + Provides train/test splits + + penalty: float + The CV score is increased by 'penalty*(number of selected basis function)' + + evolve_params: casm.learn.evolve.EvolutionaryParams + A EvolutionaryParams object containing parameters. + + constraints: casm.learn.evolve.Constraints + A Constraints object containing constraints on allowed individuals. + + children: func + A function with signature 'offspring = func(indiv)'. + + toolbox: deap.base.Toolbox + Contains methods used by deap during evolution. + + verbose: boolean + Print information to stdout. + + hall: deap.tools.HallOfFame + A Hall Of Fame containing the optimal solutions, as judged by CV score. + + pop_begin: iterable of individuals + The initial population of individual solutions. + + pop_end: iterable of individuals + The final population of individual solutions. + """ + + def __init__(self, estimator, hidden_layer_sizes, cv=None, + evolve_params_kwargs=dict(), constraints_kwargs=dict(), + children=single_flip_children, + verbose=True): + """ + Arguments + --------- + + estimator: estimator object implementing 'fit' + The estimator specified by the input settings. + + scoring: string, callable or None, optional, default=None + A string or a scorer callable object / function with signature + scorer(estimator, X, y). The parameter is passed to sklearn.model_selection.cross_val_score, + has a default=None which uses estimator.score(). + + cv: cross-validation generator or an iterable, optional, default=None + Provides train/test splits. The parameter is passed to sklearn.model_selection.cross_val_score, + has a default=None which uses KFold cross-validation with k=3. + + penalty: float, optional, default=0.0 + The CV score is increased by 'penalty*(number of selected basis function)' + + evolve_params_kwargs: dict, optional, default=dict() + Keyword arguments for construction of EvolutionaryParams object containing + parameters. See casm.learn.evolve.EvolutionaryParams for possibilities. + + constraints_kwargs: dict, optional, default=dict() + Keyword arguments for construction of Constraints object containing + constraints on allowed individuals. See casm.learn.evolve.Constraints for + possibilities. + + children: func, optional, default=single_flip_children + A function with signature 'offspring = func(indiv)'. + + verbose: boolean, optional, default=True + Print information to stdout. + + """ + # add a count of the number of individuals fully optimized + stats = default_stats() + stats_indiv = deap.tools.Statistics(key=lambda indiv: indiv.parent) + stats_indiv.register("parents", sum) + mstats = deap.tools.MultiStatistics(fitness=stats, optimization=stats_indiv) + + super(MLPRegressor, self).__init__( + eaPopulationBestFirst, + estimator, + hidden_layer_sizes=hidden_layer_sizes, + cv=cv, + evolve_params_kwargs=evolve_params_kwargs, + constraints_kwargs=constraints_kwargs, + verbose=verbose, + alg_args=[ + attrgetter("pop"), + attrgetter("toolbox"), + lambda obj: obj.evolve_params.n_generation], + alg_kwargs={ + "verbose": attrgetter("verbose"), + "halloffame": attrgetter("halloffame"), + "stats": attrgetter("stats") + }, + stats=mstats) + + self.hidden_layer_sizes = hidden_layer_sizes + + + def fit(self, X, y): + + clfa = MLPRegressor(estimator=self.estimator, hidden_layer_sizes=50, verbose=True) + +# a = clfa.fit(X, y) + + + return clfa.fit(X,y) diff --git a/python/casm/casm/learn/feature_selection.py b/python/casm/casm/learn/feature_selection.py index 85e47e4d4..c606f3660 100644 --- a/python/casm/casm/learn/feature_selection.py +++ b/python/casm/casm/learn/feature_selection.py @@ -5,7 +5,8 @@ from casm.learn.fit import make_fitting_data, make_estimator, make_selector, \ add_individual_detail, print_halloffame import casm.learn.model_selection -from casm.learn.evolve import GeneticAlgorithm, IndividualBestFirst, PopulationBestFirst +from casm.learn.evolve import GeneticAlgorithm, MLPRegressor, IndividualBestFirst, PopulationBestFirst + def fit_and_select(input, save=True, verbose=True, read_existing=True, hall=None): """ @@ -57,7 +58,7 @@ def fit_and_select(input, save=True, verbose=True, read_existing=True, hall=None # feature selection selector = make_selector(input, estimator, scoring=fdata.scoring, cv=fdata.cv, penalty=fdata.penalty, verbose=verbose) - + if not hasattr(selector, "get_halloffame") and not hasattr(selector, "get_support"): raise Exception("Selector has neither 'get_halloffame' nor 'get_support'") @@ -113,5 +114,3 @@ def fit_and_select(input, save=True, verbose=True, read_existing=True, hall=None hall.update([indiv]) return (fdata, estimator, selector) - - diff --git a/python/casm/casm/learn/fit.py b/python/casm/casm/learn/fit.py index 6115aba19..d3c59974d 100644 --- a/python/casm/casm/learn/fit.py +++ b/python/casm/casm/learn/fit.py @@ -9,6 +9,7 @@ import sklearn.linear_model import sklearn.model_selection +import sklearn.neural_network import sklearn.metrics import random, re, time, os, types, json, pickle, copy, uuid, shutil, tempfile import numpy as np @@ -1569,7 +1570,8 @@ def check_input(name): cv_kwargs = copy.deepcopy(specs["cv"]["kwargs"]) # get cv method (required user input) - cv_method = _find_method([sklearn.model_selection, casm.learn.model_selection], specs["cv"]["method"]) + cv_method = _find_method([sklearn.model_selection, sklearn.neural_network, + casm.learn.model_selection], specs["cv"]["method"]) cv = cv_method(**cv_kwargs) if verbose: @@ -1606,7 +1608,7 @@ def check_input(name): return fdata -def make_estimator(input, verbose = True): +def make_estimator(input, verbose=True): """ Construct estimator object from input settings. @@ -1637,14 +1639,18 @@ def make_estimator(input, verbose = True): # get kwargs (default: fit_intercept=False) kwargs = copy.deepcopy(input["estimator"]["kwargs"]) + print(kwargs) + # Use casm version of LinearRegression if input["estimator"]["method"] == "LinearRegression": estimator_method = casm.learn.linear_model.LinearRegressionForLOOCV + elif input["estimator"]["method"] == "MLPRegressor": + kwargs = copy.deepcopy(input["estimator"]["kwargs"]) + estimator_method = _find_method([sklearn.neural_network], input["estimator"]["method"]) else: + kwargs = False estimator_method = _find_method([sklearn.linear_model], input["estimator"]["method"]) - estimator = estimator_method(**kwargs) - - + estimator = estimator_method() return estimator @@ -1714,8 +1720,12 @@ def make_selector(input, estimator, scoring=None, cv=None, penalty=0.0, verbose= # if "verbose" in allowed_kwargs: # kwargs["verbose"] = verbose + print('sldiufbv') + print(selector_method) + sig = signature(selector_method.__init__) -# print("sig.parameters.keys():", str(sig.parameters.keys())) + print("sig.parameters.keys():", str(sig.parameters.keys())) + if "cv" in sig.parameters: kwargs["cv"] = cv diff --git a/python/casm/casm/project/io.py b/python/casm/casm/project/io.py index 11e6754e0..275924ca3 100644 --- a/python/casm/casm/project/io.py +++ b/python/casm/casm/project/io.py @@ -1,10 +1,16 @@ -"""casm.project file io""" -from __future__ import (absolute_import, division, print_function, unicode_literals) -from builtins import * +from __future__ import absolute_import, division, print_function, unicode_literals import json import six -from casm.misc import compat, noindent +from casm.misc import noindent + + +class ProjectIOError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg def write_eci(proj, eci, fit_details=None, clex=None, verbose=False): """ @@ -26,14 +32,16 @@ def write_eci(proj, eci, fit_details=None, clex=None, verbose=False): clex: ClexDescription instance, optional, default=proj.settings.default_clex Specifies where to write the ECI + + verbose: be verbose? """ - dir = proj.dir + wdir = proj.dir if clex is None: clex = proj.settings.default_clex # read basis.json - filename = dir.basis(clex) + filename = wdir.basis(clex) with open(filename, 'rb') as f: j = json.loads(f.read().decode('utf-8')) #print(json.dumps(j, indent=2)) @@ -47,7 +55,7 @@ def write_eci(proj, eci, fit_details=None, clex=None, verbose=False): # pretty printing for entry in j["site_functions"]: - if entry["basis"] != None: + if entry["basis"] is not None: basis = entry["basis"] for key, val in basis.items(): basis[key] = noindent.NoIndent(val) @@ -58,7 +66,7 @@ def write_eci(proj, eci, fit_details=None, clex=None, verbose=False): sites[i] = noindent.NoIndent(sites[i]) # write eci.json - filename = dir.eci(clex) + filename = wdir.eci(clex) if verbose: print("Writing:", filename, "\n") @@ -67,3 +75,223 @@ def write_eci(proj, eci, fit_details=None, clex=None, verbose=False): # refresh proj to reflect new eci proj.refresh(clear_clex=True) + + +def read_project_settings(filename): + """Returns a JSON object reading JSON files containing settings for Quantum Espresso PBS jobs. + + Returns: + settings = a JSON object containing the settings file contents + This can be accessed like a dict: settings["account"], etc. + ** All values are expected to be 'str' type. ** + + The required keys are: + "queue": queue to submit job in + "ppn": processors (cores) per node to request + "atom_per_proc": max number of atoms per processor (core) <--- is this a vasp parsing field or pbs? + "walltime": walltime to request (ex. "48:00:00") + + The optional keys are: + "account": account to submit job under (default None) + "pmem": string for requested memory (default None) + "priority": requested job priority (default "0") + "message": when to send messages about jobs (ex. "abe", default "a") + "email": where to send messages (ex. "me@fake.com", default None) + "qos": quality of service, 'qos' option (ex. "fluxoe") + "dft_cmd": DFT execution command (NO default) + "ncpus": number of cpus (cores) to run on (default $PBS_NP) + "run_limit": number of quantum espresso runs until "not_converging" (default 10) + "nrg_convergence": converged if last two runs complete and differ in energy by less + "move": files to move at the end of a run (ex. [ ".wfc"], default []) + "copy": files to copy from run to run ( default [infilename]) + "remove": files to remove at the end of a run ( default [".wfc",".igk",".save"] + "compress": files to compress at the end of a run (ex. [outfilename], default []) + "backup": files to compress to backups at the end of a run, used in conjunction with move (ex. [".wfc"]) + "encut": [START, STOP, STEP] values for converging ecutwfc to within nrg_convergence (ex. ["450", "Auto", "10"], + default ["Auto", "Auto", "10"] where "Auto" is either the largest ENMAX in all + UPFS called in SPECIES for START, or 2.0 * largest ENMAX for STOP) + "kpoints": [start, stop, step] values for converging KPOINTS to within nrg_convergence (ex. ["5", "50", "1"], + default ["5", "Auto", "1"] <---- Needs to be adjusted for grid convergence + "extra_input_files": extra input files to be copied from the settings directory, e.g., OCCUPATIONS file. + "initial" : location of infile with tags for the initial run, if desired + "final" : location of infile with tags for the final run, if desired + "err_types" : list of errors to check for. Allowed entries are "IbzkptError" + and "SubSpaceMatrixError". Default: ["SubSpaceMatrixError"] <---- STILL NEED TO IMPLEMENT + """ + try: + with open(filename, 'rb') as file: + settings = json.loads(file.read().decode('utf-8')) + except IOError as e: + print("Error reading settings file:", filename) + raise e + +# if settings['software'] != 'vasp': +# raise IOError('ONLY VASP in experimental status...') + + required = ["queue", "ppn", "walltime", "software", "run_cmd"] + + optional = ["account", "pmem", "priority", "message", "email", "qos", "npar", "ncore", "kpar", + "ncpus", "run_limit", "nrg_convergence", "prerun", "postrun", + "encut", "kpoints", "extra_input_files", "move", "copy", "remove", "compress", + "backup", "initial", "final", "strict_kpoints", "err_types", "preamble", + "infilename", "outfilename", "atom_per_proc", "nodes", "is_slab", "fix_pos", "basis"] + + for key in required: + if key not in settings: + raise ProjectIOError(key + "' missing from: '" + filename + "'") + + for key in optional: + if key not in settings: + if key.lower() in ["extra_input_files", "remove", "compress", "backup", "move", "copy"]: + settings[key] = [] + else: + settings[key] = None + + # TODO: Define a COMMON DEFAULT_XXXX_LIST for all codes and append / handle these, don't use separates for codes + + if settings['software'] == 'vasp': + from casm.vasp.io.io import DEFAULT_VASP_COPY_LIST as DEFAULT_COPY_LIST + from casm.vasp.io.io import DEFAULT_VASP_MOVE_LIST as DEFAULT_MOVE_LIST + from casm.vasp.io.io import DEFAULT_VASP_REMOVE_LIST as DEFAULT_REMOVE_LIST + DEFAULT_GZIP_LIST = [] + elif settings['software'] == 'qe': + from casm.quantumespresso.qeio import DEFAULT_QE_COPY_LIST as DEFAULT_COPY_LIST + from casm.quantumespresso.qeio import DEFAULT_QE_MOVE_LIST as DEFAULT_MOVE_LIST + from casm.quantumespresso.qeio import DEFAULT_QE_REMOVE_LIST as DEFAULT_REMOVE_LIST + DEFAULT_GZIP_LIST = [] + elif settings['software'] == 'aims': + from casm.aims.io.io import DEFAULT_AIMS_COPY_LIST as DEFAULT_COPY_LIST + from casm.aims.io.io import DEFAULT_AIMS_MOVE_LIST as DEFAULT_MOVE_LIST + from casm.aims.io.io import DEFAULT_AIMS_REMOVE_LIST as DEFAULT_REMOVE_LIST + from casm.aims.io.io import DEFAULT_AIMS_GZIP_LIST as DEFAULT_GZIP_LIST + else: + raise ProjectIOError('Software needs to be defined in relax.json.') + + if type(settings["remove"]) == list: + if 'default' in settings["remove"]: + settings["remove"] += DEFAULT_REMOVE_LIST + elif type(settings["remove"]) == str: + if settings["remove"].lower() == 'default': + settings["remove"] = DEFAULT_REMOVE_LIST + else: + settings["remove"] = [settings["remove"]] + + if "priority" not in settings: + settings["priority"] = 0 + if "extra_input_files" not in settings: + settings["extra_input_files"] = [] + if "strict_kpoints" not in settings: + settings["strict_kpoints"] = False + + for k in settings.keys(): + if k not in required: + if k not in optional: + print("unknown key '" + k + "' found in: '" + filename + "'") + raise ProjectIOError("unknown key '" + k + "' found in: '" + filename + "'") + + settings['DEF_CP'] = DEFAULT_COPY_LIST + settings['DEF_MV'] = DEFAULT_MOVE_LIST + settings['DEF_RM'] = DEFAULT_REMOVE_LIST + settings['DEF_GZ'] = DEFAULT_GZIP_LIST + + return settings + + +def read_band_settings(filename): + """Returns a JSON object reading JSON files containing settings for band structure plots. + + Returns: + settings = a JSON object containing the settings file contents + This can be accessed like a dict: settings["account"], etc. + ** All values are expected to be 'str' type. ** + """ + try: + with open(filename, 'rb') as file: + settings = json.loads(file.read().decode('utf-8')) + except IOError as e: + print("Error reading settings file:", filename) + raise e + + band_args = ["band_subdiv", "band_bs_projection", "band_dos_projection", "band_vb_energy_range", + "band_cb_energy_range", "band_fixed_cb_energy", "band_egrid_interval", "band_font", + "band_axis_fontsize", "band_tick_fontsize", "band_legend_fontsize", "band_bs_legend", + "band_dos_legend", "band_rgb_legend", "band_fig_size", "band_plot_name", "band_plot_dpi", + "band_kpt_dens"] + + for key in band_args: + if key not in settings: + if key == 'band_kpt_dens': + settings['band_kpt_dens'] = 1000 + if key == 'band_subdiv': + settings['band_subdiv'] = 100 + if key == 'band_bs_projection': + settings['band_bs_projection'] = 'elements' + if key == 'band_dos_projection': + settings['band_dos_projection'] = 'elements' + if key == 'band_vb_energy_range': + settings['band_vb_energy_range'] = 4 + if key == 'band_cb_energy_range': + settings['band_cb_energy_range'] = 4 + if key == 'band_fixed_cb_energy': + settings['band_fixed_cb_energy'] = False + if key == 'band_egrid_interval': + settings['band_egrid_interval'] = 1 + if key == 'band_font': + settings['band_font'] = 'Times New Roman' + if key == 'band_axis_fontsize': + settings['band_axis_fontsize'] = 20 + if key == 'band_tick_fontsize': + settings['band_tick_fontsize'] = 15 + if key == 'band_legend_fontsize': + settings['band_legend_fontsize'] = 14 + if key == 'band_bs_legend': + settings['band_bs_legend'] = 'best' + if key == 'band_dos_legend': + settings['band_dos_legend'] = 'best' + if key == 'band_rgb_legend': + settings['band_rgb_legend'] = True + if key == 'band_fig_size': + settings['band_fig_size'] = (11, 8.5) + if key == 'band_plot_name': + settings['band_plot_name'] = 'BSDOS.png' + if key == 'band_plot_dpi': + settings['band_plot_dpi'] = 170 + + return settings + + +def read_feff_settings(filename): + """Returns a JSON object reading JSON files containing settings for XAS computations. + + Returns: + settings = a JSON object containing the settings file contents + This can be accessed like a dict: settings["account"], etc. + ** All values are expected to be 'str' type. ** + """ + try: + with open(filename, 'rb') as file: + settings = json.loads(file.read().decode('utf-8')) + except IOError as e: + print("Error reading settings file:", filename) + raise e + + feff_args = ["feff_cmd", "feff_nkpts", "feff_radius", "feff_user_tags", + "feff_plot_sigma", "feff_plot_use_omega", "feff_base_dir", + "feff_base_mpi"] + + for key in feff_args: + if key not in settings: + if key == 'feff_cmd': + settings['feff_cmd'] = 'feff' + if key == 'feff_nkpts': + settings['feff_nkpts'] = 1000 + if key == 'feff_radius': + settings['feff_radius'] = 12 + if key == 'feff_user_tags': + settings['feff_user_tags'] = {} + if key == 'feff_plot_sigma': + settings['feff_plot_sigma'] = 1.0 + if key == 'feff_plot_use_omega': + settings['feff_plot_use_omega'] = False + + return settings diff --git a/python/casm/casm/scripts/casm_bands.py b/python/casm/casm/scripts/casm_bands.py new file mode 100644 index 000000000..b22160c8a --- /dev/null +++ b/python/casm/casm/scripts/casm_bands.py @@ -0,0 +1,140 @@ +from __future__ import absolute_import, division, print_function, unicode_literals + +import os +import sys +import argparse + +from casm.project.io import read_project_settings +from casm.project import Project, Selection +from casm.vasp.bands import Bands as VaspBand + + +class CasmBandsError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +""" +casm-bands --run / --submit / --setup / --report +""" + +configs_help = """ +CASM selection file or one of 'CALCULATED', 'ALL', or 'MASTER' (Default) +""" + +run_help = """ +Run calculation for all selected configurations. +""" + +submit_help = """ +Submit calculation for all selected configurations. +""" + +setup_help = """ +Setup calculation for all selected configurations. +""" + +plot_help = """ +Plots band structure and DoS for selected configurations. +""" + + +def main(argv=None): + if argv is None: + argv = sys.argv[1:] + + parser = argparse.ArgumentParser(description='Submit calculations for CASM') + parser.add_argument('-c', '--configs', help=configs_help, type=str, default="MASTER") + parser.add_argument('--run', help=run_help, action="store_true", default=False) + parser.add_argument('--submit', help=submit_help, action="store_true", default=False) + parser.add_argument('--setup', help=setup_help, action="store_true", default=False) + parser.add_argument('--plot', help=plot_help, action="store_true", default=False) + args = parser.parse_args(argv) + + args.path = os.getcwd() + + proj = Project(os.path.abspath(args.path)) + sel = Selection(proj, args.configs, all=False) + if sel.data["configname"] is not None: + configname = sel.data["configname"][0] + casm_settings = proj.settings + if casm_settings is None: + raise CasmBandsError("Not in a CASM project. The file '.casm' directory was not found.") + else: + raise CasmBandsError("Not in CASM project.") + + casm_directories = proj.dir + print(" Reading relax.json settings file") + sys.stdout.flush() + + setfile = casm_directories.settings_path_crawl("relax.json", configname, casm_settings.default_clex) + if setfile is None: + raise CasmBandsError("Could not find \"relax.json\" in \"settings\" directory") + else: + print("Using " + str(setfile) + " as computation settings...") + + bandfile = casm_directories.settings_path_crawl("bands.json", configname, casm_settings.default_clex) + if bandfile is None: + raise CasmBandsError("Could not find \"bands.json\" in \"settings\" directory") + else: + print("Using " + str(bandfile) + " as band structure settings...") + + settings = read_project_settings(setfile) + print("DFT software is:", settings['software']) + + if settings['software'] != 'vasp': + raise CasmBandsError('This is currently ONLY VASP capable.') + + band_calculator = None + + if args.setup: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmBandsError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmBandsError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + band_calculator = VaspBand(proj.dir.configuration_dir(configname)) + band_calculator.setup_chg() + + elif args.submit: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmBandsError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmBandsError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + band_calculator = VaspBand(proj.dir.configuration_dir(configname)) + band_calculator.submit() + + elif args.run: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmBandsError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmBandsError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + band_calculator = VaspBand(proj.dir.configuration_dir(configname)) + band_calculator.run() + + elif args.plot: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmBandsError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmBandsError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + band_calculator = VaspBand(proj.dir.configuration_dir(configname)) + band_calculator.plot_bandos(plot_dir=os.path.abspath(os.path.join(proj.dir.configuration_dir(configname), + 'calctype.default', 'band_structure'))) + + +if __name__ == "__main__": + main() diff --git a/python/casm/casm/scripts/casm_calc.py b/python/casm/casm/scripts/casm_calc.py index 083da6263..6cde6ed9f 100755 --- a/python/casm/casm/scripts/casm_calc.py +++ b/python/casm/casm/scripts/casm_calc.py @@ -1,22 +1,32 @@ -from __future__ import (absolute_import, division, print_function, unicode_literals) -from builtins import * +from __future__ import absolute_import, division, print_function, unicode_literals -import argparse -import json -from os import getcwd -from os.path import join, abspath +import os import sys -import six +import json +import argparse +import casm.project.io from casm.misc import compat, noindent from casm.project import Project, Selection -import casm.qewrapper -from casm.vaspwrapper import Relax +from casm.qewrapper import Relax as QERelax +from casm.vaspwrapper import Relax as VaspRelax +from casm.aimswrapper.relax import Relax as AimsRelax + + +class CasmCalcError(Exception): + def __init__(self, msg): + self.msg = msg -# casm-calc --configs selection -# --software "quantumespresso" "vasp" -# --scheduler "pbs" -# --run / --submit / --setup / --report + def __str__(self): + return self.msg + + +""" +casm-calc --configs selection + --software "quantumespresso" "vasp" "aims" + --scheduler "pbs" + --run / --submit / --setup / --report +""" configs_help = """ CASM selection file or one of 'CALCULATED', 'ALL', or 'MASTER' (Default) @@ -46,98 +56,110 @@ Report calculation results (print calc.properties.json file) for all selected configurations. """ -def main(argv = None): - if argv is None: - argv = sys.argv[1:] + +def main(argv=None): + if argv is None: + argv = sys.argv[1:] - parser = argparse.ArgumentParser(description = 'Submit calculations for CASM') - parser.add_argument('-c', '--configs', help=configs_help, type=str, default="MASTER") - parser.add_argument('--path', help=path_help, type=str, default=None) - parser.add_argument('-m','--method', help=method_help, type=str, default=None) - parser.add_argument('--run', help=run_help, action="store_true", default=False) - parser.add_argument('--submit', help=submit_help, action="store_true", default=False) - parser.add_argument('--setup', help=setup_help, action="store_true", default=False) - parser.add_argument('--report', help=report_help, action="store_true", default=False) - args = parser.parse_args(argv) + parser = argparse.ArgumentParser(description='Submit calculations for CASM') + parser.add_argument('-c', '--configs', help=configs_help, type=str, default="MASTER") + parser.add_argument('--path', help=path_help, type=str, default=None) + parser.add_argument('-m', '--method', help=method_help, type=str, default=None) + parser.add_argument('--run', help=run_help, action="store_true", default=False) + parser.add_argument('--submit', help=submit_help, action="store_true", default=False) + parser.add_argument('--setup', help=setup_help, action="store_true", default=False) + parser.add_argument('--report', help=report_help, action="store_true", default=False) + args = parser.parse_args(argv) - if args.path is None: - args.path = getcwd() + if args.path is None: + args.path = os.getcwd() - try: - proj = Project(abspath(args.path)) - sel = Selection(proj, args.configs, all=False) - if sel.data["configname"] is not None: - configname=sel.data["configname"][0] - casm_settings=proj.settings - if casm_settings == None: - raise casm.qewrapper.QEWrapperError("Not in a CASM project. The file '.casm' directory was not found.") - casm_directories=proj.dir - print(" Reading relax.json settings file") - sys.stdout.flush() - setfile = casm_directories.settings_path_crawl("relax.json",configname,casm_settings.default_clex) - if setfile == None: - raise casm.qewrapper.QEWrapperError("Could not find \"relax.json\" in an appropriate \"settings\" directory") - sys.stdout.flush() - else: - print("Using "+str(setfile)+" as settings...") - settings = casm.qewrapper.read_settings(setfile) - if settings["software"] is None: - settings["software"]="vasp" - software=settings["software"] - print("Relevant software is:", software) - if args.setup: - sel.write_pos() - for configname in sel.data["configname"]: - if software == "quantumespresso": - relaxation = casm.qewrapper.Relax(proj.dir.configuration_dir(configname)) + try: + proj = Project(os.path.abspath(args.path)) + sel = Selection(proj, args.configs, all=False) + if sel.data["configname"] is not None: + configname = sel.data["configname"][0] + casm_settings = proj.settings + if casm_settings is None: + raise CasmCalcError("Not in a CASM project. The file '.casm' directory was not found.") else: - relaxation = Relax(proj.dir.configuration_dir(configname)) - relaxation.setup() - - elif args.submit: - sel.write_pos() - for configname in sel.data["configname"]: - if software == "quantumespresso": - relaxation = casm.qewrapper.Relax(proj.dir.configuration_dir(configname)) + raise CasmCalcError("Not in CASM project.") + + casm_directories = proj.dir + print(" Reading relax.json settings file") + sys.stdout.flush() + setfile = casm_directories.settings_path_crawl("relax.json", configname, casm_settings.default_clex) + + if setfile is None: + raise CasmCalcError("Could not find \"relax.json\" in \"settings\" directory") else: - relaxation = Relax(proj.dir.configuration_dir(configname)) - relaxation.submit() + print("Using " + str(setfile) + " as settings...") + + settings = casm.project.io.read_project_settings(setfile) + print("DFT software is:", settings['software']) + + relaxation = None + + if args.setup: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + relaxation = QERelax(proj.dir.configuration_dir(configname)) + elif settings['software'] == "aims": + relaxation = AimsRelax(proj.dir.configuration_dir(configname)) + elif settings['software'] == 'vasp': + relaxation = VaspRelax(proj.dir.configuration_dir(configname)) + relaxation.setup() - elif args.run: - sel.write_pos() - for configname in sel.data["configname"]: - if software == "quantumespresso": - relaxation = casm.qewrapper.Relax(proj.dir.configuration_dir(configname)) - else: - relaxation = Relax(proj.dir.configuration_dir(configname)) - relaxation.run() + elif args.submit: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + relaxation = QERelax(proj.dir.configuration_dir(configname)) + elif settings['software'] == "aims": + relaxation = AimsRelax(proj.dir.configuration_dir(configname)) + elif settings['software'] == 'vasp': + relaxation = VaspRelax(proj.dir.configuration_dir(configname)) + relaxation.submit() + + elif args.run: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + relaxation = QERelax(proj.dir.configuration_dir(configname)) + elif settings['software'] == "aims": + relaxation = AimsRelax(proj.dir.configuration_dir(configname)) + elif settings['software'] == 'vasp': + relaxation = VaspRelax(proj.dir.configuration_dir(configname)) + relaxation.run() - elif args.report: - for configname in sel.data["configname"]: - configdir = proj.dir.configuration_dir(configname) - clex = proj.settings.default_clex - calcdir = proj.dir.calctype_dir(configname, clex) - finaldir = join(calcdir, "run.final") - try: - if software == "quantumespresso": - if settings["outfilename"] is None: - print("WARNING: No output file specified in relax.json using default outfilename of std.out") - settings["outfilename"]="std.out" - outfilename = settings["outfilename"] - output = casm.qewrapper.Relax.properties(finaldir,outfilename) - else: - output = Relax.properties(finaldir) - calc_props = proj.dir.calculated_properties(configname, clex) - print("writing:", calc_props) - with open(calc_props, 'wb') as f: - f.write(six.u(json.dumps(output, cls=noindent.NoIndentEncoder, indent=2)).encode('utf-8')) - #compat.dump(json, output, calc_props, 'w', cls=noindent.NoIndentEncoder, indent=4, sort_keys=True) - except: - print(("Unable to report properties for directory {}.\n" - "Please verify that it contains a completed calculation.".format(configdir))) - except Exception as e: - print(e) - sys.exit(1) + elif args.report: + for configname in sel.data["configname"]: + configdir = proj.dir.configuration_dir(configname) + clex = proj.settings.default_clex + calcdir = proj.dir.calctype_dir(configname, clex) + finaldir = os.path.join(calcdir, "run.final") + try: + if settings['software'] == "quantumespresso": + if settings["outfilename"] is None: + print("WARNING: No output file specified in relax.json using std.out") + settings["outfilename"] = "std.out" + outfilename = settings["outfilename"] + output = casm.qewrapper.Relax.properties(finaldir, outfilename) + else: + output = VaspRelax.properties(finaldir) + except IOError: + print(("Unable to report properties for directory {}.\n" + "Please verify that it contains a completed calculation.".format(configdir))) + + calc_props = proj.dir.calculated_properties(configname, clex) + print("writing:", calc_props) + # TODO: Fix this line to work properly + compat.dump(json, output, calc_props, 'w', cls=noindent.NoIndentEncoder, indent=4, sort_keys=True) + + except ValueError as e: + raise CasmCalcError(e) + if __name__ == "__main__": - main() + main() diff --git a/python/casm/casm/scripts/casm_feff.py b/python/casm/casm/scripts/casm_feff.py new file mode 100644 index 000000000..5f73b12c8 --- /dev/null +++ b/python/casm/casm/scripts/casm_feff.py @@ -0,0 +1,143 @@ +from __future__ import absolute_import, division, print_function, unicode_literals + +import os +import sys +import argparse + +from casm.project.io import read_project_settings, read_feff_settings +from casm.project import Project, Selection +from casm.feff import Feff, check_consistent_settings + + +class CasmFeffError(Exception): + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +""" +casm-feff --run / --submit / --setup / --plot +""" + +configs_help = """ +CASM selection file or one of 'CALCULATED', 'ALL', or 'MASTER' (Default) +""" + +run_help = """ +Run calculation for all selected configurations. +""" + +submit_help = """ +Submit calculation for all selected configurations. +""" + +setup_help = """ +Setup calculation for all selected configurations. +""" + +plot_help = """ +Plots XANES for selected configurations. +""" + + +def main(argv=None): + if argv is None: + argv = sys.argv[1:] + + parser = argparse.ArgumentParser(description='Submit calculations for CASM') + parser.add_argument('-c', '--configs', help=configs_help, type=str, default="MASTER") + parser.add_argument('--run', help=run_help, action="store_true", default=False) + parser.add_argument('--submit', help=submit_help, action="store_true", default=False) + parser.add_argument('--setup', help=setup_help, action="store_true", default=False) + parser.add_argument('--plot', help=plot_help, action="store_true", default=False) + args = parser.parse_args(argv) + + args.path = os.getcwd() + + proj = Project(os.path.abspath(args.path)) + sel = Selection(proj, args.configs, all=False) + if sel.data["configname"] is not None: + configname = sel.data["configname"][0] + casm_settings = proj.settings + if casm_settings is None: + raise CasmFeffError("Not in a CASM project. The file '.casm' directory was not found.") + else: + raise CasmFeffError("Not in CASM project.") + + casm_directories = proj.dir + print(" Reading relax.json settings file") + sys.stdout.flush() + + setfile = casm_directories.settings_path_crawl("relax.json", configname, casm_settings.default_clex) + if setfile is None: + raise CasmFeffError("Could not find \"relax.json\" in \"settings\" directory") + else: + print("Using " + str(setfile) + " as computation settings...") + settings = read_project_settings(setfile) + + fefffile = casm_directories.settings_path_crawl("feff.json", configname, casm_settings.default_clex) + if fefffile is None: + raise CasmFeffError("Could not find \"feff.json\" in \"settings\" directory") + else: + print("Using " + str(fefffile) + " for FEFF settings...") + feff_settings = read_feff_settings(fefffile) + + print("DFT software is:", settings['software']) + + if settings['software'] != 'vasp': + raise CasmFeffError('This is currently ONLY VASP capable.') + + check_consistent_settings(settings, feff_settings) + + feff_calculator = None + + if args.setup: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmFeffError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmFeffError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + feff_calculator = Feff(proj.dir.configuration_dir(configname)) + feff_calculator.setup() + + elif args.submit: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmFeffError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmFeffError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + feff_calculator = Feff(proj.dir.configuration_dir(configname)) + feff_calculator.submit() + + elif args.run: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmFeffError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmFeffError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + feff_calculator = Feff(proj.dir.configuration_dir(configname)) + feff_calculator.run() + + elif args.plot: + sel.write_pos() + for configname in sel.data["configname"]: + if settings['software'] == "quantumespresso": + raise CasmFeffError('QE not implemented, use VASP') + elif settings['software'] == "aims": + raise CasmFeffError('FHI-aims not implemented, use VASP') + elif settings['software'] == 'vasp': + feff_calculator = Feff(proj.dir.configuration_dir(configname)) + feff_calculator.plot_feff(plot_dir=os.path.abspath(os.path.join(proj.dir.configuration_dir(configname), + 'calctype.default', 'xanes'))) + + +if __name__ == "__main__": + main() diff --git a/python/casm/casm/vasp/bands.py b/python/casm/casm/vasp/bands.py new file mode 100644 index 000000000..c97eb28b1 --- /dev/null +++ b/python/casm/casm/vasp/bands.py @@ -0,0 +1,526 @@ +import re +import os +import sys +import time +import json +import matplotlib +import subprocess + +import numpy as np +import shutil as sh + +from casm.project import DirectoryStructure, ProjectSettings +from casm.project.io import read_project_settings, read_band_settings + +from pymatgen.io.vasp import Kpoints, Procar +from pymatgen.io.vasp.outputs import Vasprun +from pymatgen.core.structure import Structure +from pymatgen.symmetry.bandstructure import HighSymmKpath +from pymatgen.electronic_structure.plotter import BSDOSPlotter + +from prisms_jobs import Job, JobDB + +matplotlib.use('Agg') + + +class BandsError: + def __init__(self, msg): + self.msg = msg + + def __str__(self): + return self.msg + + +class Bands(object): + """ + Computes band structure with pymatgen standard settings + - k-point density is 1000 / atom + - high symmetry path is determined from cell geometry + """ + + def __init__(self, rundir=None): + if rundir is None: + raise BandsError('Can not create from nothing-directory') + + print("Constructing a VASP Band object") + sys.stdout.flush() + + print(" Reading CASM settings") + self.casm_directories = DirectoryStructure(rundir) + self.casm_settings = ProjectSettings(rundir) + if self.casm_settings is None: + raise BandsError("Not in a CASM project. The '.casm' directory was not found.") + + # fixed to default_clex for now + self.clex = self.casm_settings.default_clex + + _res = os.path.split(rundir) + self.configname = os.path.split(_res[0])[1] + "/" + _res[1] + + print(" Reading DFT and plot settings for configuration: " + self.configname) + sys.stdout.flush() + + setfile = self.casm_directories.settings_path_crawl("relax.json", self.configname, self.clex) + if setfile is None: + raise BandsError("Could not find relax.json in settings directory") + else: + print(" Read DFT settings from: " + setfile) + self.settings = read_project_settings(setfile) + + bandfile = self.casm_directories.settings_path_crawl("bands.json", self.configname, self.clex) + if bandfile is None: + raise BandsError("Could not find bands.json in settings directory") + else: + print(" Read band settings from: " + bandfile) + self.band_settings = read_band_settings(bandfile) + + self.submit_dir = rundir + self.band_dir = os.path.abspath(os.path.join(rundir, 'calctype.default', 'band_structure')) + self.contcar_dir = os.path.abspath(os.path.join(rundir, 'calctype.default', 'run.final')) + self.new_incar = [] + self.status_file = os.path.abspath(os.path.join(rundir, 'calctype.default', 'status.json')) + + print(" Run directory: %s" % self.band_dir) + sys.stdout.flush() + + def setup_chg(self): + """ Create VASP input files and take CONTCAR from last converged run """ + + print('Getting VASP input files for charge computation in directory: ') + print(' %s' % self.band_dir) + + if not os.path.isdir(os.path.join(self.band_dir)): + os.mkdir(os.path.join(self.band_dir)) + + s = Structure.from_file(os.path.join(self.contcar_dir, 'CONTCAR')) + + s.to(filename=os.path.join(self.band_dir, 'POSCAR'), fmt='POSCAR') + Kpoints.automatic_density(s, int(self.band_settings['band_kpt_dens']), + force_gamma=True).write_file(os.path.join(self.band_dir, 'KPOINTS')) + sh.copyfile(os.path.join(self.contcar_dir, 'POTCAR'), os.path.join(self.band_dir, 'POTCAR')) + self.manage_tags_chg(os.path.join(self.contcar_dir, 'INCAR')) + with open(os.path.join(self.band_dir, 'INCAR'), 'w') as f: + for line in self.new_incar: + f.write(line) + + def setup_band(self): + """ Create VASP input files and take CONTCAR from last converged run """ + + print('Getting VASP input files for band structure computation in directory: ') + print(' %s' % self.band_dir) + + if not os.path.isdir(os.path.join(self.band_dir)): + os.mkdir(os.path.join(self.band_dir)) + + s = Structure.from_file(os.path.join(self.contcar_dir, 'CONTCAR')) + irr_bri_zone = HighSymmKpath(s) + + s.to(filename=os.path.join(self.band_dir, 'POSCAR'), fmt='POSCAR') + Kpoints.automatic_linemode(self.band_settings['band_subdiv'], + irr_bri_zone).write_file(os.path.join(self.band_dir, 'KPOINTS')) + sh.copyfile(os.path.join(self.contcar_dir, 'POTCAR'), os.path.join(self.band_dir, 'POTCAR')) + self.manage_tags_bands(os.path.join(self.contcar_dir, 'INCAR')) + with open(os.path.join(self.band_dir, 'INCAR'), 'w') as f: + for line in self.new_incar: + f.write(line) + + def exec_chg(self, jobdir=None, stdout="chg.out", stderr="chg.err", command=None, ncpus=None, + poll_check_time=5.0, err_check_time=60.0): + """ Run selected DFT software using subprocess. + + The 'command' is executed in the directory 'jobdir'. + + Args: + jobdir: directory to run. If jobdir is None, the current directory is used. + stdout: filename to write to. If stdout is None, "std.out" is used. + stderr: filename to write to. If stderr is None, "std.err" is used. + command: (str or None) FHI-aims execution command + If command != None: then 'command' is run in a subprocess + Else, if ncpus == 1, then command = "aims" + Else, command = "mpirun -np {NCPUS} aims" + ncpus: (int) if '{NCPUS}' is in 'command' string, then 'ncpus' is substituted in the command. + if ncpus==None, $PBS_NP is used if it exists, else 1 + poll_check_time: how frequently to check if the vasp job is completed + err_check_time: how frequently to parse vasp output to check for errors + """ + print("Begin charge density run:") + sys.stdout.flush() + + if jobdir is None: + jobdir = os.getcwd() + + currdir = os.getcwd() + os.chdir(jobdir) + + if ncpus is None: + if "PBS_NP" in os.environ: + ncpus = os.environ["PBS_NP"] + elif "SLURM_NPROCS" in os.environ: + ncpus = os.environ["SLURM_NPROCS"] + else: + ncpus = 1 + + if command is None: + if ncpus == 1: + command = self.settings['run_cmd'] + else: + command = "mpirun -np {NCPUS} " + self.settings['run_cmd'] + + if re.search("NCPUS", command): + command = command.format(NCPUS=str(ncpus)) + + print(" jobdir:", jobdir) + print(" exec:", command) + sys.stdout.flush() + + err = None + sout = open(os.path.join(jobdir, stdout), 'w') + serr = open(os.path.join(jobdir, stderr), 'w') + + p = subprocess.Popen(command.split(), stdout=sout, stderr=serr) + + # wait for process to end, and periodically check for errors + last_check = time.time() + while p.poll() is None: + time.sleep(poll_check_time) + if time.time() - last_check > err_check_time: + last_check = time.time() + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + # FreezeErrors are fatal and usually not helped with abort_scf + if "FreezeError" in err.keys(): + print(" DFT for CHG run seems frozen, killing job") + sys.stdout.flush() + p.kill() + + # close output files + sout.close() + serr.close() + + os.chdir(currdir) + + print("Run ended") + sys.stdout.flush() + + # check finished job for errors + if err is None: + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + print(" Found errors:", end='') + for e in err: + print(e, end='') + print("\n") + + return err + + def exec_band(self, jobdir=None, stdout="band.out", stderr="band.err", command=None, ncpus=None, + poll_check_time=5.0, err_check_time=60.0): + """ Run selected DFT software using subprocess. + + The 'command' is executed in the directory 'jobdir'. + + Args: + jobdir: directory to run. If jobdir is None, the current directory is used. + stdout: filename to write to. If stdout is None, "std.out" is used. + stderr: filename to write to. If stderr is None, "std.err" is used. + command: (str or None) FHI-aims execution command + If command != None: then 'command' is run in a subprocess + Else, if ncpus == 1, then command = "aims" + Else, command = "mpirun -np {NCPUS} aims" + ncpus: (int) if '{NCPUS}' is in 'command' string, then 'ncpus' is substituted in the command. + if ncpus==None, $PBS_NP is used if it exists, else 1 + poll_check_time: how frequently to check if the vasp job is completed + err_check_time: how frequently to parse vasp output to check for errors + """ + print("Begin band structure run:") + sys.stdout.flush() + + if jobdir is None: + jobdir = os.getcwd() + + currdir = os.getcwd() + os.chdir(jobdir) + + if ncpus is None: + if "PBS_NP" in os.environ: + ncpus = os.environ["PBS_NP"] + elif "SLURM_NPROCS" in os.environ: + ncpus = os.environ["SLURM_NPROCS"] + else: + ncpus = 1 + + if command is None: + if ncpus == 1: + command = self.settings['run_cmd'] + else: + command = "mpirun -np {NCPUS} " + self.settings['run_cmd'] + + if re.search("NCPUS", command): + command = command.format(NCPUS=str(ncpus)) + + print(" jobdir:", jobdir) + print(" exec:", command) + sys.stdout.flush() + + err = None + sout = open(os.path.join(jobdir, stdout), 'w') + serr = open(os.path.join(jobdir, stderr), 'w') + + p = subprocess.Popen(command.split(), stdout=sout, stderr=serr) + + # wait for process to end, and periodically check for errors + last_check = time.time() + while p.poll() is None: + time.sleep(poll_check_time) + if time.time() - last_check > err_check_time: + last_check = time.time() + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + # FreezeErrors are fatal and usually not helped with abort_scf + if "FreezeError" in err.keys(): + print(" DFT run seems frozen, killing job") + sys.stdout.flush() + p.kill() + + # close output files + sout.close() + serr.close() + + os.chdir(currdir) + + print("Run ended") + sys.stdout.flush() + + # check finished job for errors + if err is None: + err = error_check(jobdir, os.path.join(jobdir, stdout)) + if err is not None: + print(" Found errors:", end='') + for e in err: + print(e, end='') + print("\n") + + return err + + def run(self): + self.setup_chg() + result = self.exec_chg(jobdir=self.band_dir) + if result is not None: + raise BandsError('Self consistent charge computation did not complete, check what happended.') + self.setup_band() + result = self.exec_band(jobdir=self.band_dir) + if result is None: + self.plot_bandos(plot_dir=self.band_dir) + else: + raise BandsError('Band structure computation did not complete, check what happended.') + + def submit(self): + """Submit a job for this band structure computation""" + print("Submitting configuration: " + self.configname) + # first, check if the job has already been submitted and is not completed + db = JobDB() + print("Calculation directory: ", self.band_dir) + sub_id = db.select_regex_id("rundir", self.band_dir) + + if sub_id is not []: + for j in sub_id: + job = db.select_job(j) + if job["jobstatus"] != "?": + print("JobID: " + job["jobid"], + " Jobstatus:" + job["jobstatus"] + " Not submitting.") + sys.stdout.flush() + return + + with open(self.status_file, 'r') as f: + state = json.load(f) + + if 'bands' in state: + if state['bands'] == 'plotted': + print(' ******** Band Structure in ' + self.configname + ' is already plotted, skipping submit.') + return + status = state['status'] + + # check the current status + if status == "complete": + print("Preparing to submit the band structure job") + sys.stdout.flush() + + # cd to configdir, submit jobs from configdir, then cd back to currdir + currdir = os.getcwd() + os.chdir(self.submit_dir) + + # construct command to be run + cmd = "" + if self.settings["prerun"] is not None: + cmd += self.settings["prerun"] + "\n" + cmd += "python -c \"from casm.vasp.bands import Bands; Bands('" + self.submit_dir + "').run()\"\n" + if self.settings["postrun"] is not None: + cmd += self.settings["postrun"] + "\n" + + print("Constructing the job") + sys.stdout.flush() + + # construct a pbs.Job + job = Job(name=self.configname + '_BANDS', + account=self.settings["account"], + nodes=int(self.settings["nodes"]), + ppn=int(self.settings['ppn']), + walltime=self.settings["walltime"], + pmem=self.settings["pmem"], + qos=self.settings["qos"], + queue=self.settings["queue"], + message=self.settings["message"], + email=self.settings["email"], + priority=self.settings["priority"], + command=cmd) + + print("Submitting") + sys.stdout.flush() + # submit the job + job.submit() + + # return to current directory + os.chdir(currdir) + + print("CASM band structure job submission complete\n") + sys.stdout.flush() + + return + + elif status == "not_converging": + print("Status:" + status + " Not submitting.") + sys.stdout.flush() + return + + elif status != "incomplete": + raise BandsError("unexpected relaxation status: '" + status) + + def manage_tags_chg(self, incar_file): + remove_tags = ['NSW', 'EDIFFG', 'IBRION', 'ISIF', 'ISMEAR', 'LCHARG'] + with open(os.path.abspath(incar_file)) as f: + for line in f: + if not any(tag in line for tag in remove_tags): + self.new_incar.append(line) + nbands = int(Procar(os.path.join(self.contcar_dir, 'PROCAR')).nbands) + self.new_incar.append('ISMEAR = 2\n') + self.new_incar.append('NBANDS = %i\n' % int(nbands * 1.5)) # use 50% more bands just to make sure + self.new_incar.append('LCHARG = .TRUE.\n') + + def manage_tags_bands(self, incar_file): + remove_tags = ['NSW', 'EDIFFG', 'IBRION', 'ISIF', 'ISMEAR', 'LCHARG'] + with open(os.path.abspath(incar_file)) as f: + for line in f: + if not any(tag in line for tag in remove_tags): + self.new_incar.append(line) + nbands = int(Procar(os.path.join(self.contcar_dir, 'PROCAR')).nbands) + self.new_incar.append('ISMEAR = 2\n') + self.new_incar.append('NBANDS = %i\n' % int(nbands * 1.5)) # use 80% more bands just to make sure + self.new_incar.append('ICHARG = 11\n') + self.new_incar.append('NEDOS = 5001\n') + self.new_incar.append('EMIN = -15\n') + self.new_incar.append('EMAX = 15\n') + self.new_incar.append('ICORELEVEL = 1') + + def plot_bandos(self, plot_dir=None): + if plot_dir is None: + raise BandsError('Need to supply a diretory to plot in') + + currdir = os.getcwd() + os.chdir(plot_dir) + + run = Vasprun('vasprun.xml') + + if not run.converged: + raise BandsError('The band structure computation of VASP did not converge for ' + self.configname + '.\n' + 'Not plotting bands / DoS from it. Check what happend and rerun...') + + dos = run.complete_dos + bst = run.get_band_structure(kpoints_filename='KPOINTS', efermi=run.efermi, line_mode=True) + + bsplot = BSDOSPlotter(bs_projection=self.band_settings['band_bs_projection'], + dos_projection=self.band_settings['band_dos_projection'], + vb_energy_range=self.band_settings['band_vb_energy_range'], + cb_energy_range=self.band_settings['band_cb_energy_range'], + fixed_cb_energy=self.band_settings['band_fixed_cb_energy'], + egrid_interval=self.band_settings['band_egrid_interval'], + font=self.band_settings['band_font'], + axis_fontsize=self.band_settings['band_axis_fontsize'], + tick_fontsize=self.band_settings['band_tick_fontsize'], + legend_fontsize=self.band_settings['band_legend_fontsize'], + bs_legend=self.band_settings['band_bs_legend'], + dos_legend=self.band_settings['band_dos_legend'], + rgb_legend=self.band_settings['band_rgb_legend'], + fig_size=self.band_settings['band_fig_size']).get_plot(bst, dos) + + bsplot.savefig(self.band_settings['band_plot_name'], dpi=self.band_settings['band_plot_dpi']) + + os.chdir(currdir) + + state = {'status': 'complete', 'bands': 'plotted'} + + with open(self.status_file, 'w') as f: + json.dump(state, f) + + +class FreezeError(object): + """DFT code appears frozen""" + def __init__(self): + self.pattern = None + + def __str__(self): + return "DFT code appears to be frozen" + + @staticmethod + def error(jobdir=None): + """ Check if code is frozen + + Args: + jobdir: job directory + Returns: + True if: + 1) no file has been modified for 5 minutes + """ + + # Check if any files modified in last 300s + most_recent = None + most_recent_file = None + for f in os.listdir(jobdir): + t = time.time() - os.path.getmtime(os.path.join(jobdir, f)) + if most_recent is None: + most_recent = t + most_recent_file = f + elif t < most_recent: + most_recent = t + most_recent_file = f + + print("Most recent file output (" + most_recent_file + "):", + np.round(most_recent, decimals=3), " seconds ago.") + sys.stdout.flush() + if most_recent < 300: + return False + + @staticmethod + def fix(): + """ Fix by killing the job and resubmitting.""" + print(" Kill job and continue...") + raise BandsError('DFT code was frozen [no outpot 5 Minutes], killed. FIXME if needed...') + + +def error_check(jobdir, stdoutfile): + """ Check vasp stdout for errors """ + err = dict() + + # Error to check line by line, only look for first of each type + sout = open(stdoutfile, 'r') + + # Error to check for once + possible = [FreezeError()] + for p in possible: + if p.error(jobdir=jobdir): + err[p.__class__.__name__] = p + + sout.close() + if len(err) == 0: + return None + else: + return err diff --git a/python/casm/casm/vasp/io/incar.py b/python/casm/casm/vasp/io/incar.py index 617decbc1..63cc3dff4 100644 --- a/python/casm/casm/vasp/io/incar.py +++ b/python/casm/casm/vasp/io/incar.py @@ -1,3 +1,4 @@ + from __future__ import (absolute_import, division, print_function, unicode_literals) from builtins import * @@ -40,7 +41,6 @@ class Incar(object): """ The INCAR class contains: tags: a dict of all INCAR settings - All input tags and associated values are stored as key-value pairs in the dicionary called 'tags'. """ def __init__(self,filename, species=None, poscar=None, sort=True): @@ -208,4 +208,3 @@ def write(self, filename): else: incar_write.write('{} = {}\n'.format(tag.upper(),self.tags[tag])) incar_write.close() - diff --git a/python/casm/casm/vasp/io/poscar.py b/python/casm/casm/vasp/io/poscar.py index 2836f6865..f415d9f66 100644 --- a/python/casm/casm/vasp/io/poscar.py +++ b/python/casm/casm/vasp/io/poscar.py @@ -22,7 +22,7 @@ class Site: self.occ_alias = alias (POTCAR) name, empty string by default self.position = np.array coordinate """ - def __init__(self, cart, position, SD_FLAG = "", occupant = "", occ_alias = ""): + def __init__(self, cart, position, SD_FLAG="", occupant="", occ_alias=""): """ Site constructor """ self.cart = cart self.SD_FLAG = SD_FLAG @@ -390,7 +390,6 @@ def _read_basis_legacy(self,file): self.basis.append(Site(cart, np.array(pos), SD_FLAG, atom_type, atom_type)) - def update(self, species): """ Set self.type_atoms_alias and self.basis[x].alias according to Species dict. @@ -412,6 +411,3 @@ def update(self, species): base.mag = float(species[base.occupant].tags['MAGMOM']) return - - - diff --git a/python/casm/casm/vasp/relax.py b/python/casm/casm/vasp/relax.py index d3d3dad9a..32238e27c 100644 --- a/python/casm/casm/vasp/relax.py +++ b/python/casm/casm/vasp/relax.py @@ -269,7 +269,7 @@ def run(self): if result is None or self.not_converging(): # Check for actions that should be taken after the initial run if len(self.rundir) == 1: - if self.settings["fine_ngx"]: + if "fine_ngx" in self.settings: outcarfile = os.path.join(self.rundir[-1], "OUTCAR") if not os.path.isfile(outcarfile): # This is an error but I'm not sure what to do about it diff --git a/python/casm/casm/vaspwrapper/__init__.py b/python/casm/casm/vaspwrapper/__init__.py index 3c002b0e2..48fc18759 100644 --- a/python/casm/casm/vaspwrapper/__init__.py +++ b/python/casm/casm/vaspwrapper/__init__.py @@ -1,13 +1,11 @@ """A wrapper for running vasp through casm""" -from casm.vaspwrapper.vaspwrapper import VaspWrapperError, read_settings, write_settings, \ - vasp_input_file_names, read_properties +from casm.vaspwrapper.vaspwrapper import VaspWrapperError, write_settings, vasp_input_file_names from casm.vaspwrapper.relax import Relax from casm.vaspwrapper.converge import Converge __all__ = [ 'Relax', 'Converge', 'VaspWrapperError', - 'read_settings', 'write_settings', 'vasp_input_file_names' ] diff --git a/python/casm/casm/vaspwrapper/converge.py b/python/casm/casm/vaspwrapper/converge.py index 7591e3722..75bb69218 100644 --- a/python/casm/casm/vaspwrapper/converge.py +++ b/python/casm/casm/vaspwrapper/converge.py @@ -23,8 +23,9 @@ from casm import vasp, wrapper from casm.misc import noindent from casm.project import DirectoryStructure, ProjectSettings -from casm.vaspwrapper import VaspWrapperError, read_settings, write_settings, \ - vasp_input_file_names +from casm.vaspwrapper import VaspWrapperError, write_settings, vasp_input_file_names + +from casm.project.io import read_project_settings ### Globals ### VALID_PROP_TYPES = ["KPOINTS", "ENCUT", "NBANDS", "SIGMA"] @@ -163,7 +164,7 @@ def __init__(self, configdir=None, auto=True, sort=True, propdir=None, prop=None else: print(" Read settings from:", setfile) - self.settings = read_settings(setfile) + self.settings = read_project_settings(setfile) # add required keys to settings if not present if not "ncore" in self.settings: diff --git a/python/casm/casm/vaspwrapper/relax.py b/python/casm/casm/vaspwrapper/relax.py index bc2504947..286b0c3b7 100644 --- a/python/casm/casm/vaspwrapper/relax.py +++ b/python/casm/casm/vaspwrapper/relax.py @@ -19,8 +19,9 @@ from casm import vasp, wrapper from casm.misc import noindent from casm.project import DirectoryStructure, ProjectSettings -from casm.vaspwrapper import VaspWrapperError, read_settings, write_settings, \ - vasp_input_file_names +from casm.vaspwrapper import VaspWrapperError, write_settings, vasp_input_file_names + +from casm.project.io import read_project_settings class Relax(object): """The Relax class contains functions for setting up, executing, and parsing a VASP relaxation. @@ -133,7 +134,7 @@ def __init__(self, configdir=None, auto=True, sort=True): else: print(" Read settings from:", setfile) - self.settings = read_settings(setfile) + self.settings = read_project_settings(setfile) # set default settings if not present if not "ncore" in self.settings: @@ -296,18 +297,18 @@ def submit(self): print("Constructing a job") sys.stdout.flush() # construct a Job - job = Job(name=wrapper.jobname(self.configname),\ - account=self.settings["account"],\ - nodes=int(math.ceil(float(N)/float(self.settings["atom_per_proc"])/float(self.settings["ppn"]))),\ - ppn=int(self.settings["ppn"]),\ - walltime=self.settings["walltime"],\ - pmem=self.settings["pmem"],\ - qos=self.settings["qos"],\ - queue=self.settings["queue"],\ - message=self.settings["message"],\ - email=self.settings["email"],\ - priority=self.settings["priority"],\ - command=cmd,\ + job = Job(name=wrapper.jobname(self.configname), + account=self.settings["account"], + nodes=self.settings["nodes"], + ppn=int(self.settings["ppn"]), + walltime=self.settings["walltime"], + pmem=self.settings["pmem"], + qos=self.settings["qos"], + queue=self.settings["queue"], + message=self.settings["message"], + email=self.settings["email"], + priority=self.settings["priority"], + command=cmd, auto=self.auto) print("Submitting") diff --git a/python/casm/casm/vaspwrapper/vaspwrapper.py b/python/casm/casm/vaspwrapper/vaspwrapper.py index 2abc50802..6768efa80 100644 --- a/python/casm/casm/vaspwrapper/vaspwrapper.py +++ b/python/casm/casm/vaspwrapper/vaspwrapper.py @@ -13,7 +13,7 @@ def __str__(self): return self.msg -def read_settings(filename): +def old_read_vasp_settings(filename): """Returns a JSON object reading JSON files containing settings for VASP PBS jobs. Returns: diff --git a/python/casm/casm/wrapper/misc.py b/python/casm/casm/wrapper/misc.py index acf60cc41..67e647467 100644 --- a/python/casm/casm/wrapper/misc.py +++ b/python/casm/casm/wrapper/misc.py @@ -1,10 +1,9 @@ -"""Functions used by wrappers""" -from __future__ import (absolute_import, division, print_function, unicode_literals) -from builtins import * +from __future__ import absolute_import, division, print_function, unicode_literals import os import re + def jobname(configname): """Return a name for a submitted job for configuration with 'configname' @@ -16,5 +15,17 @@ def jobname(configname): """ return configname.replace(os.sep, '.') + def remove_chars(val, chars): - return re.sub(' +', ' ', re.sub(chars,'',str(val).strip())) + return re.sub(' +', ' ', re.sub(chars, '', str(val).strip())) + + +def confname_as_jobname(configname): + """ + Returns only the configuration name (SCEL... etc) i.e. for job name in queue + :param + configname: Name of the configuration + :return + configuration name without path + """ + return str(configname.split('/')[-2]) + '/' + str(configname.split('/')[-1]) diff --git a/python/casm/requirements.txt b/python/casm/requirements.txt index b03df3cdb..97d1d59e1 100644 --- a/python/casm/requirements.txt +++ b/python/casm/requirements.txt @@ -8,3 +8,4 @@ scipy sh six tornado +atomate \ No newline at end of file diff --git a/python/casm/setup.py b/python/casm/setup.py index f08c21b3c..b17ff8389 100644 --- a/python/casm/setup.py +++ b/python/casm/setup.py @@ -6,7 +6,7 @@ # get console_scripts def script_str(file): name = os.path.splitext(os.path.split(file)[1])[0] - if name in ['casm_calc', 'casm_learn', 'casm_plot']: + if name in ['casm_calc', 'casm_learn', 'casm_plot', 'casm_bands', 'casm_feff']: return name.replace('_','-') + '=casm.scripts.' + name + ':main' else: return name.replace('_','.') + '=casm.scripts.' + name + ':main' @@ -33,7 +33,8 @@ def script_str(file): 'scikit-learn', 'scipy', 'sh', - 'tornado' + 'tornado', + 'atomate' ], classifiers=[ 'Development Status :: 4 - Beta', diff --git a/src/casm/clex/Configuration.cc b/src/casm/clex/Configuration.cc index f70852adc..cc06766a5 100644 --- a/src/casm/clex/Configuration.cc +++ b/src/casm/clex/Configuration.cc @@ -469,14 +469,15 @@ namespace CASM { //********************************************************************************* bool Configuration::read_calc_properties(jsonParser &parsed_props) const { - //std::cout << "begin Configuration::read_calculated()" << std::endl; + // std::cout << "begin Configuration::read_calculated()" << std::endl; bool success = true; /// properties.calc.json: contains calculated properties /// For default clex calctype only fs::path filepath = calc_properties_path(); - //std::cout << "filepath: " << filepath << std::endl; + // std::cout << "filepath: " << filepath << std::endl; parsed_props = jsonParser(); if(fs::exists(filepath)) { + jsonParser json(filepath); //Record file timestamp @@ -484,7 +485,7 @@ namespace CASM { std::vector props = get_primclex().settings().properties(); for(Index i = 0; i < props.size(); i++) { - //std::cout << "checking for: " << props[i] << std::endl; + // std::cout << "checking for: " << props[i] << std::endl; if(json.contains(props[i])) { // normal by #prim cells for some properties @@ -550,6 +551,7 @@ namespace CASM { else success = false; + std::cout << "RESULT: " << success << std::endl; return success; }