From 467e6af1b161a45cb8f225fde78c62caff7cc66a Mon Sep 17 00:00:00 2001 From: Tara Maria Boland Date: Wed, 13 Apr 2022 14:10:44 -0700 Subject: [PATCH] Electronic (#2) * adding electronic properites workflow * first commit of new workflow * introduce new code: automate routine for computing electronic and bader * update vasp from prev for dos and bader * update workflow package imports, add new electronic workflow, move functions to logical modules * create parsers, now in workflow testing phase * fix workflow run bugs * update icharg value * modify run_calc file and workflow to do a double kpoint run for chgcar * update gridfs for dos * update parser output * update electronic workflow * fix bug in average_z_sep * fix bug in average_z_sep * add correct name to 2d-substrate heterostructs based on stacked wyckoff sites * fix bugs * update parse output * update run calc and core to use custom custodian handler * debugging run_calc and appending cdd analysis to same db entry as dos,bader run * modify run_calc to set npar. custodian currently does not search for the correct slurm variable to set npar. so i fixed it * updated ToDb: dos->element projected with 0.05 eV res and default [4,6] eV window around efermi. add user option for window * update obj_id variable pass * fix auto_npar * push obj_id to cdd calculation spec * add flags to depends to trigger slurm_npar * remove parent cdd_combined from iso NSCF, clarify error message in workflow, add task_label to additional_fields not tags * update cdd error message * fix bug in vasp interface set * fix logic error for vis_static; structure now updates to the isolated structures * add logger to CDD parser * delete random } * clean up unneeded comments * include amin. add option to include dipole correction in nscf calculation. fix cmdlrelaxset's ediff value. reformat name for electronicFW if not cdd. * change out final structure from previous run is obtained. * remove amin * fix pushing obj_id for combined system so the cdd analysis can be pushed to that obj_id * remove conflicting isym tag * clean up formatting of documents * change task_name to task_label * update parser * fix parse output * fix parse output * fix parse output * fix entry of dos into grid fs * don't push None to fw spec * fix brackets for database pushing * cleaning up store_data for cdd update * add logger info to processing document clean up debug code, add exception errors * fix logic error in slurm_npar * update docs, add get_dos utility, update parser to put completedos into database. machine learning needs that object * fix Electronic to DB parser. add z proj chg for NSCF structures * update decompress to describe behaviour of function * fix stupid error * add complete dos to import for utils * add todo --- hetero2d/firetasks/heteroiface_tasks.py | 56 +-- hetero2d/firetasks/parse_outputs.py | 459 +++++++++++++++------ hetero2d/firetasks/run_calc.py | 269 ++++++++---- hetero2d/firetasks/write_inputs.py | 87 ++-- hetero2d/fireworks/core.py | 216 +++++++--- hetero2d/io/CMDLRelaxSet.yaml | 2 +- hetero2d/io/VaspInterfaceSet.py | 233 +++++++++-- hetero2d/io/__init__.py | 1 + hetero2d/manipulate/heterotransmuter.py | 183 ++++---- hetero2d/manipulate/layersolver.py | 6 +- hetero2d/manipulate/utils.py | 527 +++++++++++++++++------- hetero2d/workflow/core.py | 475 +++++++++++++++------ 12 files changed, 1758 insertions(+), 756 deletions(-) diff --git a/hetero2d/firetasks/heteroiface_tasks.py b/hetero2d/firetasks/heteroiface_tasks.py index 4b5e089..6424311 100644 --- a/hetero2d/firetasks/heteroiface_tasks.py +++ b/hetero2d/firetasks/heteroiface_tasks.py @@ -10,8 +10,7 @@ from __future__ import division, print_function, unicode_literals, absolute_import -import fnmatch, glob, gzip, json, os, re - +import glob, os, re from six.moves import range from copy import deepcopy from importlib import import_module @@ -19,12 +18,15 @@ from monty.json import jsanitize from monty.serialization import dumpfn -from atomate.utils.utils import get_logger, env_chk +from pymatgen import Structure + from fireworks.core.firework import FiretaskBase, FWAction from fireworks.utilities.fw_utilities import explicit_serialize -from pymatgen import Structure -from hetero2d.manipulate.utils import get_mongo_client +from atomate.utils.utils import get_logger, env_chk + +from hetero2d.manipulate.utils import get_mongo_client + __author__ = 'Tara M. Boland' __copyright__ = "Copyright 2020, CMD Lab" @@ -33,44 +35,6 @@ logger = get_logger(__name__) - -def find(pattern, path): - """ - Finds a file on a given search path that matches the - pattern. - - Args: - pattern (str): String to match on. - path (str): File path to search over. - - Returns: - a list of files matching the criteria - """ - result = [] - for root, dirs, files in os.walk(path): - for name in files: - if fnmatch.fnmatch(name, pattern): - result.append(name) - return result - - -def get_FWjson(): - """ - Helper function which reads the FW.json file in the - local directory and returns the FW.json file. - """ - fw_file = find('FW.json*', '.')[0] - - if fw_file == 'FW.json.gz': - with gzip.GzipFile('FW.json.gz', 'r') as p: - fw_json_file = json.loads(p.read().decode('utf-8')) - elif fw_file == 'FW.json': - with open('FW.json', 'rb') as p: - fw_json_file = json.load(p) - - return fw_json_file - - # Helper Functions def _update_spec(additional_spec): """ @@ -174,7 +138,6 @@ def _update_spec(additional_spec): "_dupefinder": {"_fw_name": "DupeFinderExact"}}) return default_spec - @explicit_serialize class TransferSlabTask(FiretaskBase): """ @@ -184,8 +147,8 @@ class TransferSlabTask(FiretaskBase): Args: label (str): The name of the structure in the additional_spec dict - struct_sub. - - Returns: + + Returns: FWAction updating the spec """ required_params = ['label'] @@ -202,7 +165,6 @@ def run_task(self, fw_spec): return FWAction(update_spec=update_spec) - @explicit_serialize class CreateHeterostructureTask(FiretaskBase): """ diff --git a/hetero2d/firetasks/parse_outputs.py b/hetero2d/firetasks/parse_outputs.py index 921f529..154a9b3 100644 --- a/hetero2d/firetasks/parse_outputs.py +++ b/hetero2d/firetasks/parse_outputs.py @@ -3,26 +3,39 @@ # Distributed under the terms of the GNU License. """ -These modules parse output from VASP runs, calculate various energetic and structural +These modules parse output from VASP runs, calculate various energetic and structural parameters and update a mongoDB with the relevant information. """ from __future__ import division, print_function, unicode_literals, absolute_import -import glob, os, re -from monty.json import jsanitize +import glob, os, re, gridfs, json, zlib, traceback +import numpy as np +from bson import ObjectId + +from monty.os.path import which +from monty.json import MontyEncoder, jsanitize from monty.serialization import dumpfn -from atomate.common.firetasks.glue_tasks import get_calc_loc -from atomate.utils.utils import get_logger, env_chk, get_meta_from_structure -from atomate.vasp.drones import VaspDrone +from pymatgen import Structure +from pymatgen.io.vasp import Potcar +from pymatgen.io.vasp.outputs import Vasprun +from pymatgen.io.vasp.sets import get_vasprun_outcar +from pymatgen.command_line.bader_caller import bader_analysis_from_path + from fireworks.core.firework import FiretaskBase, FWAction from fireworks.utilities.fw_serializers import DATETIME_HANDLER from fireworks.utilities.fw_utilities import explicit_serialize -from pymatgen import Structure -from hetero2d.firetasks.heteroiface_tasks import get_FWjson -from hetero2d.manipulate.utils import tag_iface, get_mongo_client +from atomate.common.firetasks.glue_tasks import get_calc_loc +from atomate.utils.utils import get_logger, env_chk, get_meta_from_structure +from atomate.vasp.drones import VaspDrone + +from hetero2d.manipulate.utils import tag_iface, get_mongo_client, get_FWjson +from hetero2d.manipulate.utils import vtotav, decompress + +#bader_exe_exists = which("bader") or which("bader.exe") +bader_exe_exists = '/home/tboland1/anaconda3/envs/cms/bin/bader' __author__ = 'Tara M. Boland' __copyright__ = "Copyright 2020, CMD Lab" @@ -31,16 +44,139 @@ logger = get_logger(__name__) +# Analysis +@explicit_serialize +class HeteroAnalysisToDb(FiretaskBase): + """ + Enter heterostructure workflow analysis into the database. Determines + what data to enter into the database based on what calculation was + performed. -def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=None, db_file=None): + Args: + db_file (str): Path to file containing the database credentials. + Supports env_chk. Default: write data to JSON file. + task_label (str): The task_label (firework name) automatically generated + by the FireWorks in Hetero2d. Used to determine what type of calculation + is being parsed. + + Optional params: + calc_dir (str): Path to dir (on current filesystem) that contains VASP + output files. Default: use current working directory. + calc_loc (str/bool): If True will set most recent calc_loc. If str + will search calc_locs for the matching name (most recent). + dos (bool/str): If True, parses the density of states assuming uniform + mode. Set to line if you are parsing bandstructures. Data stored + in GridFS. Default set to False. + bader (bool): If True, bader analysis is performed for the current + directory. The bader.exe must exist on the path. Default set to + False. + cdd (bool): If True, the charge density difference is performed. Default + set to False. + Adsorption_Energy (bool): If True, the adsorption formation energy + analysis is calculated for the run the results are stored with the + structure in the database under the Adsorption_Energy key. + Binding_Energy (bool): If True, the binding energy will be calculated + and stored in the database under the Binding Energy key. + Formation_Energy (bool): If True, the formation energy will be calculated + and stored with structure in the database under the Formation_Energy + key. + additional_fields (dict): Dict of additional fields to add to the database. """ - Insert a new doc for the 2d, 3d2d, bulk, substrate slab, 2d-subs generator, or 2d-subs - configuration into the database. + required_params = ["db_file", "task_label"] + optional_params = ["dos", "bader", "cdd", "Adsorption_Energy", "Binding_Energy", + "Formation_Energy", "additional_fields"] + + def run_task(self, fw_spec): + logger.info("Launching HeteroAnalysisToDb") + #try: + db_file = env_chk('>>db_file<<', fw_spec) + task_label = self.get("task_label", None) # get task_label + # determine what data to insert into the database + additional_fields = self.get("additional_fields", None) # update additional_fields + dos, bader, cdd = [self.get(i, False) for i in ['dos','bader','cdd']] + ## remove this after manual updates are done + #except: + # db_file = self.get('db_file', None) or env_chk('>>db_file<<', fw_spec) + # task_label = self.get("task_label", None) # get task_label + # additional_fields = self.get("additional_fields", None) # update additional_fields + # dos, bader, cdd = [additional_fields.get(i, False) for i in ['dos','bader','cdd']] + + # define input to push to mod_spec and stored data + if re.search("Optimization", task_label): + logger.info("Processing TaskDoc for Optimization Run.") + stored_data, mod_spec = {'analysis_info': {}}, [{'_push': {'analysis_info': {}}}] + info = {} + [info.update(element) for element in fw_spec.get('analysis_info', [{}])] + else: + stored_data, mod_spec = None, None + + # Bulk Substrate Optimization Analysis # + if re.search("Bulk Structure Optimization", task_label): + logger.info("PASSING PARAMETERS TO TASKDOC for Bulk") + struct_bulk, N_bulk, E_bulk = HeteroTaskDoc(self, fw_spec, + task_label, 'Bulk', additional_fields, + db_file) + + # Oriented Substrate Slab Optimization Analysis # + if re.search("Slab Structure Optimization", task_label): + logger.info("PASSING PARAMETERS TO TASKDOC for Slab") + struct_sub, N_sub, E_sub = HeteroTaskDoc(self, fw_spec, + task_label, 'Substrate', additional_fields, + db_file) + info.update({'E_sub': E_sub, 'N_sub': N_sub}) + stored_data["analysis_info"].update(info) + mod_spec[0]['_push']['analysis_info'].update(info) + + # 2D Structure Optimization Analysis # + if re.search("[^3D]2D Structure Optimization", task_label): + logger.info("PASSING PARAMETERS TO TASKDOC for 2D") + struct_2D, N_2D, E_2D = HeteroTaskDoc(self, fw_spec, task_label, '2D', + additional_fields, db_file) + stored_data["analysis_info"].update({'N_2d': struct_2D.num_sites, + 'E_2d': E_2D}) + mod_spec[0]['_push']['analysis_info'].update({'N_2d': struct_2D.num_sites, + 'E_2d': E_2D}) + + # 3D2D Structure Optimization Analysis # + if re.search("3D2D Structure Optimization", task_label): + logger.info("PASSING PARAMETERS TO TASKDOC for 3D2D") + struct_3D2D, N_3D2D, E_3D2D = HeteroTaskDoc(self, fw_spec, task_label, + '3D2D', additional_fields, db_file) + stored_data["analysis_info"].update({'N_3d2d': N_3D2D, + 'E_3d2d': E_3D2D}) + mod_spec[0]['_push']['analysis_info'].update({'N_3d2d': N_3D2D, + 'E_3d2d': E_3D2D}) + + # 2d on substrate Optimization Analysis # + if re.search("Heterostructure Optimization:", task_label): + logger.info("PASSING PARAMETERS TO TASKDOC for 2D_on_Substrate") + struct_2Dsub, N_2Dsub, E_2Dsub = HeteroTaskDoc(self, fw_spec, task_label, + "2D_on_Substrate", + additional_fields, + db_file) + + # Density of States and Bader Analysis # + if True in [dos, bader, cdd]: + logger.info("PASSING PARAMETERS TO TASKDOC for Electronic Property") + parse_vasp = False if cdd else True + obj_id = DosBaderTaskDoc(self, fw_spec, task_label, "DosBader", dos, + bader, cdd, parse_vasp, additional_fields, db_file) + if obj_id: + stored_data = {'obj_id': obj_id} + mod_spec = [{'_push': {'obj_id': obj_id}}] + return FWAction(stored_data=stored_data, mod_spec=mod_spec) + +def HeteroTaskDoc(self, fw_spec, task_label, task_collection, + additional_fields=None, db_file=None): + """ + Insert a new doc for the 2d, 3d2d, bulk, substrate slab, 2d-subs + generator, or 2d-subs configuration into the database. Args: + self (self): The self parameter for HeteroAnalysisToDb. fw_spec (dict): A dictionary containing all the information linked to this firework. - task_name (str): The name of the firework being analyzed. + task_label (str): The name of the firework being analyzed. task_collection (str): The name of the task collection you want to push the data to. additional_fields (dict): a dictionary containing additional @@ -51,8 +187,6 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N Returns: Analyzed structure, number of sites, and energy """ - dumpfn(fw_spec, 'fw_spec.json') - # get directory info calc_dir = os.getcwd() if "calc_dir" in self: @@ -76,7 +210,7 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N final_struct = Structure.from_dict(task_doc["calcs_reversed"][0]["output"]['structure']) # get the actual initial POSCAR for database f = glob.glob("POSCAR.orig*")[0] - # this takes the poscar and preps it for database doc + # this takes the poscar and preps it for database doc init_struct = Structure.from_file(f, False) cif = final_struct.to(fmt='cif') N = final_struct.num_sites @@ -84,12 +218,14 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N info = {} [info.update(element) for element in fw_spec.get('analysis_info', [{}])] - # standard database information + # standard database information heterostructure_dict = {'compound': task_doc['formula_pretty'], - 'dir_name': task_doc['dir_name'], 'fw_id': fw_id, 'task_label': task_name, + 'dir_name': task_doc['dir_name'], 'fw_id': fw_id, + 'task_label': task_label, 'final_energy': E, 'initial_structure': init_struct.as_dict(), 'final_structure': final_struct.as_dict(), 'cif': cif, - 'analysis_data': info, "metadata": get_meta_from_structure(structure=final_struct)} + 'analysis_data': info, + "metadata": get_meta_from_structure(structure=final_struct)} ########################################## if task_collection == "2D_on_Substrate": @@ -107,7 +243,7 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N # additional energetic information # NOTE: never pull the *_Energy tags from fw_spec only - # pull them from self. + # pull them from self. Formation_Energy = self.get("Formation_Energy", None) Binding_Energy = self.get("Binding_Energy", None) Adsorption_Energy = self.get("Adsorption_Energy", None) @@ -127,7 +263,7 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N if not isinstance(struct_sub, Structure): struct_sub = Structure.from_dict(struct_sub) N_iface = final_struct.num_sites - # calc E_ads + # calc E_ads align_info = heterostructure_dict['tags']['alignment_info'] n_2d, n_sub = align_info['fu_2d'], align_info['fu_sub'] E_bind = (n_2d * info['E_2d'] + n_sub * info['E_sub'] - E) / (n_2d * info['N_2d']) @@ -146,7 +282,7 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N if not isinstance(struct_sub, Structure): struct_sub = Structure.from_dict(struct_sub) N_iface = final_struct.num_sites - # calc E_ads + # calc E_ads align_info = heterostructure_dict['tags']['alignment_info'] n_2d, n_sub = align_info['fu_2d'], align_info['fu_sub'] E_ads = E_form - (n_2d * info['E_2d'] + n_sub * info['E_sub'] - E) / (n_2d * info['N_2d']) @@ -171,110 +307,189 @@ def HeteroTaskDoc(self, fw_spec, task_name, task_collection, additional_fields=N col = db[task_collection] col.insert_one(heterostructure_dict) conn.close() - dumpfn(heterostructure_dict, 'heterostructure_doc.json') return final_struct, N, E - -# Analysis -@explicit_serialize -class HeteroAnalysisToDb(FiretaskBase): +def DosBaderTaskDoc(self, fw_spec, task_label, task_collection, dos, bader, + cdd, parse_vasp, additional_fields=None, db_file=None): """ - Enter heterostructure workflow analysis into the database. + Insert a new task document into task_collection specified by db_file + that contains dos, bader, and charge density difference for a set of + calculations. - Args: - db_file (str): path to file containing the database credentials. - Supports env_chk. Default: write data to JSON file. - wf_name (str): The name of the workflow that this analysis is part of. - - Other Parameters: - Adsorption_Energy (bool): If set this will perform adsorption energy - analysis for the run and send the results to the Adsorption_Energy - collection. - Binding_Energy (bool): If set the binding energy will be calculated - and sent to the Binding Energy collection. - Formation_Energy (bool): If set the formation energy will be calculated - and sent to the Formation_Energy collection. - calc_dir (str): The calculation directory to parse. - calc_loc (str): The location to the directory to parse. + Args: + fw_spec (dict): A dictionary containing all the information + linked to this firework. + task_label (str): The name of the firework being analyzed. + task_collection (str): The name of the task collection you + want to push the data to. + dos (bool/str): If True, parses the density of states assuming uniform + mode. Set to line if you are parsing bandstructures. Data stored + in GridFS. Default set to False. + bader (bool): If True, bader analysis is performed for the current + directory. The bader.exe must exist on the path. Default set to + False. + cdd (bool): If True, the charge density difference is performed. Default + set to False. + parse_vasp (bool): If True, parses the vasprun to add extra information + to the database. + additional_fields (dict): A dictionary containing additional + information you want added to the database. + db_file (str): a string representation for the location of + the database file. """ - required_params = ["db_file", "wf_name", "task_label"] - optional_params = ["Adsorption_Energy", "Binding_Energy", "Formation_Energy", 'additional_fields'] - - def run_task(self, fw_spec): - logger.info("Starting HeteroAnalysisToDb") - wf_name = self.get('wf_name') - db_file = env_chk('>>db_file<<', fw_spec) - - # determine what data to insert into the database - task_label = self.get("task_label", None) - additional_fields = self.get("additional_fields", None) - - # initialize data arrays - stored_data = {'analysis_info': {}} - mod_spec = [{'_push': {'analysis_info': {}}}] - analysis = None - - # ensure that the analysis_info dictionary is not - # a list of dictionaries after updating - info = {} - [info.update(element) for element in fw_spec.get('analysis_info', [{}])] - - print('### PRINTING TASK_LABEL', task_label) - ################################################# - # Bulk Substrate Optimization Analysis # - if re.search("Bulk Structure Optimization", task_label): - logger.info("PASSING PARAMETERS TO TASKDOC: Bulk") - struct_bulk, N_bulk, E_bulk = HeteroTaskDoc(self, fw_spec, - task_label, 'Bulk', additional_fields, - db_file) - ################################################# - - ################################################# - # Oriented Substrate Slab Optimization Analysis # - if re.search("Slab Structure Optimization", task_label): - logger.info("PASSING PARAMETERS TO TASKDOC: Slab") - struct_sub, N_sub, E_sub = HeteroTaskDoc(self, fw_spec, - task_label, 'Substrate', additional_fields, - db_file) - info.update({'E_sub': E_sub, 'N_sub': N_sub}) - stored_data["analysis_info"].update(info) - mod_spec[0]['_push']['analysis_info'].update(info) - ################################################# - - ################################################# - # 2D Structure Optimization Analysis # - if re.search(" 2D Structure Optimization", task_label): - logger.info("PASSING PARAMETERS TO TASKDOC: 2D") - struct_2D, N_2D, E_2D = HeteroTaskDoc(self, fw_spec, - task_label, '2D', additional_fields, - db_file) - stored_data["analysis_info"].update({'N_2d': struct_2D.num_sites, - 'E_2d': E_2D}) - mod_spec[0]['_push']['analysis_info'].update({'N_2d': struct_2D.num_sites, - 'E_2d': E_2D}) - ################################################# - - ################################################# - # 3D2D Structure Optimization Analysis # - if re.search("3D2D Structure Optimization", task_label): - logger.info("PASSING PARAMETERS TO TASKDOC: 3D2D") - struct_3D2D, N_3D2D, E_3D2D = HeteroTaskDoc(self, fw_spec, - task_label, '3D2D', additional_fields, - db_file) - stored_data["analysis_info"].update({'N_3d2d': N_3D2D, - 'E_3d2d': E_3D2D}) - mod_spec[0]['_push']['analysis_info'].update({'N_3d2d': N_3D2D, - 'E_3d2d': E_3D2D}) - ################################################# - - ################################################# - # 2d on substrate Optimization Analysis # - if re.search("Heterostructure Optimization:", task_label): - logger.info("PASSING PARAMETERS TO TASKDOC: 2D_on_Substrate") - struct_2Dsub, N_2Dsub, E_2Dsub = HeteroTaskDoc(self, fw_spec, - task_label, "2D_on_Substrate", additional_fields, - db_file) - ################################################# + # get directory info + if "calc_dir" in self: # passed calc dir + calc_dir = self["calc_dir"] + elif self.get("calc_loc"): # find the calc_loc in fw_spec + calc_dir = get_calc_loc(self["calc_loc"], fw_spec["calc_locs"])["path"] - return FWAction(stored_data=stored_data, mod_spec=mod_spec) + ## delete this after done updating bader + #elif fw_spec.get('_launch_dir'): + # #print(fw_spec['_launch_dir']) + # calc_dir = fw_spec['_launch_dir'] + #elif fw_spec.get('launches'): + # if fw_spec['launches'][0].get('launch_dir'): + # #print(fw_spec['launches'][0]['launch_dir']) + # calc_dir = fw_spec['launches'][0]['launch_dir'] + + # the main taskdoc, the cdd taskdoc, and the dos taskdoc + store_doc, cdd_dict, dos_dict = {}, {}, {} # store data doc + store_doc.update({'fw_id': get_FWjson()['fw_id']}) # get fw_id + # TASKDOC: Parse Vasprun + if parse_vasp: + logger.info("PARSING DIRECTORY with VaspDrone: {}".format(calc_dir)) + drone = VaspDrone() + task_doc = drone.assimilate(calc_dir) + vrun, outcar = get_vasprun_outcar('.') + doc_keys = ['dir_name', 'run_stats', 'chemsys', 'formula_reduced_abc', 'completed_at', + 'nsites', 'composition_unit_cell', 'composition_reduced', 'formula_pretty', + 'elements', 'nelements', 'input', 'last_updated', 'custodian', 'orig_inputs', + 'output'] + store_doc.update({key: task_doc[key] for key in doc_keys}) + + # TASKDOC: DOS processing + if dos: + logger.info("Processing DOS") + try: + # element projected dos from complete dos + dos_dict = vrun.complete_dos.as_dict() + store_doc['get_dos'] = True + except Exception: + raise ValueError("No valid dos data exist") + + # TASKDOC: Bader processing + if bader: + logger.info("Processing Bader. Executable path: {}".format(str(bader_exe_exists))) + ba = bader_analysis_from_path(path=calc_dir) + structure = Structure.from_dict(task_doc["calcs_reversed"][0]["output"]['structure']) + potcar = Potcar.from_file(filename=os.path.join(calc_dir, 'POTCAR.gz')) + nelectrons = {p.element: p.nelectrons for p in potcar} + ba['species'] = [s.species_string for s in structure.sites] + ba['nelectrons'] = [nelectrons[s] for s in ba['species']] + ba['zcoords'] = [s.coords[2] for s in structure.sites] + store_doc.update(ba) + + # TASKDOC: Z-projected Charge Density + if re.search('NSCF:', task_label) and cdd != True: + logger.info("Computing Z-projected Charge Density: {}".format(calc_dir)) + try: + # compute the projected charge density difference, assume uncompressed + store_data.update({'z_proj_chg': vtotav('CHGCAR')}) + except: + # decompress CHGCAR + decompress(glob.glob("CHGCAR*")[0], 'CHGCAR_projz') + # compute the projected charge density difference + store_data.update({'z_proj_chg': vtotav('CHGCAR_projz')}) + os.remove('CHGCAR_projz') + + # TASKDOC: Charge Density Difference processing + if cdd: + logger.info("Computing Charge Density Difference: {}".format(calc_dir)) + calc_locs = fw_spec['calc_locs'] + + # get the path for each of the structures for CDD + iso1_path = [loc['path'] for loc in calc_locs if re.search('ISO 1 NSCF:', loc['name'])][0] + iso2_path = [loc['path'] for loc in calc_locs if re.search('ISO 2 NSCF:', loc['name'])][0] + comb_path = [loc['path'] for loc in calc_locs if re.search('Combined NSCF:', loc['name'])][0] + # decompress the files + decompress(glob.glob(os.path.join(iso1_path, "CHGCAR*"))[0], 'CHGCAR_1') + decompress(glob.glob(os.path.join(iso2_path, "CHGCAR*"))[0], 'CHGCAR_2') + decompress(glob.glob(os.path.join(comb_path, "CHGCAR*"))[0], 'CHGCAR_comb') + + # compute the projected charge density difference + chg1, chg2, chg_comb = vtotav('CHGCAR_1'), vtotav('CHGCAR_2'), vtotav('CHGCAR_comb') + chg_cdd = list(np.array(chg_comb['chg_density']) - \ + (np.array(chg1['chg_density']) + np.array(chg2['chg_density']))) + # remove all the files you created + [os.remove(i) for i in ['CHGCAR_1','CHGCAR_2','CHGCAR_comb']] + + # put it in a dictionary + cdd_dict = {'chg_density_diff': {'distance': chg1['distance'], 'chg_density': chg_cdd}} + + # TASKDOC: Add additional information to task doc, clean it up + if additional_fields: + for key, value in additional_fields.items(): + store_doc[key] = value + electronic_dict = jsanitize(store_doc) # sanitize database info + + # connect to database & insert electronic_dict + conn, database = get_mongo_client(db_file, db_type=fw_spec.get('db_type', None)) + db = conn[database] + col = db[task_collection] + + # MAIN TASKDOC: Insert Data into Database + if any([dos, bader, parse_vasp]): + try: + t_id = col.insert(electronic_dict) + logger.info('TASKDOC: ElectronicFW taskdoc upload.') + except: + traceback.print_exc() + raise ValueError('Failed to insert electronic_dict to the DB') + ## DOS TASKDOC: DOS into the GridFS in DB ## + if dos: + # LOOKUP: DosBader._id=dos_fs.files._id=dos_fs.chunks.files_id + + # Putting object_id in the metadata subdocument as per mongo specs + m_data = {"compression": "zlib", "task_label": task_label} + # always perform the string conversion when inserting directly to gridfs + d = json.dumps(dos_dict, cls=MontyEncoder) + d = zlib.compress(d.encode(), True) + # connect to gridFS + fs = gridfs.GridFS(db, "dos_fs") + try: + fs_id = fs.put(d, _id=t_id, metadata=m_data) + # insert into gridfs + col.update_one({"task_id": t_id}, {"$set": {"dos_compression": "zlib"}}) + col.update_one({"task_id": t_id}, {"$set": {"dos_fs_id": fs_id}}) + logger.info('DOS inserted into GridFS.') + except: + raise ValueError('Failed to insert DOS into GridFS') + ## CDD TASKDOC: Append to the Combined system task doc ## + if cdd: + # get the workflow to update CDD with + obj_id = fw_spec.get('obj_id', None) + if obj_id: + logger.info('CDD: Updating {} with charge density difference'.format(obj_id[0])) + # update the workflow or write to the directory + try: + col.update_one({"_id": ObjectId(obj_id[0])}, {"$set": cdd_dict } ) + except Exception: + with open('cdd_dict.json', 'w') as f: + f.write(json.dumps(cdd_dict, default=DATETIME_HANDLER)) + raise ValueError('Failed to insert CDD into previous task_doc') + else: + logger.info('### WARNING CDD NOT INSERTED INTO PREVIOUS TASK DOC: No ObjectId found.') + with open('cdd_dict.json', 'w') as f: + f.write(json.dumps(cdd_dict, default=DATETIME_HANDLER)) + # dump analysis to directory in case something happens. just re-upload this + # file to the db. + analysis = {**cdd_dict, **dos_dict, **electronic_dict} + with open('electronic_property.json', 'w') as f: + f.write(json.dumps(analysis, default=DATETIME_HANDLER)) + # return ObjectId if the Combined system + if re.search('Combined NSCF:', task_label): + return str(t_id) + else: + return None + conn.close() diff --git a/hetero2d/firetasks/run_calc.py b/hetero2d/firetasks/run_calc.py index a9723cb..bb12591 100644 --- a/hetero2d/firetasks/run_calc.py +++ b/hetero2d/firetasks/run_calc.py @@ -3,27 +3,33 @@ # Distributed under the terms of the GNU License. """ -This module defines tasks that support running vasp in various ways. -Under active development for Hetero2d. Currently not in use as the default -behaviour from atomate is sufficient. +This module defines tasks that support running vasp jobs for a specific purpose. +Under active development for Hetero2d. Simple modifications to atomate functions +to modify the default behaviour. """ -import os -import shlex +import os, re, shlex, numpy as np + +from pymatgen.core.structure import Structure +from pymatgen.io.vasp.inputs import Incar, Kpoints, Poscar, VaspInput -from atomate.utils.utils import env_chk, get_logger -from atomate.vasp.config import HALF_KPOINTS_FIRST_RELAX from custodian import Custodian +from custodian.vasp.jobs import VaspJob +from custodian.vasp.validators import VasprunXMLValidator, VaspFilesValidator from custodian.vasp.handlers import VaspErrorHandler, AliasingErrorHandler, \ MeshSymmetryErrorHandler, UnconvergedErrorHandler, MaxForceErrorHandler, \ FrozenJobErrorHandler, NonConvergingErrorHandler, PositiveEnergyErrorHandler, \ WalltimeHandler, StdErrHandler, DriftErrorHandler -from custodian.vasp.jobs import VaspJob, VaspNEBJob -from custodian.vasp.validators import VasprunXMLValidator, VaspFilesValidator + from fireworks import explicit_serialize, FiretaskBase, FWAction + +from atomate.utils.utils import env_chk, get_logger + from monty.os.path import zpath from monty.serialization import loadfn +from hetero2d.manipulate.utils import slurm_set_npar + __author__ = 'Tara M. Boland ' __credits__ = 'Shyue Ping Ong ' __copyright__ = "Copyright 2020, CMD Lab" @@ -32,39 +38,184 @@ logger = get_logger(__name__) +@explicit_serialize +class ElectronicJob(VaspJob): + """ + A subclass of VaspJob created to run a double kpoint run for calculating bader + and DOS. Just runs whatever is in the directory. But conceivably can be a complex + processing of inputs etc. with initialization. + """ + + def __init__(self, vasp_cmd, output_file="vasp.out", stderr_file="std_err.txt", + suffix="", final=True, backup=True, auto_npar=False, auto_gamma=True, + settings_override=None, gamma_vasp_cmd=None, copy_magmom=False, auto_continue=False): + """ + This constructor is necessarily complex due to the need for + flexibility. For standard kinds of runs, it's often better to use one + of the static constructors. The defaults are usually fine too. + + Args: + vasp_cmd (str): Command to run vasp as a list of args. For example, + if you are using mpirun, it can be something like + ["mpirun", "pvasp.5.2.11"] + output_file (str): Name of file to direct standard out to. + Defaults to "vasp.out". + stderr_file (str): Name of file to direct standard error to. + Defaults to "std_err.txt". + suffix (str): A suffix to be appended to the final output. E.g., + to rename all VASP output from say vasp.out to + vasp.out.relax1, provide ".relax1" as the suffix. + final (bool): Indicating whether this is the final vasp job in a + series. Defaults to True. + backup (bool): Whether to backup the initial input files. If True, + the INCAR, KPOINTS, POSCAR and POTCAR will be copied with a + ".orig" appended. Defaults to True. + auto_npar (bool): Whether to automatically tune NPAR to be sqrt( + number of cores) as recommended by VASP for DFT calculations. + Generally, this results in significant speedups. Defaults to + True. Set to False for HF, GW and RPA calculations. + auto_gamma (bool): Whether to automatically check if run is a + Gamma 1x1x1 run, and whether a Gamma optimized version of + VASP exists with ".gamma" appended to the name of the VASP + executable (typical setup in many systems). If so, run the + gamma optimized version of VASP instead of regular VASP. You + can also specify the gamma vasp command using the + gamma_vasp_cmd argument if the command is named differently. + settings_override ([dict]): An ansible style list of dict to + override changes. For example, to set ISTART=1 for subsequent + runs and to copy the CONTCAR to the POSCAR, you will provide:: + [{"dict": "INCAR", "action": {"_set": {"ISTART": 1}}}, + {"file": "CONTCAR", + "action": {"_file_copy": {"dest": "POSCAR"}}}] + gamma_vasp_cmd (str): Command for gamma vasp version when + auto_gamma is True. Should follow the list style of + subprocess. Defaults to None, which means ".gamma" is added + to the last argument of the standard vasp_cmd. + copy_magmom (bool): Whether to copy the final magmom from the + OUTCAR to the next INCAR. Useful for multi-relaxation runs + where the CHGCAR and WAVECAR are sometimes deleted (due to + changes in fft grid, etc.). Only applies to non-final runs. + auto_continue (bool): Whether to automatically continue a run + if a STOPCAR is present. This is very useful if using the + wall-time handler which will write a read-only STOPCAR to + prevent VASP from deleting it once it finishes + """ + self.vasp_cmd = vasp_cmd + self.output_file = output_file + self.stderr_file = stderr_file + self.final = final + self.backup = backup + self.suffix = suffix + self.settings_override = settings_override + self.auto_npar = auto_npar + self.auto_gamma = auto_gamma + self.gamma_vasp_cmd = gamma_vasp_cmd + self.copy_magmom = copy_magmom + self.auto_continue = auto_continue + super().__init__(vasp_cmd, output_file, stderr_file, final, backup, suffix, settings_override, + auto_npar, auto_gamma, gamma_vasp_cmd, copy_magmom, auto_continue) + + @classmethod + def double_kpoints_run(cls, vasp_cmd, auto_npar=True, half_kpts_first=True, + auto_continue=False, slurm_npar=False): + """ + Returns a list of two jobs the first is to obtain the CHGCAR and WAVECAR + with a low kp number using ICHARG = 1 with increased bands, charge density + grid, and nedos and the second calcaluation with ICHARG = 11 with increased + kp grid. + + Args: + vasp_cmd (str): Command to run vasp as a list of args. For example, + if you are using mpirun, it can be something like + ["mpirun", "pvasp.5.2.11"] + auto_npar (bool): Whether to automatically tune NPAR to be sqrt( + number of cores) as recommended by VASP for DFT calculations. + Generally, this results in significant speedups. Defaults to + True. + half_kpts_first (bool): Whether to halve the kpoint grid for the + first run to obtain convergence chgcar file. Speeds up calculation + time considerably for bader analysis. Defaults to True. + slurm_npar (bool): Use hetero2d to set npar. Use this if you are on + a slurm based computing resource. auto_npar does not work for + slurm computing arch. Default set to false. If set to true remove + auto_npar from my_fworker.yaml. + + Returns: + List of two jobs corresponding to an AFLOW style run. + """ + incar_orig = Incar.from_file('INCAR') + incar1 = {"ICHARG": 1, "LAECHG": False, "NEDOS": 301} + incar2 = {key: incar_orig.get(key) for key in incar1.keys() if key in incar_orig.keys()} + if slurm_npar: + npar = slurm_set_npar() + incar1['NPAR'] = npar + incar2['NPAR'] = npar + # set the override settings + settings_overide_1 = [{"dict": "INCAR", "action": {"_set": incar1}}] + settings_overide_2 = [{"dict": "INCAR", "action": {"_set": incar2}}] + if half_kpts_first and os.path.exists("KPOINTS") and os.path.exists("POSCAR"): + kpts = Kpoints.from_file("KPOINTS") + orig_kpts_dict = kpts.as_dict() + # lattice vectors with length < 8 will get >1 KPOINT + kpts.kpts = np.round(np.maximum(np.array(kpts.kpts) / 2, 1)).astype(int).tolist() + low_kpts_dict = kpts.as_dict() + settings_overide_1.append({"dict": "KPOINTS", "action": {"_set": low_kpts_dict}}) + settings_overide_2.append({"dict": "KPOINTS", "action": {"_set": orig_kpts_dict}}) + + return [ + VaspJob( + vasp_cmd, + final=False, + suffix=".kpoints1", + auto_npar=auto_npar, + auto_continue=auto_continue, + settings_override=settings_overide_1, + ), + VaspJob( + vasp_cmd, + final=True, + backup=False, + suffix="", + auto_npar=auto_npar, + auto_continue=auto_continue, + settings_override=settings_overide_2, + ), + ] @explicit_serialize -class RunVaspCustodian(FiretaskBase): +class RunElectronicCustodian(FiretaskBase): """ - Run VASP using custodian "on rails", i.e. in a simple way that supports most common options. + Modified version of RunVaspCustodian that runs VASP using custodian "on rails" but + taylored to run double runs for computing electronic structure properties. Convergences + the chgcar using small kp grid with bader level ngif grid density then an increased grid + density. Args: vasp_cmd (str): the name of the full executable for running VASP. Supports env_chk. Other Parameters: - job_type: (str) - choose from "normal" (default), "double_relaxation_run" (two consecutive - jobs), "full_opt_run" (multiple optimizations), and "neb" - handler_group: (str or [ErrorHandler]) - group of handlers to use. See handler_groups dict in the code for - the groups and complete list of handlers in each group. Alternatively, you can - specify a list of ErrorHandler objects. - max_force_threshold: (float) - if >0, adds MaxForceErrorHandler. Not recommended for - nscf runs. + job_type: (str) - "double_kpoints_run" (two consecutive jobs) - default. + handler_group: (str or [ErrorHandler]) - group of handlers to use. See handler_groups dict + in the code for the groups and complete list of handlers in each group. Alternatively, + you can specify a list of ErrorHandler objects. scratch_dir: (str) - if specified, uses this directory as the root scratch dir. Supports env_chk. gzip_output: (bool) - gzip output (default=T) max_errors: (int) - maximum # of errors to fix before giving up (default=5) - ediffg: (float) shortcut for setting EDIFFG in special custodian jobs auto_npar: (bool) - use auto_npar (default=F). Recommended set to T for single-node jobs only. Supports env_chk. gamma_vasp_cmd: (str) - cmd for Gamma-optimized VASP compilation. Supports env_chk. wall_time (int): Total wall time in seconds. Activates WallTimeHandler if set. - half_kpts_first_relax (bool): Use half the k-points for the first relaxation + half_kpts_first (bool): Use half the k-points first to converge charge density. + slurm_npar (bool): Use to set npar when using multiple nodes. Use this to get the + total number of tasks on a slurm queue system. auto_npar only works with 1 node. + Default set to false. """ required_params = ["vasp_cmd"] - optional_params = ["job_type", "handler_group", "max_force_threshold", "scratch_dir", - "gzip_output", "max_errors", "ediffg", "auto_npar", "gamma_vasp_cmd", - "wall_time", "half_kpts_first_relax"] + optional_params = ["job_type", "handler_group", "scratch_dir", "gzip_output", + "max_errors", "auto_npar", "gamma_vasp_cmd", "slurm_npar", + "wall_time", "half_kpts_first"] def run_task(self, fw_spec): @@ -90,94 +241,38 @@ def run_task(self, fw_spec): vasp_cmd = shlex.split(vasp_cmd) # initialize variables - job_type = self.get("job_type", "normal") + job_type = self.get("job_type", "double_kpoints_run") scratch_dir = env_chk(self.get("scratch_dir"), fw_spec) gzip_output = self.get("gzip_output", True) - max_errors = self.get("max_errors", CUSTODIAN_MAX_ERRORS) + max_errors = self.get("max_errors", 5) auto_npar = env_chk(self.get("auto_npar"), fw_spec, strict=False, default=False) + slurm_npar = self.get("slurm_npar", False) gamma_vasp_cmd = env_chk(self.get("gamma_vasp_cmd"), fw_spec, strict=False, default=None) if gamma_vasp_cmd: gamma_vasp_cmd = shlex.split(gamma_vasp_cmd) # construct jobs - if job_type == "normal": - jobs = [VaspJob(vasp_cmd, auto_npar=auto_npar, gamma_vasp_cmd=gamma_vasp_cmd)] - elif job_type == "double_relaxation_run": - jobs = VaspJob.double_relaxation_run(vasp_cmd, auto_npar=auto_npar, - ediffg=self.get("ediffg"), - half_kpts_first_relax=self.get("half_kpts_first_relax", - HALF_KPOINTS_FIRST_RELAX)) - elif job_type == "metagga_opt_run": - jobs = VaspJob.metagga_opt_run(vasp_cmd, auto_npar=auto_npar, - ediffg=self.get("ediffg"), - half_kpts_first_relax=self.get("half_kpts_first_relax", - HALF_KPOINTS_FIRST_RELAX)) - - elif job_type == "full_opt_run": - jobs = VaspJob.full_opt_run(vasp_cmd, auto_npar=auto_npar, - ediffg=self.get("ediffg"), - max_steps=9, - half_kpts_first_relax=self.get("half_kpts_first_relax", - HALF_KPOINTS_FIRST_RELAX)) - elif job_type == "neb": - # TODO: @shyuep @HanmeiTang This means that NEB can only be run (i) in reservation mode - # and (ii) when the queueadapter parameter is overridden and (iii) the queue adapter - # has a convention for nnodes (with that name). Can't the number of nodes be made a - # parameter that the user sets differently? e.g., fw_spec["neb_nnodes"] must be set - # when setting job_type=NEB? Then someone can use this feature in non-reservation - # mode and without this complication. -computron - nnodes = int(fw_spec["_queueadapter"]["nnodes"]) - - # TODO: @shyuep @HanmeiTang - I am not sure what the code below is doing. It looks like - # it is trying to override the number of processors. But I tried running the code - # below after setting "vasp_cmd = 'mpirun -n 16 vasp'" and the code fails. - # (i) Is this expecting an array vasp_cmd rather than String? If so, that's opposite to - # the rest of this task's convention and documentation - # (ii) can we get rid of this hacking in the first place? e.g., allowing the user to - # separately set the NEB_VASP_CMD as an env_variable and not rewriting the command - # inside this. - # -computron - - # Index the tag "-n" or "-np" - index = [i for i, s in enumerate(vasp_cmd) if '-n' in s] - ppn = int(vasp_cmd[index[0] + 1]) - vasp_cmd[index[0] + 1] = str(nnodes * ppn) - - # Do the same for gamma_vasp_cmd - if gamma_vasp_cmd: - index = [i for i, s in enumerate(gamma_vasp_cmd) if '-n' in s] - ppn = int(gamma_vasp_cmd[index[0] + 1]) - gamma_vasp_cmd[index[0] + 1] = str(nnodes * ppn) - - jobs = [VaspNEBJob(vasp_cmd, final=False, auto_npar=auto_npar, - gamma_vasp_cmd=gamma_vasp_cmd)] + if job_type == "double_kpoints_run": + jobs = ElectronicJob.double_kpoints_run(vasp_cmd, auto_npar=auto_npar, + slurm_npar=slurm_npar, half_kpts_first=self.get("half_kpts_first", True)) else: raise ValueError("Unsupported job type: {}".format(job_type)) # construct handlers - handler_group = self.get("handler_group", "default") if isinstance(handler_group, str): handlers = handler_groups[handler_group] else: handlers = handler_group - if self.get("max_force_threshold"): - handlers.append(MaxForceErrorHandler(max_force_threshold=self["max_force_threshold"])) - if self.get("wall_time"): handlers.append(WalltimeHandler(wall_time=self["wall_time"])) - - if job_type == "neb": - validators = [] # CINEB vasprun.xml sometimes incomplete, file structure different - else: - validators = [VasprunXMLValidator(), VaspFilesValidator()] - + validators = [VasprunXMLValidator(), VaspFilesValidator()] c = Custodian(handlers, jobs, validators=validators, max_errors=max_errors, scratch_dir=scratch_dir, gzipped_output=gzip_output) - c.run() if os.path.exists(zpath("custodian.json")): stored_custodian_data = {"custodian": loadfn(zpath("custodian.json"))} return FWAction(stored_data=stored_custodian_data) + diff --git a/hetero2d/firetasks/write_inputs.py b/hetero2d/firetasks/write_inputs.py index c2f052b..643115b 100644 --- a/hetero2d/firetasks/write_inputs.py +++ b/hetero2d/firetasks/write_inputs.py @@ -5,29 +5,24 @@ """ These classes write the vasp input sets used to control VASP tags. """ - from __future__ import division, print_function, unicode_literals, absolute_import -import glob, shutil, os, numpy as np +import os, numpy as np +from importlib import import_module -from pymatgen import Lattice, Structure, Specie -from pymatgen.core.surface import Slab, SlabGenerator +from pymatgen import Structure +from pymatgen.io.vasp import Poscar from pymatgen.alchemy.materials import TransformedStructure from pymatgen.alchemy.transmuters import StandardTransmuter -from pymatgen.io.vasp import Incar, Poscar, Potcar, PotcarSingle from pymatgen.transformations.site_transformations import AddSitePropertyTransformation -from importlib import import_module -from monty.serialization import dumpfn, loadfn - -from atomate.utils.utils import get_logger - -from fireworks import Firework from fireworks.core.firework import FiretaskBase from fireworks.utilities.fw_utilities import explicit_serialize -# mp_interfaces +from atomate.utils.utils import get_logger + +from hetero2d.io import CMDLElectronicSet from hetero2d.manipulate.utils import set_sd_flags -from hetero2d.manipulate.heterotransmuter import hetero_interfaces + __author__ = 'Tara M. Boland' __email__ = 'tboland1@asu.edu' @@ -36,6 +31,48 @@ logger = get_logger(__name__) +@explicit_serialize +class WriteVaspElectronicFromPrev(FiretaskBase): + """ + Writes input files to perform bader analysis, density of states, and charge + density difference calculations. + + Args: + dedos (float): Automatically set nedos using the total energy range + which will be divided by the energy step dedos. Default 0.05 eV. + grid_density (float): Distance between grid points for the NGXF,Y,Z grids. + Defaults to 0.03 Angs; NGXF,Y,Z are ~2x > default. For charge + density difference calculations the parent grid density is used for all + children fireworks. + dos (bool): If True, sets INCAR tags for high quality site-orbital projected + density of states. Defaults to True. + bader (bool): If True, sets INCAR tags to generate bader analysis files. + cdd (bool): If True, ensures the grid density matches between the parent + and child Fireworks. Default set to False. + + Other Parameters: + **kwargs (keyword arguments): Dict of any keyword arguments supported by the + CMDLRelaxSet.from_prev_calc(). + + """ + required_params = ["dedos", "grid_density", "dos", "bader", "cdd"] + + optional_params = ["prev_calc_dir", "small_gap_multiply", "nbands_factor", "electronic_set_overrides"] + + def run_task(self, fw_spec): + # get previous calculation information and increase accuracy + vis_orig = CMDLElectronicSet.from_prev_calc(prev_calc_dir=self.get("prev_calc_dir", "."), + dedos=self.get("dedos", 0.05), + grid_density=self.get("grid_density", 0.03), + dos=self.get("dos", True ), + bader=self.get("bader", True), + cdd=self.get("cdd", False), + nbands_factor=self.get("nbands_factor", 1), + small_gap_multiply=self.get("small_gap_multiply", None), + **self.get("electronic_set_overrides", {})) + + vis_dict = vis_orig.as_dict() # make changes to vis in dict format + vis_orig.write_input(".") @explicit_serialize class WriteSlabStructureIOSet(FiretaskBase): @@ -46,17 +83,17 @@ class WriteSlabStructureIOSet(FiretaskBase): only the last structure in the list is used. Args: - structure (Structure): input structure. - transformations (list): list of names of transformation classes as defined - in the modules in pymatgen.transformations + structure (Structure): Input structure. + transformations (list): List of names of transformation classes as defined + in the modules in pymatgen.transformations. vasp_input_set (VaspInputSet): VASP input set. Other Parameters: - transformation_params (list): list of dicts where each dict specifies the + transformation_params (list): List of dicts where each dict specifies the input parameters to instantiate the transformation class in the transformations list. - override_default_vasp_params (dict): additional user input settings. - prev_calc_dir (str): path to previous calculation if using structure + override_default_vasp_params (dict): Additional user input settings. + prev_calc_dir (str): Path to previous calculation if using structure from another calculation. Returns: @@ -68,7 +105,7 @@ class WriteSlabStructureIOSet(FiretaskBase): def run_task(self, fw_spec): """ - Execute the transformations and write the input files for the calculation + Execute the transformations and write the input files for the calculation. """ transformations = [] transformation_params = self.get("transformation_params", @@ -108,10 +145,6 @@ def run_task(self, fw_spec): vis = vis_orig.__class__.from_dict(vis_dict) vis.write_input(".") - # uncomment if you need to debug - # dumpfn(transmuter.transformed_structures[-1], "transformations.json") - - @explicit_serialize class WriteHeteroStructureIOSet(FiretaskBase): """ @@ -124,9 +157,9 @@ class WriteHeteroStructureIOSet(FiretaskBase): parameters include all pymatgen structure objects. Other Parameters: - vasp_input_set (VaspInputSet): VASP input set for the transformed 2d on sub-slab - override_default_vasp_params (dict): additional user input settings. - user_incar_settings (dict): additional incar settings to add to the calculation. + vasp_input_set (VaspInputSet): VASP input set for the transformed 2d on sub-slab. + override_default_vasp_params (dict): Additional user input settings. + user_incar_settings (dict): Additional incar settings to add to the calculation. Returns: None diff --git a/hetero2d/fireworks/core.py b/hetero2d/fireworks/core.py index b760b3b..484cacd 100644 --- a/hetero2d/fireworks/core.py +++ b/hetero2d/fireworks/core.py @@ -9,25 +9,27 @@ """ from __future__ import division, print_function, unicode_literals, absolute_import +from monty.os.path import which -from atomate.utils.utils import get_logger -from atomate.vasp.config import HALF_KPOINTS_FIRST_RELAX, RELAX_MAX_FORCE, VASP_CMD, DB_FILE -from atomate.common.firetasks.glue_tasks import PassCalcLocs -from atomate.vasp.firetasks.run_calc import RunVaspCustodian -from atomate.vasp.firetasks.glue_tasks import CopyVaspOutputs -from atomate.vasp.firetasks.write_inputs import WriteVaspFromIOSet +from pymatgen import Structure from fireworks import Firework, FileWriteTask from fireworks.utilities.fw_utilities import get_slug -from pymatgen import Structure -from pymatgen.core.surface import Slab +from atomate.utils.utils import get_logger +from atomate.vasp.firetasks.run_calc import RunVaspCustodian +from atomate.vasp.firetasks.glue_tasks import CopyVaspOutputs +from atomate.vasp.firetasks.write_inputs import WriteVaspFromIOSet +from atomate.common.firetasks.glue_tasks import PassCalcLocs +from atomate.vasp.config import HALF_KPOINTS_FIRST_RELAX, RELAX_MAX_FORCE, VASP_CMD, DB_FILE -from hetero2d.firetasks.heteroiface_tasks import CreateHeterostructureTask +from hetero2d.io import CMDLInterfaceSet +from hetero2d.firetasks.run_calc import RunElectronicCustodian from hetero2d.firetasks.parse_outputs import HeteroAnalysisToDb -from hetero2d.firetasks.write_inputs import WriteHeteroStructureIOSet, WriteSlabStructureIOSet +from hetero2d.firetasks.heteroiface_tasks import CreateHeterostructureTask +from hetero2d.firetasks.write_inputs import WriteHeteroStructureIOSet, WriteSlabStructureIOSet,\ + WriteVaspElectronicFromPrev -from hetero2d.io.VaspInterfaceSet import CMDLInterfaceSet __author__ = 'Tara M. Boland' __copyright__ = "Copyright 2020, CMD Lab" @@ -36,17 +38,16 @@ logger = get_logger(__name__) -# handler group for 2d-substrates -from custodian.vasp.handlers import VaspErrorHandler, MeshSymmetryErrorHandler, UnconvergedErrorHandler, \ - NonConvergingErrorHandler, PositiveEnergyErrorHandler, FrozenJobErrorHandler, StdErrHandler -# Set up Custodian Error Handlers +# custom error handler +from custodian.vasp.handlers import VaspErrorHandler, MeshSymmetryErrorHandler,\ + UnconvergedErrorHandler, NonConvergingErrorHandler, PositiveEnergyErrorHandler, \ + FrozenJobErrorHandler, StdErrHandler subset = list(VaspErrorHandler.error_msgs.keys()) subset.remove('brions') handler = [VaspErrorHandler(errors_subset_to_catch=subset), MeshSymmetryErrorHandler(), UnconvergedErrorHandler(), NonConvergingErrorHandler(), PositiveEnergyErrorHandler(), FrozenJobErrorHandler(), StdErrHandler()] - class HeteroOptimizeFW(Firework): def __init__(self, spec, structure, name="Structure Optimization", vasp_input_set=None, user_incar_settings=None, vasp_cmd=VASP_CMD, @@ -60,26 +61,26 @@ def __init__(self, spec, structure, name="Structure Optimization", Args: structure (Structure): Input structure. - spec (dict): The specification parameters used to control the + spec (dict): The specification parameters used to control the workflow and pass variables. Other Parameters: name (str): Name for the Firework. user_incar_settings (dict): Input settings to update the settings for the vasp input set. - vasp_input_set (VaspInputSet): input set to use. Defaults to + vasp_input_set (VaspInputSet): input set to use. Defaults to CMDLInterfaceSet() if None. vasp_cmd (str): Command to run vasp. ediffg (float): Shortcut to set ediffg in certain jobs. - db_file (str): Path to file specifying db credentials to place output + db_file (str): Path to file specifying db credentials to place output parsing. force_gamma (bool): Force gamma centered kpoint generation. job_type (str): custodian job type (default "double_relaxation_run"). - max_force_threshold (float): max force on a site allowed at end; otherwise, + max_force_threshold (float): max force on a site allowed at end; otherwise, reject job. - auto_npar (bool or str): whether to set auto_npar. defaults to env_chk: + auto_npar (bool or str): whether to set auto_npar. defaults to env_chk: ">>auto_npar<<". - half_kpts_first_relax (bool): whether to use half the kpoints for the first + half_kpts_first_relax (bool): whether to use half the kpoints for the first relaxation. parents ([Firework]): Parents of this particular Firework. kwargs: Other kwargs that are passed to Firework.__init__. @@ -87,26 +88,26 @@ def __init__(self, spec, structure, name="Structure Optimization", name = "{}: {}".format(structure.composition.reduced_formula, name) user_incar_settings = user_incar_settings or None vasp_input_set = vasp_input_set or CMDLInterfaceSet(structure, - auto_dipole=True, user_incar_settings=user_incar_settings, + auto_dipole=True, + user_incar_settings=user_incar_settings, vdw='optB88', iface=True) - + if vasp_input_set.incar["ISIF"] in (0, 1, 2, 7) and job_type == "double_relaxation": warnings.warn( "A double relaxation run might not be appropriate with ISIF {}".format( vasp_input_set.incar["ISIF"])) t = [WriteVaspFromIOSet(structure=structure, vasp_input_set=vasp_input_set), - RunVaspCustodian(handler_group=handler, vasp_cmd=vasp_cmd, + RunVaspCustodian(handler_group=handler, vasp_cmd=vasp_cmd, job_type=job_type, ediffg=ediffg, max_force_threshold=max_force_threshold, auto_npar=auto_npar, wall_time=spec.get('wall_time', None), half_kpts_first_relax=half_kpts_first_relax), PassCalcLocs(name=name), - HeteroAnalysisToDb(wf_name=spec['wf_name'], - db_file=db_file, + HeteroAnalysisToDb(db_file=db_file, task_label=name, additional_fields={ "tags": spec['tags'] })] - super(HeteroOptimizeFW, self).__init__(t, spec=spec, parents=parents, + super(HeteroOptimizeFW, self).__init__(t, spec=spec, parents=parents, name=name, **kwargs) class SubstrateSlabFW(Firework): @@ -116,41 +117,41 @@ def __init__(self, spec, structure, slab_params, vasp_input_set=None, parents=None, **kwargs): """ Apply the transformations to the bulk structure, write the slab - set corresponding to the transformed structure, and run vasp. Note - that if a transformation yields many structures from one, only the + set corresponding to the transformed structure, and run vasp. Note + that if a transformation yields many structures from one, only the last structure in the list is used. By default all structures will have selective dynamics tags as one of the transformations applied to the input structure. - + Args: name (string): Name for the Firework. - spec (dict): The specification parameters used to control the + spec (dict): The specification parameters used to control the workflow and pass variables. structure (Structure): Bulk input structure to apply transformations onto. - slab_params (dict): A dictionary containing a list of transformations - and transformation_params to generate the substrate slab + slab_params (dict): A dictionary containing a list of transformations + and transformation_params to generate the substrate slab structure for the hetero_interface. Example: slab_params = {'transformations': ['SlabTransformation', 'AddSitePropertyTransformation'],'transformation_params': [{},{}]}. - Definitions: transformations (list): list of names of transformation + Definitions: transformations (list): list of names of transformation classes as defined in the modules in pymatgen.transformations. - transformation_params (list): list of dicts where each dict specify - the input parameters to instantiate the transformation + transformation_params (list): list of dicts where each dict specify + the input parameters to instantiate the transformation class in the transformations list. Other Parameters: user_incar_settings (dict): VASP INCAR settings to override default settings for the vasp input set. - vasp_input_set (VaspInputSet): VASP input set, used to write the - input set for the transmuted structure. Defaults to + vasp_input_set (VaspInputSet): VASP input set, used to write the + input set for the transmuted structure. Defaults to CMDLInterfaceSet. vasp_cmd (string): Command to run vasp. - copy_vasp_outputs (bool): Whether to copy outputs from previous + copy_vasp_outputs (bool): Whether to copy outputs from previous run. Defaults to True. prev_calc_dir (str): Path to a previous calculation to copy from. db_file (string): Path to file specifying db credentials. - parents (Firework): Parents of this particular Firework. FW or + parents (Firework): Parents of this particular Firework. FW or list of FWs. kwargs: Other kwargs that are passed to Firework.__init__. """ @@ -160,7 +161,8 @@ class in the transformations list. # vasp input settings user_incar_settings = user_incar_settings or None vasp_input_set = vasp_input_set or CMDLInterfaceSet(structure, - auto_dipole=True, user_incar_settings=user_incar_settings, + auto_dipole=True, + user_incar_settings=user_incar_settings, vdw='optB88', iface=True) # slab settings @@ -191,39 +193,38 @@ class in the transformations list. raise ValueError("Must specify structure or previous calculation") t.append(RunVaspCustodian(vasp_cmd=vasp_cmd, handler_group=handler)) t.append(PassCalcLocs(name=fw_name)) - t.append(HeteroAnalysisToDb(wf_name=spec['wf_name'], - db_file=db_file, + t.append(HeteroAnalysisToDb(db_file=db_file, task_label=fw_name, additional_fields={ "tags": spec['tags'], "transmuter_input": {"transformations": transformations, "transformation_params": transformation_params} })) - super(SubstrateSlabFW, self).__init__(t, spec=spec, parents=parents, + super(SubstrateSlabFW, self).__init__(t, spec=spec, parents=parents, name=fw_name, **kwargs) class GenHeteroStructuresFW(Firework): def __init__(self, spec, structure, heterotransformation_params, - vasp_input_set=None, user_incar_settings=None, prev_calc_dir=None, + vasp_input_set=None, user_incar_settings=None, prev_calc_dir=None, name="Generate Heterostructures", vasp_cmd=">>vasp_cmd<<", db_file=">>db_file<<", copy_vasp_outputs=False, parents=None, **kwargs): """ - Apply transformations to 2D material and substrate slab to generate + Apply transformations to 2D material and substrate slab to generate fireworks that create and relax hetero_interface structures using - vdW corrections. Note: If the transformation produces many from + vdW corrections. Note: If the transformation produces many from one, all structures are simulated. Args: spec (dict): Specification of the job to run. structure (Structure): 2D material to align onto the substrate. - hetero_transformation_params (dict): dictionary containing the + hetero_transformation_params (dict): dictionary containing the input to control the hetero2d.manipulate.hetero_transmuter modules. transformations (list): list of transformations to apply to the structure as defined in the modules in hetero2d.manipulate module. transformation_params (list): list of dicts where each dict specifies - the input parameters to instantiate the transformation class in the + the input parameters to instantiate the transformation class in the transformations list. Example: h_params={'transformations': ['hetero_interfaces'], 'transformation_params':[{hetero_interface parameters dictionary}]}. @@ -231,9 +232,9 @@ def __init__(self, spec, structure, heterotransformation_params, Other Parameters: user_incar_settings (dict): Input settings to update the settings for the vasp input set. - vasp_input_set (VaspInputSet): VASP input set, used to write the + vasp_input_set (VaspInputSet): VASP input set, used to write the input set for the transmuted structure. Defaults to CMDLInterfaceSet. - name (string): Name for the Firework. Default is "Generate + name (string): Name for the Firework. Default is "Generate HeteroStructures {2D}-on-{Substrate} {unique_id}". vasp_cmd (string): Command to run vasp. copy_vasp_outputs (bool): Whether to copy outputs from the previous run. @@ -243,6 +244,9 @@ def __init__(self, spec, structure, heterotransformation_params, parents (Firework): Parents of this particular Firework. FW or list of FWS. kwargs: Other kwargs that are passed to Firework.__init__. """ + #TODO: @tboland1: add a check to make sure the final structure has a large + # enough periodic z separation between the top layer of the 2d and the bottom layer + # of the substrate. - @tboland1 struct_sub = spec["struct_sub"] # only used to set correct composition # FW Name @@ -254,7 +258,8 @@ def __init__(self, spec, structure, heterotransformation_params, # VASP Input user_incar_settings = user_incar_settings or None vasp_input_set = vasp_input_set or CMDLInterfaceSet(structure, - auto_dipole=True, user_incar_settings=user_incar_settings, + auto_dipole=True, + user_incar_settings=user_incar_settings, vdw='optB88', iface=True) # Create hetero_interface Structure Firetask @@ -262,7 +267,7 @@ def __init__(self, spec, structure, heterotransformation_params, if prev_calc_dir: t.append(CopyVaspOutputs(calc_dir=prev_calc_dir, contcar_to_poscar=True)) # creates hetero_interface structures and then relaxes structures - t.append(CreateHeterostructureTask(spec=spec, + t.append(CreateHeterostructureTask(spec=spec, structure=structure, heterotransformation_params=heterotransformation_params, vasp_input_set=vasp_input_set, @@ -271,7 +276,7 @@ def __init__(self, spec, structure, heterotransformation_params, elif copy_vasp_outputs: t.append(CopyVaspOutputs(calc_loc=True, contcar_to_poscar=True)) # creates hetero_interface structures and then relaxes structures - t.append(CreateHeterostructureTask(spec=spec, + t.append(CreateHeterostructureTask(spec=spec, structure=structure, heterotransformation_params=heterotransformation_params, vasp_input_set=vasp_input_set, @@ -279,7 +284,7 @@ def __init__(self, spec, structure, heterotransformation_params, name=fw_name)) elif structure: # creates hetero_interface structures and then relaxes structures - t.append(CreateHeterostructureTask(spec=spec, + t.append(CreateHeterostructureTask(spec=spec, structure=structure, heterotransformation_params=heterotransformation_params, vasp_input_set=vasp_input_set, @@ -287,7 +292,6 @@ def __init__(self, spec, structure, heterotransformation_params, super(GenHeteroStructuresFW, self).__init__(t, spec=spec, parents=parents, name=fw_name, **kwargs) - class HeteroStructuresFW(Firework): def __init__(self, spec, structure, name=None, transformation=None, vasp_input_set=None, user_incar_settings=None, @@ -295,8 +299,8 @@ def __init__(self, spec, structure, name=None, transformation=None, copy_vasp_outputs=False, prev_calc_dir=None, parents=None, **kwargs): """ - Relax the hetero_structures generated by CreateHeterostructureTask - and perform various energetic analyses to store in the database. + Relax the hetero_structures generated by CreateHeterostructureTask + and perform various energetic analyses to store in the database. Args: spec (dict): The specification dictionary used to control the @@ -308,23 +312,23 @@ def __init__(self, spec, structure, name=None, transformation=None, Other Parameters: transformation (dict): The list of dictionaries containing the input used to create the hetero_interfaces. Defaults to none. - name (str): The name for this firework. Defaults to wf_name + + name (str): The name for this firework. Defaults to wf_name + heterostructure optimization (used as the fw_name). transformation (dict): transformation parameters used to create the heterostructure. user_incar_settings (dict): Input settings to update the settings for the vasp input set. - vasp_input_set (VaspInputSet): VASP input set for the transformed + vasp_input_set (VaspInputSet): VASP input set for the transformed 2d on sub-slab. Default CMDLInterfaceSet. - prev_calc_dir: path to previous calculation if using structure + prev_calc_dir: path to previous calculation if using structure from another calculation. """ fw_name = "{}: {}".format("Heterostructure Optimization", name) # Update VASP input params user_incar_settings = user_incar_settings or None - vasp_input_set = vasp_input_set or CMDLInterfaceSet(structure, - auto_dipole=False, iface=True, vdw='optB88', + vasp_input_set = vasp_input_set or CMDLInterfaceSet(structure, auto_dipole=False, + iface=True, vdw='optB88', user_incar_settings=user_incar_settings) # Analysis Flags for HeteroAnalysisToDb @@ -363,8 +367,7 @@ def __init__(self, spec, structure, name=None, transformation=None, raise ValueError("Must specify structure or previous calculation") t.append(RunVaspCustodian(vasp_cmd=vasp_cmd, handler_group=handler, wall_time=wall_time)) t.append(PassCalcLocs(name=fw_name)) - t.append(HeteroAnalysisToDb(wf_name=spec['wf_name'], - db_file=db_file, + t.append(HeteroAnalysisToDb(db_file=db_file, task_label=fw_name, Formation_Energy=formation_energy, Binding_Energy=binding_energy, @@ -376,3 +379,84 @@ def __init__(self, spec, structure, name=None, transformation=None, super(HeteroStructuresFW, self).__init__(t, spec=spec, parents=parents, name=fw_name, **kwargs) + +class ElectronicFW(Firework): + def __init__(self, structure, name="Electronic", dedos=0.05, grid_density=0.03, + tags={}, dos=True, bader=True, cdd=False, parents=None, prev_calc_dir=None, + vasp_cmd=VASP_CMD, db_file=DB_FILE, copy_vasp_outputs=True, + electronic_set_overrides=None, **kwargs): + """ + Performs standard NonSCF Calculation Firework to generate density of states + and bader files. Proceeds to parse the store results in the database. + + If cdd is set to True, calculations to generate 1D charge density difference + for the given structure are performed. The NonSCFFW for each structure has + the NGiF grid densities set to 0.03 A between each point by default. + + Args: + structure (Structure): Input structure - only used to set the name of + the FW. + + Other Parameters: + name (str): Name for the Firework. Defaults to compsition-Electronic. + dedos (float): Automatically set nedos using the total energy range + which will be divided by the energy step dedos. Default 0.05 eV. + grid_density (float): Distance between grid points for the NGXF,Y,Z + grids. Defaults to 0.03 Angstroms; NGXF,Y,Z are ~2x > default. + tags (dict): A dictionary listing the tags for the firework. + dos (bool): If True, peforms high quality site-orbital projected density + of states. Use electronic_set_overrides to change default params. + Defaults to True. + bader (bool): If True, peforms high quality bader analysis for the given + structure. + cdd (bool): Whether to compute the z projected charge density difference + for the given structure. Default set to False. + parents (Firework): Parents of this particular Firework. FW or list of + FWS. + prev_calc_dir (str): Path to a previous calculation to copy from. + vasp_cmd (str): Command to run vasp. + db_file (str): Path to file specifying db credentials. + copy_vasp_outputs (bool): Whether to copy outputs from previous run. + Defaults to True copying the CHGCAR and OUTCAR. + electronic_set_overrides (dict): Arguments listed in "from_prev_calc" + method of the CMDLElectronicSet that are not explicitly listed + as inputs in this FW. This dictionary allows a user to modify the + default settings of the input set. Valid keys are force_gamma, + small_gap_multiply, nbands_factor, dos_around_fermi, slurm_npar, + auto_dipole, and **kwargs. + **kwargs: Other kwargs that are passed to Firework.__init__. + """ + fw_name = "{}-{}".format(structure.composition.reduced_formula, name) + electronic_set_overrides = electronic_set_overrides or {} + slurm_npar = electronic_set_overrides.get('slurm_npar', False) + + t = [] + if prev_calc_dir: + t.append(CopyVaspOutputs(calc_dir=prev_calc_dir, additional_files=["CHGCAR"])) + elif parents: + t.append(CopyVaspOutputs(calc_loc=True, additional_files=["CHGCAR"])) + else: + raise ValueError("Must specify previous calculation for ElectronicFW") + + t.append(WriteVaspElectronicFromPrev(prev_calc_dir=".", + grid_density=grid_density, + dedos=dedos, + dos=dos, + bader=bader, + cdd=cdd, + **electronic_set_overrides)) + t.append(RunElectronicCustodian(vasp_cmd=vasp_cmd, + handler_group=handler, + slurm_npar=slurm_npar, + auto_npar=">>auto_npar<<")) + t.append(PassCalcLocs(name=name)) + t.append(HeteroAnalysisToDb(db_file=db_file, + task_label=fw_name, + dos=dos, + bader=bader, + cdd=cdd, + additional_fields={ + "task_label": fw_name, + "tags": tags + })) + super(ElectronicFW, self).__init__(t, parents=parents, name=fw_name, **kwargs) diff --git a/hetero2d/io/CMDLRelaxSet.yaml b/hetero2d/io/CMDLRelaxSet.yaml index ffd302d..8b4ef8c 100644 --- a/hetero2d/io/CMDLRelaxSet.yaml +++ b/hetero2d/io/CMDLRelaxSet.yaml @@ -4,7 +4,7 @@ INCAR: ALGO: FAST ENCUT: 520 - EDIFF: 1e-6 + EDIFF: 1e-4 EDIFFG: -0.02 IBRION: 2 ICHARG: 1 diff --git a/hetero2d/io/VaspInterfaceSet.py b/hetero2d/io/VaspInterfaceSet.py index 8a1de18..dece226 100644 --- a/hetero2d/io/VaspInterfaceSet.py +++ b/hetero2d/io/VaspInterfaceSet.py @@ -5,17 +5,18 @@ """ This class instantiates the default VASP input set to simulate 2D-substrate slab hetero_interfaces. """ - from __future__ import division, unicode_literals, print_function -import numpy as np -import os -import warnings +import os, warnings, six, numpy as np +from math import ceil, floor from monty.json import MSONable from monty.serialization import loadfn + from pymatgen import Structure -from pymatgen.io.vasp.sets import DictSet +from pymatgen.io.vasp.inputs import Incar, Kpoints, Poscar, Potcar +from pymatgen.io.vasp.sets import DictSet, get_vasprun_outcar + __author__ = "Tara M. Boland" __copyright__ = "Copyright 2020, CMD Lab" @@ -24,16 +25,14 @@ __date__ = "Jan 06, 2020" MODULE_DIR = os.path.dirname(os.path.abspath(__file__)) -warnings.filterwarnings("ignore") +warnings.filterwarnings("ignore") def _load_yaml_config(fname): config = loadfn(os.path.join(MODULE_DIR, "%s.yaml" % fname)) - config["INCAR"].update(loadfn(os.path.join(MODULE_DIR, - "VASPIncarBase.yaml"))) + config["INCAR"].update(loadfn(os.path.join(MODULE_DIR, "VASPIncarBase.yaml"))) return config - class CMDLRelaxSet(DictSet): """ Implementation of VaspInputSet utilizing parameters in the public 2D Materials Synthesis Database. @@ -41,24 +40,22 @@ class CMDLRelaxSet(DictSet): CONFIG = _load_yaml_config("CMDLRelaxSet") def __init__(self, structure, **kwargs): - super(CMDLRelaxSet, self).__init__( - structure, CMDLRelaxSet.CONFIG, **kwargs) + super(CMDLRelaxSet, self).__init__(structure, CMDLRelaxSet.CONFIG, **kwargs) self.kwargs = kwargs - class CMDLInterfaceSet(CMDLRelaxSet): """ - Class for writing a set of interface vasp runs, including the 2D, - substrate slabs, and hetero_interface (oriented along the - c direction) as well as the bulk (3D) phase of the 2D material and the - un-oriented unit cells, to ensure the same K_POINTS, POTCAR and INCAR criterion. + Class for writing a set of interface vasp runs, including the 2D, substrate + slabs, and hetero_interface (oriented along the c direction) as well as the + bulk (3D) phase of the 2D material and the un-oriented unit cells, to ensure + the same K_POINTS, POTCAR and INCAR criterion. Args: - structure (Structure): The structure for the calculation + structure (Structure): The structure for the calculation. Other Parameters: - k_product (int): default to 20, k_point number x length for a & b directions, - also for c direction in bulk calculations + k_product (int): Default to 20, k_point number x length for a & b directions, + also for c direction in bulk calculations. iface (bool): Set to False for initial structure calculations to fully relax the structures to calculate the adsorption energy of a 2D material on the substrate slab. Defaults to True. @@ -72,18 +69,15 @@ class CMDLInterfaceSet(CMDLRelaxSet): False. user_incar_settings (dict): A way to override the default settings for the incar. - kwargs: - Other kwargs supported by :class:`DictSet`. + kwargs: Other kwargs supported by :class:`DictSet`. Returns: Vasp input set object """ - def __init__(self, structure, k_product=20, iface=True, vdw='optB88', auto_dipole=False, set_mix=False, sort_structure=False, - user_incar_settings=None, **kwargs): + user_incar_settings={}, **kwargs): super(CMDLInterfaceSet, self).__init__(structure, **kwargs) - if sort_structure: structure = structure.get_sorted_structure() @@ -92,9 +86,9 @@ def __init__(self, structure, k_product=20, iface=True, vdw='optB88', self.iface = iface self.vdw = vdw.lower() if vdw is not None else None self.auto_dipole = auto_dipole - self.kwargs = kwargs self.set_mix = set_mix self.user_incar_settings = user_incar_settings + self.kwargs = kwargs iface_incar = {} if self.iface: @@ -112,6 +106,7 @@ def __init__(self, structure, k_product=20, iface=True, vdw='optB88', iface_incar["IDIPOL"] = 3 iface_incar["LDIPOL"] = True iface_incar["DIPOL"] = ' '.join(str(i) for i in center_of_mass) + if self.vdw: vdw_par = loadfn(os.path.join(MODULE_DIR, "vdW_parameters.yaml")) try: @@ -123,24 +118,19 @@ def __init__(self, structure, k_product=20, iface=True, vdw='optB88', if user_incar_settings: iface_incar.update(user_incar_settings) + + # modify configuration dictionary self._config_dict["INCAR"].update(iface_incar) - if self._config_dict["INCAR"]['ISPIN'] == 1: self._config_dict["INCAR"].pop("MAGMOM", '') - + @property def kpoints(self): """ - k_product (default to 20) is k_point number * length for a & b - directions, also for c direction in bulk calculations. Automatic - mesh & Gamma is the default setting. + k_product (default to 20) is the number of k-points * length for a & b + directions, also for c direction in bulk calculations. Results in + k_product k-points/Angstrom. Defaults to automatic mesh & Gamma. """ - - # To get input sets, the input structure has to has the same number - # of required parameters as a Structure object (ie. 4). Slab - # attributes aren't going to affect the VASP inputs anyways so - # converting the slab into a structure should not matter - kpt = super(CMDLInterfaceSet, self).kpoints kpt.comment = "Automatic mesh" kpt.style = 'Gamma' @@ -158,12 +148,171 @@ def kpoints(self): return kpt - def as_dict(self, verbosity=2): +class CMDLElectronicSet(CMDLRelaxSet): + """ + Class for writing vasp inputs for DOS, bader, and charge density difference + runs from previous VASP jobs. Typically, you would use the classmethod + from_prev_calc to initialize from a previous SCF run. + + Args: + structure (Structure): The structure for the calculation. + + Other Parameters: + reciprocal_density (int): For static calculations, we usually set the + reciprocal density by volume. This is a convenience arg to change + that, rather than using user_kpoints_settings. Defaults to 100, + which is ~50% more than that of standard relaxation calculations. + force_gamma (bool): Force gamma k-point mesh. + prev_incar (Incar/string): Incar file from previous run. + **kwargs: kwargs supported by CMDLRelaxSet. + + Returns: + Vasp input set object + """ + def __init__(self, structure, reciprocal_density=200, force_gamma=False, + prev_incar=None, **kwargs): + super(CMDLElectronicSet, self).__init__(structure, **kwargs) + if isinstance(prev_incar, six.string_types): + prev_incar = Incar.from_file(prev_incar) + + self.reciprocal_density = reciprocal_density + self.force_gamma = force_gamma + self.prev_incar = prev_incar + self.kwargs = kwargs + + @property + def incar(self): """ - A JSON serializable dict representation of an object. + :return: Incar """ + parent_incar = super(CMDLElectronicSet, self).incar # get inputset parent incar + # if no prev use parent + incar = Incar(self.prev_incar) if self.prev_incar is not None else Incar(parent_incar) + + # remove incar tags that are not needed for NonSCF + remove = ["EDIFFG", "ISIF", "LREAL", "POTIM", "KPOINT_BSE", "NELMDL", "MAGMOM", + "AMIX_MAG", "BMIX_MAG", "AMIX", "BMIX", "IMIX", "AMIN", "NELMIN"] + [ incar.pop(tag) for tag in remove if tag in incar.keys() ] + + # enforce tags for NonSCF calcs: ICHARG should be 11 + incar.update({"IBRION": -1, "ISTART": 1, "LWAVE": False, "NSW": 0, "ISYM": 0, + "ICHARG": 11}) + incar.update(self.kwargs.get("user_incar_settings", {})) + return incar + + @property + def kpoints(self): + """ + Generate a dense k-point grid for dos/bader analysis. + """ + return Kpoints.automatic_density_by_vol(self.structure, self.reciprocal_density, + force_gamma=self.force_gamma) + + @classmethod + def from_prev_calc(cls, prev_calc_dir, dedos=0.05, grid_density=0.03, + dos=True, bader=True, cdd=False, small_gap_multiply=None, + nbands_factor=1, dos_around_fermi=[4,6], auto_dipole=False, + **kwargs): + """ + Generate a set of Vasp input files for ElectronicFW calculations from a + directory of previous directory. + + Args: + prev_calc_dir (str): The directory containing the outputs( + vasprun.xml and OUTCAR) from the previous vasp run. + + Other Parameters: + dedos (float): Automatically set nedos using the total energy range + which will be divided by the energy step dedos. Default 0.05 eV. + grid_density (float): Distance between grid points for the NGXF,Y,Z grids. + Defaults to 0.03 Angs; NGXF,Y,Z are ~2x > default. For charge + density difference calculations the parent grid density is used for all + children fireworks. + dos (bool): If True, sets INCAR tags for high quality site-orbital projected + density of states. Defaults to True. + bader (bool): If True, sets INCAR tags to generate bader analysis files. + cdd (bool): If True, ensures the grid density matches between the parent + and child Fireworks. Default set to False. + small_gap_multiply ([float, float]): If the gap is less than + 1st index, multiply the default reciprocal_density by the 2nd + index. + nbands_factor (float): Multiplicative factor for NBANDS. + dos_around_fermi (bool/list): The element projected density of states is + calculated around the fermi level. Default range is [efermi-4, efermi+6]. + If you want a different range supply a list i.e. [4,6] or False to compute + the entire dos range. + auto_dipole (bool): Whether to set dipole corrections. Defaults to + False. + **kwargs: All kwargs supported by CMDLRelaxSet, other than prev_incar + and structure which are determined from the prev_calc_dir. + """ + vasprun, outcar = get_vasprun_outcar(prev_calc_dir) + + incar = vasprun.incar + structure = vasprun.final_structure + + # Turn off spin when magmom for every site is smaller than 0.02. + if outcar and outcar.magnetization: + site_magmom = np.array([i["tot"] for i in outcar.magnetization]) + ispin = 2 if np.any(site_magmom[np.abs(site_magmom) > 0.02]) else 1 + elif vasprun.is_spin: + ispin = 2 + else: + ispin = 1 + + # set nbands factor + nbands = int(np.ceil(vasprun.parameters["NBANDS"] * nbands_factor)) + incar.update({"ISPIN": ispin, "NBANDS": nbands}) + + # set dipole corrections + if auto_dipole: + weights = [s.species_and_occu.weight for s in structure] + center_of_mass = np.average(structure.frac_coords, + weights=weights, axis=0) + incar["IDIPOL"] = 3 + incar["LDIPOL"] = True + incar["DIPOL"] = ' '.join(str(i) for i in center_of_mass) + + # multiply the reciprocal density if needed: + if small_gap_multiply: + gap = vasprun.eigenvalue_band_properties[0] + if gap <= small_gap_multiply[0]: + kwargs["reciprocal_density"] = kwargs.get("reciprocal_density", 200) * small_gap_multiply[1] + + # check if the previous calc was metallic or insulating + gap = vasprun.complete_dos.get_gap() + if gap > 0.0: # if insulator/semi -5 + incar.update({"ISMEAR":-5}) + incar.pop("SIGMA") + else: # ismear 0; 0.05 smear + incar.update({"ISMEAR": 0, "SIGMA": 0.05}) + + # dos settings + if dos: + # automatic setting of nedos using the energy range and the energy step dedos + if dos_around_fermi: + fermi = vasprun.efermi + emin, emax = floor(fermi - dos_around_fermi[0]), ceil(fermi + dos_around_fermi[1]) + nedos = ceil(abs(emin - emax) / dedos) # compute dos spacing + incar.update({"LORBIT": 11, "NEDOS": nedos + 1 if nedos % 2 == 0 else nedos, + "EMIN": emin, "EMAX": emax}) + else: + emax, emin = max(vasprun.complete_dos.energies), min(vasprun.complete_dos.energies) + nedos = ceil(abs(emin - emax) / dedos) # compute dos spacing + incar.update({"LORBIT": 11, "NEDOS": nedos + 1 if nedos % 2 == 0 else nedos}) + # bader analysis settings + if bader: + # get new NGiF grid spacing + a, b, c = structure.lattice.abc + ngxf, ngyf, ngzf = [int(ceil(i/grid_density/10) * 10) for i in [a,b,c]] + + # modification to the incar + bader_settings = {"LCHARG": True, "LAECHG": True, "NGXF": ngxf, "NGYF": ngyf, + "NGZF": ngzf} + incar.update(bader_settings) + + # charge density difference settings + if cdd: + cdd_settings = {} - d = MSONable.as_dict(self) - if verbosity == 1: - d.pop("structure", None) - return d + return cls(structure=structure, prev_incar=incar, **kwargs) diff --git a/hetero2d/io/__init__.py b/hetero2d/io/__init__.py index 04d77e7..45ede1b 100644 --- a/hetero2d/io/__init__.py +++ b/hetero2d/io/__init__.py @@ -1 +1,2 @@ __author__ = 'Tara M. Boland ' +from .VaspInterfaceSet import CMDLElectronicSet, CMDLInterfaceSet diff --git a/hetero2d/manipulate/heterotransmuter.py b/hetero2d/manipulate/heterotransmuter.py index 855eadd..c5ef336 100644 --- a/hetero2d/manipulate/heterotransmuter.py +++ b/hetero2d/manipulate/heterotransmuter.py @@ -8,19 +8,19 @@ atomate architecture, to fix minor bugs in the original code, and return interface matching criteria. """ +import numpy as np, sys from six.moves import range from copy import deepcopy -import numpy as np, sys from mpinterfaces.transformations import reduced_supercell_vectors, get_r_list, get_angle, \ - get_mismatch, get_area, get_uniq_layercoords + get_mismatch, get_area from pymatgen import Structure, Lattice from pymatgen.analysis.structure_matcher import StructureMatcher from pymatgen.analysis.structure_analyzer import SpacegroupAnalyzer from pymatgen.core.surface import Slab -from hetero2d.manipulate.utils import slab_from_struct, show_struct_ase, get_fu +from hetero2d.manipulate.utils import slab_from_struct, get_fu from hetero2d.manipulate.layersolver import LayerSolver __author__ = "Tara M. Boland, Arunima Singh" @@ -70,10 +70,9 @@ def rotate_to_acute(structure): return new_structure - def aligned_hetero_structures(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, nlayers_2d=3, nlayers_sub=2, r1r2_tol=0.02, max_angle_diff=1, - separation=3.4): + separation=3.4, symprec=False): """ Given the 2 slab structures and the alignment parameters, return slab structures with lattices that are aligned with respect to each @@ -99,7 +98,11 @@ def aligned_hetero_structures(struct_2d, struct_sub, max_mismatch=0.01, max_area 3 layers. separation (float): separation between the substrate and the 2d material. Defaults to 3.4 angstroms. - + symprec (bool/float): Perform symmetry matching to the specified + tolerance for symmetry finding between the aligned and + original 2D structure. Enable if you notice the 2D lattice is + disformed. Defaults to False. + Returns: aligned_sub, aligned_2d, alignment_info """ @@ -285,25 +288,27 @@ def aligned_hetero_structures(struct_2d, struct_sub, max_mismatch=0.01, max_area mat2d.modify_lattice(lmap) ## ensure that the symmetry of the 2d and aligned 2d agree - sg_mat2d = SpacegroupAnalyzer(mat2d) - sg_align = SpacegroupAnalyzer(struct_2d) + if symprec: + sg_mat2d = SpacegroupAnalyzer(mat2d, symprec=symprec) + sg_align = SpacegroupAnalyzer(struct_2d, symprec=symprec) - m2d = {} - align = {} - [m2d.update({key:sg_mat2d.get_symmetry_dataset()[key]}) - for key in ['number','hall_number','international','hall','pointgroup']] - [align.update({key:sg_align.get_symmetry_dataset()[key]}) - for key in ['number','hall_number','international','hall','pointgroup']] - - # compare mat2d and aligned lattice - is_equal = m2d == align - - if is_equal == True: - return substrate, mat2d, alignment_info + m2d = {} + align = {} + [m2d.update({key:sg_mat2d.get_symmetry_dataset()[key]}) + for key in ['number','hall_number','international','hall','pointgroup']] + [align.update({key:sg_align.get_symmetry_dataset()[key]}) + for key in ['number','hall_number','international','hall','pointgroup']] + # compare mat2d and aligned lattice + is_equal = m2d == align + + if is_equal == True: + return substrate, mat2d, alignment_info + else: + print('SpacegroupAnalyzer failed.\n') + return None, None, None + else: - print('SpacegroupAnalyzer failed.\n') - return None, None, None - + return substrate, mat2d, alignment_info def generate_all_configs(mat2d, substrate, nlayers_2d=3, nlayers_substrate=2, separation=3.4): """ @@ -311,7 +316,9 @@ def generate_all_configs(mat2d, substrate, nlayers_2d=3, nlayers_substrate=2, se unique (Wyckoff) sites of the mat2d and substrate. The unique sites are iteratively matched between the mat2d and substrate stacking the unique sites on top of each other separated by the separation distance parameter. This subsequently generates all possible 2d/substrate heterostructure - configurations stacked over high symmetry points. All unique structures are returned. + configurations stacked over high symmetry points. All unique structures are returned. To identify + the wyckoff sites stacked on top of each other the name is attached to each structure under + site_properties. Args: mat2d (Structure): Lattice and symmetry-matched 2D material @@ -332,10 +339,16 @@ def generate_all_configs(mat2d, substrate, nlayers_2d=3, nlayers_substrate=2, se sys.exit() # unique site coordinates in the substrate top layers - coords_uniq_sub = get_uniq_layercoords(substrate, nlayers_substrate, top=True) + coords_uniq_sub, uniq_sub_idx = get_uniq_layercoords(substrate, nlayers_substrate, top=True) + # add the wyckoff site to each site in substrate + spg = SpacegroupAnalyzer(substrate) + substrate.add_site_property('wyckoffs', spg.get_symmetry_dataset()['wyckoffs']) # unique site coordinates in the 2D material bottom layers - coords_uniq_2d = get_uniq_layercoords(mat2d, nlayers_2d, top=False) + coords_uniq_2d, uniq_2d_idx = get_uniq_layercoords(mat2d, nlayers_2d, top=False) + # add the wyckoff site to each site in mat2d + spg = SpacegroupAnalyzer(mat2d) + mat2d.add_site_property('wyckoffs', spg.get_symmetry_dataset()['wyckoffs']) # set separation distance betwn 2d and substrate substrate_top_z = np.max(np.array([site.coords for site in substrate])[:, 2]) @@ -351,12 +364,16 @@ def generate_all_configs(mat2d, substrate, nlayers_2d=3, nlayers_substrate=2, se # an interface structure for each parallel shift # interface = 2D material + substrate hetero_interfaces = [] - shifted_mat2ds = [] - for coord_i in coords_uniq_sub: - for coord_j in coords_uniq_2d: - interface = substrate.copy() # in new code - new_2d = deepcopy(mat2d) # shfited 2d w/o substrate - new_2d.remove_sites(range(0, new_2d.num_sites)) + for coord_i, sub_idx in zip(coords_uniq_sub, uniq_sub_idx): + for coord_j, td_idx in zip(coords_uniq_2d, uniq_2d_idx): + # create interface structure to append the shifted 2d material + interface = substrate.copy() + interface.remove_site_property('wyckoffs') # remove wyckoffs + + # get the wyckoff site of the 2 sites being aligned on top of each other + td_wyckoff = mat2d[td_idx].properties['wyckoffs'] + sub_wyckoff = substrate[sub_idx].properties['wyckoffs'] + name='2D-'+td_wyckoff+' '+'Sub-'+sub_wyckoff # stacking of 2d on sub wyckoff # set the x,y coors for 2d shift_parallel = coord_i - coord_j @@ -371,21 +388,14 @@ def generate_all_configs(mat2d, substrate, nlayers_2d=3, nlayers_substrate=2, se new_coords = new_coords + origin + shift_net interface.append(site.specie, new_coords, coords_are_cartesian=True) - new_2d.append(site.specie, new_coords, - coords_are_cartesian=True) # add for name + interface.add_site_property('name', [name] * interface.num_sites) hetero_interfaces.append(interface) - shifted_mat2ds.append(new_2d) # shifted pure 2d - - # get the names for the hetero_interfaces - for td, iface in zip(shifted_mat2ds, hetero_interfaces): - name = iface_name(td, substrate) - iface.add_site_property('name', [name] * iface.num_sites) return hetero_interfaces - -def hetero_interfaces(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, - nlayers_2d=3, nlayers_sub=2, r1r2_tol=0.02, max_angle_diff=1, separation=3.4): +def hetero_interfaces(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, nlayers_2d=3, + nlayers_sub=2, r1r2_tol=0.02, max_angle_diff=1, separation=3.4, + symprec=False): """ The given the 2D material and the substrate slab, the 2 slabs are combined to generate all possible unique structures for the @@ -418,7 +428,11 @@ def hetero_interfaces(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, to 2 layers with bottom layer frozen. separation (float): Separation distance between struct_sub and struct_2d. Default 3.4 Angstroms. - + symprec (bool/float): Perform symmetry matching to the specified + tolerance for symmetry finding between the aligned and + original 2D structure. Enable if you notice the 2D lattice is + disformed. Defaults to False. + Returns: Unique hetero_structures list, last entry contains lattice alignment information. Site properties contain iface name. @@ -432,7 +446,8 @@ def hetero_interfaces(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, max_area=max_area, max_mismatch=max_mismatch, max_angle_diff=max_angle_diff, - r1r2_tol=r1r2_tol) + r1r2_tol=r1r2_tol, + symprec=symprec) # exit if the aligned_hetero_structures returns None due to bad symmetry if sub_aligned is None: @@ -449,7 +464,8 @@ def hetero_interfaces(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, max_area=max_area, max_mismatch=max_mismatch, max_angle_diff=max_angle_diff, - r1r2_tol=r1r2_tol) + r1r2_tol=r1r2_tol, + symprec=symprec) # exit if the aligned_hetero_structures returns None due to bad symmetry if sub_aligned is None: @@ -477,50 +493,43 @@ def hetero_interfaces(struct_2d, struct_sub, max_mismatch=0.01, max_area=200, hetero_interfaces.append(alignment_info) return hetero_interfaces - -def iface_name(mat2d, substrate): - """ - Helper function used to generate a unique interface name for a set of interfaces. The substrate's name - will always be the same but the 2d materials indices shit and these are changed from one configuration to - the next. +def get_uniq_layercoords(struct, nlayers, top=True): """ - # Gather Layer Data from LayerSolver - sub_data = LayerSolver(structure=substrate) - td_data = LayerSolver(structure=mat2d) - - # number of 2d layers to get the bottom layer - numl_2d = td_data['num_layers'] + Returns the coordinates and indices of unique sites in the top or bottom + nlayers of the given structure. - # top layer of the substrate - top_layer = sub_data['Layer0'] - # bottom layer of the mat2d - bottom_layer = td_data['Layer' + str(numl_2d - 1)] - - # Generate Substrate Name - wyc_s = np.array(top_layer['wyckoffs']) # get unique wyckoffs - wyc_lyr_s, idx_s, multi_s = np.unique(wyc_s, return_index=True, - return_counts=True) # get unique layer info - - # pull the unique atom species strings - s_s = [top_layer['sites'][idx].species_string - for idx in idx_s] - - # create the substrates name - sub_name = ','.join([str(s) + '(' + str(w) + ')' - for s, w in zip(s_s, wyc_lyr_s)]) - - # Generate Two-D Name - wyc_t = np.array(bottom_layer['wyckoffs']) - wyc_lyr_t, idx_t, multi_t = np.unique(wyc_t, return_index=True, - return_counts=True) - s_t = [bottom_layer['sites'][idx].species_string - for idx in idx_t] - td_name = ';'.join([str(s) + '(' + str(w) + ')' - for s, w in zip(s_t, wyc_lyr_t)]) - - name = sub_name + td_name - return name + Args: + struct (Structure): Input structure. + nlayers (int): Number of layers. + top (bool): Top or bottom layers, default is top layer. + Returns: + numpy array of unique coordinates, corresponding list of atom indices + """ + coords = np.array([site.coords for site in struct]) + z = coords[:, 2] + z = np.around(z, decimals=4) + zu, zuind = np.unique(z, return_index=True) + z_nthlayer = z[zuind[-nlayers]] + zfilter = (z >= z_nthlayer) + if not top: + z_nthlayer = z[zuind[nlayers - 1]] + zfilter = (z <= z_nthlayer) + # site indices in the layers + indices_layers = np.argwhere(zfilter).ravel() + sa = SpacegroupAnalyzer(struct) + symm_data = sa.get_symmetry_dataset() + # equivalency mapping for the structure + # i'th site in the struct equivalent to eq_struct[i]'th site + eq_struct = symm_data["equivalent_atoms"] + # equivalency mapping for the layers + eq_layers = eq_struct[indices_layers] + # site indices of unique atoms in the layers + __, ueq_layers_indices = np.unique(eq_layers, return_index=True) + # print(ueq_layers_indices) + indices_uniq = indices_layers[ueq_layers_indices] + # coordinates of the unique atoms in the layers + return coords[indices_uniq], indices_uniq def remove_duplicates(uv_list, tm_list): """ diff --git a/hetero2d/manipulate/layersolver.py b/hetero2d/manipulate/layersolver.py index aa332b8..a465b2d 100644 --- a/hetero2d/manipulate/layersolver.py +++ b/hetero2d/manipulate/layersolver.py @@ -6,17 +6,17 @@ Analyze an existing structure and parse it based on atomic layers. The cutoff can be changed but defaults to 0.5 Angstroms. Returns a LayeredStructure dictionary. """ - +import numpy as np from copy import deepcopy from operator import itemgetter -import numpy as np from pymatgen.core.structure import Structure from pymatgen.symmetry.analyzer import SpacegroupAnalyzer + __author__ = "Sydney Olson" __copyright__ = "Copyright 2020, CMD Lab" -__maintainer__ = "Sydney Olson, Tara M. Boland" +__maintainer__ = "Tara M. Boland" __email__ = "snolson1@asu.edu, tboland1@asu.edu" __date__ = "June 5, 2020" diff --git a/hetero2d/manipulate/utils.py b/hetero2d/manipulate/utils.py index eaddb07..74a96fc 100644 --- a/hetero2d/manipulate/utils.py +++ b/hetero2d/manipulate/utils.py @@ -7,29 +7,36 @@ the user. """ +import re, shlex, os, sys, json, gzip, numpy as np, math, pymongo, gridfs, zlib from copy import deepcopy -from monty.serialization import loadfn -import numpy as np, math, pymongo +from bson import ObjectId +from monty.serialization import loadfn, dumpfn -from pymatgen import Structure, Lattice, Specie +from ase import Atom +from ase.visualize import view +from ase.calculators.vasp import VaspChargeDensity + +from pymatgen import Structure, Lattice, Specie, Element +from pymatgen.io.ase import AseAtomsAdaptor from pymatgen.core.surface import Slab from pymatgen.analysis.local_env import CrystalNN from pymatgen.analysis.structure_matcher import StructureMatcher from pymatgen.symmetry.analyzer import SpacegroupAnalyzer -from pymatgen.io.ase import AseAtomsAdaptor +from pymatgen.electronic_structure.dos import Dos,CompleteDos -from ase.visualize import view -from ase import Atom +from atomate.vasp.powerups import get_fws_and_tasks from hetero2d.manipulate.layersolver import LayerSolver + __author__ = "Tara M. Boland, Arunima Singh" __copyright__ = "Copyright 2020, CMD Lab" __maintainer__ = "Tara M. Boland" __email__ = "tboland1@asu.edu" __date__ = "June 5, 2020" -# Connection utilities +############################ +### Connection utilities ### def get_mongo_client(db_file, db_type=None): """ Connect to the database using mongoclient. Has multiple @@ -64,10 +71,133 @@ def get_mongo_client(db_file, db_type=None): + database_name + "?authSource=admin" + "&retryWrites=true" + "&w=majority" \ + "&wtimeout=300000" client = pymongo.MongoClient(connection_url) - return client, database_name -# Post processing functions +def get_dos(client, taskdoc): + ''' + Uses the object_id from the task document to search for the dos in gridfs + and loads the DOS object into a PMG DOS object. The key get_dos must be + true in the task doc or the object does not exist in the gridfs. + + Args: + client (MongoClient Object): The MongoClient connected to the database + you are extracting information from to get the dos. + taskdoc (dict): The task doc you want to get the DOS for. + + Returns: + DOS obect + ''' + # TASKDOC LOOKUP: DosBader._id=dos_fs.files._id=dos_fs.chunks.files_id + fs_id = ObjectId(taskdoc['_id']) + + # connect to dos_fs + fs = gridfs.GridFS(client, 'dos_fs') + bs_json = zlib.decompress(fs.get(fs_id).read()) + obj_dict = json.loads(bs_json.decode()) + + try: + dos_obj = CompleteDos.from_dict(obj_dict) + except: + dos_obj = { Element(key): Dos.from_dict(value) for key, value in obj_dict.items()} + + return dos_obj + +############################################ +### Post/Pre processing helper functions ### +def vtotav(chgcar_file, axis='z'): + ''' + A script which averages a CHGCAR or LOCPOT file in one axis to + make a 1D curve. User must specify filename and axis on command + line. Depends on ase (converted to 3.6 from ase's vtotav.py + python 2.7). + + Args: + chgcar (str): Path to the VASP CHGCAR file + axis (str): Axis to project the charge density onto. + + Returns: + { 'distance' (Ang.), 'chg_density' (projected e/A^3) } + ''' + if not os.path.isfile(chgcar_file): + print("\n** ERROR: Input file %s was not found." % chgcar_file) + sys.exit() + + # Specify the axis to make average in. Default is Z. + allowed = "xyzXYZ" + if allowed.find(axis) == -1 or len(axis) !=1: + axis = 'Z' if allowed.find(axis) == -1 or len(axis)!=1 else axis + axis = axis.upper() if axis.islower() else axis + + # Open geometry and density class objects + #----------------------------------------- + vasp_charge = VaspChargeDensity(filename = chgcar_file) + potl = vasp_charge.chg[-1] + atoms = vasp_charge.atoms[-1] + del vasp_charge + + # Read in lattice parameters and scale factor + #--------------------------------------------- + cell = atoms.cell + + # Find length of lattice vectors + #-------------------------------- + latticelength = np.dot(cell, cell.T).diagonal() + latticelength = latticelength**0.5 + + # Read in potential data + #------------------------ + ngridpts = np.array(potl.shape) + x_grid ,y_grid ,z_grid = ngridpts[0],ngridpts[1],ngridpts[2] + totgridpts = ngridpts.prod() # Total number of points + sys.stdout.flush() + + # Perform average + #----------------- + if axis=="X": + idir = 0 + a = 1 + b = 2 + elif axis=="Y": + a = 0 + idir = 1 + b = 2 + else: + a = 0 + b = 1 + idir = 2 + a = (idir+1)%3 + b = (idir+2)%3 + + # At each point, sum over other two indices + average = np.zeros(ngridpts[idir],np.float) + for ipt in range(ngridpts[idir]): + if axis=="X": + average[ipt] = potl[ipt,:,:].sum() + elif axis=="Y": + average[ipt] = potl[:,ipt,:].sum() + else: + average[ipt] = potl[:,:,ipt].sum() + + # Scale chgcar by size of area element in the plane, + # gives unit e/Ang. I.e. integrating the resulting + # CHG_dir file should give the total charge. + area = np.linalg.det([ (cell[a,a], cell[a,b] ), + (cell[b,a], cell[b,b])]) + dA = area/(ngridpts[a]*ngridpts[b]) + average *= dA + + # Get the average + #----------------- + sys.stdout.flush() + proj_chg = {'distance': [], 'chg_density': []} + xdiff = latticelength[idir]/float(ngridpts[idir]-1) + for i in range(ngridpts[idir]): + x = i*xdiff + proj_chg['distance'].append(round(x, 8)) + proj_chg['chg_density'].append(round(average[i], 8)) + return proj_chg + +## Structure Analyzers ## def average_z_sep(structure, iface_idx, initial=None): ''' Get the average z separation distance between 2 @@ -92,11 +222,11 @@ def average_z_sep(structure, iface_idx, initial=None): ''' # get the 2d sub separation dist (bot & top atoms) # define layers for top and bottom - nlayers = iface_idx['nlayers_2d'] + nlayers = iface_idx['num_2d_layer'] z = ['sub_layer_1','2d_layer_1','2d_layer_'+str(nlayers)] # get the indices from iface that define the interface - idx_subs, idx_tD_top, idx_tD = [iface_idx[z] for i in z] + idx_subs, idx_tD_top, idx_tD = [iface_idx[i] for i in z] coords_f = structure.cart_coords # compute the initial 2d film width @@ -126,23 +256,6 @@ def average_z_sep(structure, iface_idx, initial=None): else: return zsep_f, delta_2d_z -def center_slab(structure): - """ - Centers the atoms in a slab structure around 0.5 - fractional height. - - Args: - structure (Structure): Structure to center. - - Returns: - Centered Structure object - """ - - center = np.average([s._fcoords[2] for s in structure.sites]) - translation = (0, 0, 0.5 - center) - structure.translate_sites(range(len(structure.sites)), translation) - return structure - def get_fu(struct_sub, struct_2d, sub_aligned, td_aligned): ''' Given superlattice structures and original unit cell structures @@ -194,6 +307,127 @@ def get_fu(struct_sub, struct_2d, sub_aligned, td_aligned): return n_2d, n_sub +def center_slab(structure): + """ + Centers the atoms in a slab structure around 0.5 + fractional height. + + Args: + structure (Structure): Structure to center. + + Returns: + Centered Structure object + """ + + center = np.average([s._fcoords[2] for s in structure.sites]) + translation = (0, 0, 0.5 - center) + structure.translate_sites(range(len(structure.sites)), translation) + return structure + +def tag_iface(structure, nlayers_2d, nlayers_sub=2): + """ + Find the atom indices in a heterostructure by specifying how many + layers there are for the 2D material. Returns a dictionary of atom + ids for each layer of the 2D material and nlayers_sub of the + substrate surface. + + Args: + structure (Structure): The hetero_interface structure. + Should be an unrelaxed structure. + nlayers_2d (int): The number of layers contained + within the 2d material. + nalyers_sub (int): The number of layers of the substrate + surface to include in the interface tags. + + Returns: + layer_indices + """ + # set structure as slab + interface = slab_from_struct(structure) + # get z coords of the sites + z_coords = interface.frac_coords[:, 2] + # round z coords to avoid precision errors from np + z_coords_round = [round_decimals_down(i, decimals=3) for i in z_coords] + # get the z-cords of atoms at interface + unique_layer_coords = np.unique(z_coords_round) + + # get all unique coords for each layer + layer_coords_dict = {} + # get 2d layer indices + for idx in list(range(1, nlayers_2d+1)): + layer_coords = unique_layer_coords[-idx] + layer_name = '2d_layer_' + str(abs(idx)) + layer_coords_dict[layer_name] = layer_coords + # get sub layer indices + for count, idx in enumerate(range(nlayers_2d+1, nlayers_2d+nlayers_sub+1)): + layer_coords = unique_layer_coords[-idx] + layer_name = 'sub_layer_' + str(count+1) + layer_coords_dict[layer_name] = layer_coords + + # get atom ids for each layer + layer_indices = {key: [] for key in layer_coords_dict.keys()} + for idx, site in enumerate(structure.sites): + site_coords = round_decimals_down(site.frac_coords[2], decimals=3) + if site_coords in layer_coords_dict.values(): + key = get_key(layer_coords_dict, site_coords) + layer_indices[key].append(idx) + + layer_indices['num_2d_layer'] = nlayers_2d + layer_indices['num_sub_layer'] = nlayers_sub + + return layer_indices + +def iface_layer_locator(structure, cutoff, iface_elements): + """ + Helper function used to locate the iface layers from + the layer solver. + + Args: + structure (Structure): input structure to find the + interface layers for. + cutoff (float): the layer separation cutoff distance. + iface_elements (list): A list of element species + that compose the top and bottom iface layers. + + Returns: + LayerSolver object, top_layer specie, bottom_layer specie + """ + sample = LayerSolver(structure, cutoff=cutoff) + flag = True + + # get the 2d sub separation dist (bot & top atoms) + for i in range(0, len(sample) - 2): + # define layer data + l1 = sample['Layer' + str(i + 1)] + l2 = sample['Layer' + str(i)] + + # sort layers by species + test_layer = sorted([l1['species'][0], + l2['species'][0]]) + + # make sure test_layer matches the iface composition + if test_layer == sorted(iface_elements): + return sample, l1, l2 + + if flag: + print('\t\t## Iface composition not found. Maybe the surface has undergone significant relaxation.') + return None, None, None + +def atomic_distance(structure, dist=2): + """ + Given a structure and an interatomic separation distance all sites which have a + spacing less than dist will be printed out. + + Args: + structure (Structure): The structure to analyze. + dist (float): A cut-off value which defines the minimum inter_atomic spacing. + """ + for idx, i in enumerate(structure.sites): + for idx1, j in enumerate(structure.sites): + if struct.distance_matrix[idx][idx1] < dist and idx != idx1: + print(idx, idx1, i, j) + +## POSCAR Modifier ## def set_sd_flags(interface=None, n_layers=2, top=True, bottom=True, lattice_dir=2): """ Set the relaxation flags for top and bottom layers of interface. @@ -231,14 +465,114 @@ def set_sd_flags(interface=None, n_layers=2, top=True, bottom=True, lattice_dir= print('sd_flags', sd_flags) return sd_flags.tolist() -def get_key(my_dict, val): - ''' - Function returns the key corresponding to a dictionary value. +## Workflow Modifier ## +def change_Tasks(original_wf, mode, fw_name_constraint=None, + task_name_constraint='VaspToDb', change_method=None): ''' - for key, value in my_dict.items(): - if val == value: - return key - return "key doesn't exist" + Change the default behavior for task_name_constraint for all fw_name_constraint + by changing, removing, or updating the input of the method. Changes must be + provided using change_method. This setting will completely overwrite the current + method. + + Args: + original_wf (Workflow): Original workflow to change. + mode (str): Mode determines how the VaspToDb method is changed for the + workflow. Options are: replace (with new vasptodb method), update + (the current settings for vasptodb method), and remove (vasptodb + method from workflow). + fw_name_constraint (str): Only apply changes to FWs where fw.name + contains this substring. Default set to None. + task_name_constraint (str): Name of the Firetasks to be tagged, + Default is 'VaspToDb'. + change_method (dict): If update mode, provide a dictionary of + valid arguements for the VaspToDb method in key: value format. + ''' + idx_list = get_fws_and_tasks(original_wf, fw_name_constraint=fw_name_constraint, + task_name_constraint=task_name_constraint) + + if mode == 'change': + for idx_fw, idx_t in idx_list: + original_wf.fws[idx_fw].tasks[idx_t].update(change_method) + + # update all vasptodb methods given the contrainsts provided + if mode == 'update': + for idx_fw, idx_t in idx_list: + original_wf.fws[idx_fw].tasks[idx_t].update(change_method) + + # remove all matching vasptodb methods given the contrainsts provided + if mode == 'remove': + for idx_fw, idx_t in idx_list: + original_wf.fws[idx_fw].tasks.pop(idx_t) + return original_wf + +## Helper Function ## +def slurm_set_npar(): + ''' + Set npar for slurm computing systems. NSLOTS is not valid. + + Returns sqrt(ncores) or npar = 2 + ''' + # set npar looking for ntasks from os.environ: nslots and slurm_cpus_per_task + # fails on all slurm systems + slurm_keys = list(os.environ.keys()) + ncore_keys = [ [key for key in slurm_keys if re.search(match, key) ] + for match in ['TASKS','NSLOTS','CPU','NODE']] + ncore_keys = list(np.unique(sum(ncore_keys, []))) + + slurm_dict = {key:os.environ[key] for key in ncore_keys} + dumpfn(slurm_dict, 'slurm_keys.json') + + for key in ncore_keys: + if re.search('SLURM_NTASKS',key): + ncores = int(os.environ[key]) + break + #elif re.search('SLURM_JOB_CPUS_ON_NODE',key): + # + # ncores = int(os.environ[key]) + # break + elif re.search('NSLOTS',key): + ncores = int(os.environ[key]) + break + npars = [] + npar_found = False + if ncores and ncores != 1: + for npar in range(int(math.sqrt(ncores)),ncores): + if ncores % npar == 0: + npars.append(npar) + npar_found = True + if npar_found: + # the first npar found is the closest to sqrt of ncores + return npars[0] + else: + # returning npar = 2 is a safe default + return 2 + +def decompress(infile, tofile): + ''' + Method to decompress a gzip file. Specify the input file and the output. + Leaves the original gzipped file compressed. + + Args: + infile (str): Absolute or relative file path of file to decompress + using gzip. + tofile (str): Absolute or relative file path of the decompressed file. + ''' + with open(infile, 'rb') as inf, open(tofile, 'w', encoding='utf8') as tof: + decom_str = gzip.decompress(inf.read()).decode('utf-8') + tof.write(decom_str) + +def get_FWjson(): + """ + Helper function which reads the FW.json file in the + local directory and returns the FW.json file. + """ + try: + with gzip.GzipFile('FW.json.gz', 'r') as p: + fw_json_file = json.loads(p.read().decode('utf-8')) + except: + with open('FW.json', 'rb') as p: + fw_json_file = json.load(p) + return fw_json_file def round_decimals_down(number: float, decimals: int = 2): """ @@ -254,7 +588,17 @@ def round_decimals_down(number: float, decimals: int = 2): factor = 10 ** decimals return math.floor(number * factor) / factor -# File Conversions Tools +def get_key(my_dict, val): + ''' + Function returns the key corresponding to a dictionary value. + ''' + for key, value in my_dict.items(): + if val == value: + return key + return "key doesn't exist" + +############################ +### Conversion Functions ### def slab_from_struct(structure, hkl=None): """ Reads a pymatgen structure object and returns a Slab object. Useful for reading @@ -318,61 +662,6 @@ def struct_from_str(string): return Structure(lattice, specie, coords) - -def tag_iface(structure, nlayers_2d, nlayers_sub=2): - """ - Find the atom indices in a heterostructure by specifying how many - layers there are for the 2D material. Returns a dictionary of atom - ids for each layer of the 2D material and nlayers_sub of the - substrate surface. - - Args: - structure (Structure): The hetero_interface structure. - Should be an unrelaxed structure. - nlayers_2d (int): The number of layers contained - within the 2d material. - nalyers_sub (int): The number of layers of the substrate - surface to include in the interface tags. - - Returns: - layer_indices - """ - # set structure as slab - interface = slab_from_struct(structure) - # get z coords of the sites - z_coords = interface.frac_coords[:, 2] - # round z coords to avoid precision errors from np - z_coords_round = [round_decimals_down(i, decimals=3) for i in z_coords] - # get the z-cords of atoms at interface - unique_layer_coords = np.unique(z_coords_round) - - # get all unique coords for each layer - layer_coords_dict = {} - # get 2d layer indices - for idx in list(range(1, nlayers_2d+1)): - layer_coords = unique_layer_coords[-idx] - layer_name = '2d_layer_' + str(abs(idx)) - layer_coords_dict[layer_name] = layer_coords - # get sub layer indices - for count, idx in enumerate(range(nlayers_2d+1, nlayers_2d+nlayers_sub+1)): - layer_coords = unique_layer_coords[-idx] - layer_name = 'sub_layer_' + str(count+1) - layer_coords_dict[layer_name] = layer_coords - - # get atom ids for each layer - layer_indices = {key: [] for key in layer_coords_dict.keys()} - for idx, site in enumerate(structure.sites): - site_coords = round_decimals_down(site.frac_coords[2], decimals=3) - if site_coords in layer_coords_dict.values(): - key = get_key(layer_coords_dict, site_coords) - layer_indices[key].append(idx) - - layer_indices['num_2d_layer'] = nlayers_2d - layer_indices['num_sub_layer'] = nlayers_sub - - return layer_indices - - def show_struct_ase(structure): """ Creates a pop up structure model for a pymatgen structure object using ase's @@ -387,9 +676,9 @@ def show_struct_ase(structure): viewer = view(structure) return viewer - -# Structure Analysis Tools -# TODO: it would be nice to return bond distance statistics to the user +######################################### +### local structure analysis function ### +# TODO: it would be nice to return bond distance statistics to the user @tboland1 class nn_site_indices(CrystalNN): """ This function returns the nearest neighbor atom id's for a @@ -622,53 +911,3 @@ def site_elements(self): nn_elms[t_site] = nn_site_elm return nn_elms - -def iface_layer_locator(structure, cutoff, iface_elements): - """ - Helper function used to locate the iface layers from - the layer solver. - - Args: - structure (Structure): input structure to find the - interface layers for. - cutoff (float): the layer separation cutoff distance. - iface_elements (list): A list of element species - that compose the top and bottom iface layers. - - Returns: - LayerSolver object, top_layer specie, bottom_layer specie - """ - sample = LayerSolver(structure, cutoff=cutoff) - flag = True - - # get the 2d sub separation dist (bot & top atoms) - for i in range(0, len(sample) - 2): - # define layer data - l1 = sample['Layer' + str(i + 1)] - l2 = sample['Layer' + str(i)] - - # sort layers by species - test_layer = sorted([l1['species'][0], - l2['species'][0]]) - - # make sure test_layer matches the iface composition - if test_layer == sorted(iface_elements): - return sample, l1, l2 - - if flag: - print('\t\t## Iface composition not found. Maybe the surface has undergone significant relaxation.') - return None, None, None - -def atomic_distance(structure, dist=2): - """ - Given a structure and an interatomic separation distance all sites which have a - spacing less than dist will be printed out. - - Args: - structure (Structure): The structure to analyze. - dist (float): A cut-off value which defines the minimum inter_atomic spacing. - """ - for idx, i in enumerate(structure.sites): - for idx1, j in enumerate(structure.sites): - if struct.distance_matrix[idx][idx1] < dist and idx != idx1: - print(idx, idx1, i, j) diff --git a/hetero2d/workflow/core.py b/hetero2d/workflow/core.py index 09b6348..4645173 100644 --- a/hetero2d/workflow/core.py +++ b/hetero2d/workflow/core.py @@ -3,28 +3,31 @@ # Distributed under the terms of the GNU License. """ -This module sets up the vasp input files for the workflow to simulate +This module sets up the vasp input files for the workflow to simulate 2D materials adsorbed on a substrate slab. """ from __future__ import division, print_function, unicode_literals, absolute_import - -from uuid import uuid4 from copy import deepcopy +import re +from pymatgen import Structure from pymatgen.analysis.structure_analyzer import SpacegroupAnalyzer -from atomate.vasp.config import VASP_CMD, DB_FILE -from atomate.vasp.powerups import add_namefile -from atomate.utils.utils import get_logger - from fireworks import Firework from fireworks.core.firework import Workflow from fireworks.user_objects.dupefinders.dupefinder_exact import DupeFinderExact +from atomate.utils.utils import get_logger +from atomate.vasp.config import VASP_CMD, DB_FILE +from atomate.vasp.powerups import add_namefile +from atomate.vasp.fireworks.core import StaticFW + +from hetero2d.manipulate.utils import center_slab, change_Tasks from hetero2d.firetasks.heteroiface_tasks import _update_spec -from hetero2d.manipulate.utils import center_slab -from hetero2d.io.VaspInterfaceSet import CMDLInterfaceSet -from hetero2d.fireworks.core import HeteroOptimizeFW, SubstrateSlabFW, GenHeteroStructuresFW +from hetero2d.firetasks.parse_outputs import HeteroAnalysisToDb +from hetero2d.io import CMDLInterfaceSet, CMDLElectronicSet +from hetero2d.fireworks.core import HeteroOptimizeFW, SubstrateSlabFW, GenHeteroStructuresFW, \ + ElectronicFW __author__ = 'Tara M. Boland' @@ -34,265 +37,257 @@ logger = get_logger(__name__) -def get_heterostructures_stabilityWF(struct_2d, struct_sub, struct_3d2d, - heterotransformation_params, slab_params, - user_additions, tags={}, bin_2d=VASP_CMD, - bin_3d=VASP_CMD, dipole=None, uis=None, - uis_2d=None, vis_2d=None, uis_3d2d=None, - vis_3d2d=None, uis_bulk=None, vis_bulk=None, - uis_trans=None, vis_trans=None, uis_iface=None, - vis_iface=None): +def get_heterostructures_stabilityWF(struct_2d, struct_sub, struct_3d2d, heterotransformation_params, + slab_params, user_additions, tags={}, bin_2d=VASP_CMD, bin_3d=VASP_CMD, dipole=None, uis=None, + uis_2d=None, vis_2d=None, uis_3d2d=None, vis_3d2d=None, uis_bulk=None, vis_bulk=None, + uis_trans=None, vis_trans=None, uis_iface=None, vis_iface=None): """ Relax reference structures to determine energetic parameters related - to the formation of 2D films adsorbed onto substrates. Heterointerfaces - are created from the relaxed 2D and substrate slabs which generate a - symmetry matched, low lattice mismatch between the 2D material and - substrate slab. The user should be familiar with the behaviour of - heterointerface function before runnings simulations (load - hetero2d.manipulate.heterotransmuter import hetero_interfaces). Default - relaxation set CMDLInterfaceSet. Specify a complete vasp input + to the formation of 2D films adsorbed onto substrates. Heterointerfaces + are created from the relaxed 2D and substrate slabs which generate a + symmetry matched, low lattice mismatch between the 2D material and + substrate slab. The user should be familiar with the behaviour of + heterointerface function before runnings simulations (load + hetero2d.manipulate.heterotransmuter import hetero_interfaces). Default + relaxation set CMDLInterfaceSet. Specify a complete vasp input set or update the vis set using the uis_i tag to change vasp behaviour. It is also possible to use the general uis to update all vis. - + Args: struct_2d (Structure): The primitive unit cell of the 2D structure to adsorb onto the substrate slab. struct_sub (Structure): The substrate (bulk phase or substrate slab) to adsorb the 2D material. - struct_3d2d (Structure): The bulk reference phase for the 2D + struct_3d2d (Structure): The bulk reference phase for the 2D material. - heterotransformation_params (list): A list of dictionaries where + heterotransformation_params (list): A list of dictionaries where the keys represent the arguments of the hetero_interfaces function - and the value is the argument excluding the structures. Dictionary - keys [{hetero_interfaces function args: value}]. See hetero_interfaces - module for args. - slab_params (dict): Same parameter format as the TransmuterFW to - create a substrate slab. - user_additions (dict): A specification to control the workflow. - See firetasks.heteroiface_tasks._update_spec for a detailed list - of parameters. - bin_2d (str): VASP run command for the VASP version compiled to restrict + and the value is the argument excluding the structures. Dictionary + keys [{hetero_interfaces function args: value}]. See hetero_interfaces + module for args. + slab_params (dict): Same parameter format as the TransmuterFW to + create a substrate slab. + user_additions (dict): A specification to control the workflow. + See firetasks.heteroiface_tasks._update_spec for a detailed list + of parameters. + bin_2d (str): VASP run command for the VASP version compiled to restrict vacuum spacing in z direction from shrinking artifically. - bin_3d (str): VASP run command for the VASP version compiled normally. - + bin_3d (str): VASP run command for the VASP version compiled normally. + Other Parameters: - dipole (bool): If True, dipole corrections will be used for all slab + dipole (bool): If True, dipole corrections will be used for all slab calculations. Defaults to True. - tags (dict): A dictionary applying tags to - general (all), 2d, + tags (dict): A dictionary applying tags to - general (all), 2d, 3d2d, bulk, transmuter, and iface fireworks. - uis (dict): A dictionary of general user incar settings you wish to + uis (dict): A dictionary of general user incar settings you wish to apply to every simulation. Defaults to None. - uis_I (dict): A dictionary of INCAR settings to override the - VASPInputSet defaults or set additional commands for the Ith + uis_I (dict): A dictionary of INCAR settings to override the + VASPInputSet defaults or set additional commands for the Ith structure where I = 2d, 3d2d, bulk, trans, and iface. - vis_I (VASPInputSet): A VASP input set to relax the Ith materials - where I = 2d, 3d2d, bulk, trans, and iface. Defaults to + vis_I (VASPInputSet): A VASP input set to relax the Ith materials + where I = 2d, 3d2d, bulk, trans, and iface. Defaults to CMDLInterfaceSet. Returns: Heterostructure workflow - """ ### Vasp INPUT SET ### # set up default spec user_additions = user_additions or {} unique_id = user_additions.get('unique_id') sub_nodes = user_additions.pop('sub_nodes', 1) - user_additions.update({'struct_sub': struct_sub, - 'struct_2d': struct_2d, 'nlayers_sub': + user_additions.update({'struct_sub': struct_sub, + 'struct_2d': struct_2d, 'nlayers_sub': heterotransformation_params[0]['nlayers_sub']}) - + # spec workflow controls spec = _update_spec(additional_spec=user_additions) # set up the vdw corrections vdw = user_additions.pop('vdw', 'optB88') # optimization flags is_bulk_optimized = spec.get('is_bulk_optimized') - is_2d_optimized = spec.get('is_2d_optimized') + is_2d_optimized = spec.get('is_2d_optimized') is_3d2d_optimized = spec.get('is_3d2d_optimized') is_sub_optimized = spec.get('is_sub_optimized') - + # Vasp Input and wf input from spec wf_name = spec["wf_name"] - if dipole == None: - dipole = True + dipole = True if dipole == None else dipole bin_2d = bin_2d or VASP_CMD bin_3d = bin_3d or VASP_CMD # 2d symmetry data sga = SpacegroupAnalyzer(struct_2d) sg_film = {} - [sg_film.update({key:sga.get_symmetry_dataset()[key]}) + [sg_film.update({key:sga.get_symmetry_dataset()[key]}) for key in ['number','hall_number','international', 'hall','pointgroup']] # bulk symmetry data sga = SpacegroupAnalyzer(struct_sub) sg_sub = {} - [sg_sub.update({key:sga.get_symmetry_dataset()[key]}) + [sg_sub.update({key:sga.get_symmetry_dataset()[key]}) for key in ['number','hall_number','international', 'hall','pointgroup']] # links dict input_fws = [] links_dict = {} - + ##################### ### OPTIMZIE BULK ### if not is_bulk_optimized: - name = "Bulk Structure Optimization: {}".format( - unique_id) + name = "Bulk Structure Optimization: {}".format(unique_id) # symmetry information sa_bulk = SpacegroupAnalyzer(struct_sub) #slab compt struct - struct_sub = sa_bulk.get_conventional_standard_structure() - + struct_sub = sa_bulk.get_conventional_standard_structure() + # VASP calculation controls uis_bulk = uis if not uis_bulk else uis_bulk - vis_bulk = vis_bulk or CMDLInterfaceSet(struct_sub, - vdw = vdw, iface = False, auto_dipole = False, + vis_bulk = vis_bulk or CMDLInterfaceSet(struct_sub, + vdw = vdw, iface = False, auto_dipole = False, user_incar_settings = uis_bulk) - + # database tags - tags_bulk = {'composition': struct_sub.composition.reduced_formula, - 'task_label': name, 'wf_name': wf_name, + tags_bulk = {'composition': struct_sub.composition.reduced_formula, + 'task_label': name, 'wf_name': wf_name, 'spacegroup': sg_sub, 'unique_id': unique_id} [tags_bulk.update(tags.get(i,{})) for i in ['bulk','general'] ] - - # make bulk spec + + # make bulk spec spec_bulk = deepcopy(spec) spec_bulk['tags'] = tags_bulk # optimize bulk structure firework - input_fws.append(HeteroOptimizeFW(spec=spec_bulk, + input_fws.append(HeteroOptimizeFW(spec=spec_bulk, structure=struct_sub, name=name, vasp_input_set=vis_bulk, vasp_cmd=bin_3d, - db_file=DB_FILE)) + db_file=DB_FILE)) # update links dict links_dict.update({'bulk':{'fws':input_fws[0]}}) - ##################### - + ##################### + ############################### - ### OPTIMIZE SUBSTRATE SLAB ### + ### OPTIMIZE SUBSTRATE SLAB ### if not is_sub_optimized: - name = "Slab Structure Optimization-{}: {}".format( - spec['orient'], unique_id) - + name = "Slab Structure Optimization-{}: {}".format(spec['orient'], + unique_id) + # VASP calculation controls uis_trans = uis if not uis_trans else uis_trans - vis_trans = vis_trans or CMDLInterfaceSet(struct_sub, - vdw = vdw, iface = True, auto_dipole = dipole, + vis_trans = vis_trans or CMDLInterfaceSet(struct_sub, + vdw = vdw, iface = True, auto_dipole = dipole, user_incar_settings = uis_trans) - + # database tags tags_trans = {'substrate_composition': struct_sub.composition.reduced_formula, - 'surface_plane': spec.get('orient'), 'unique_id': unique_id, - 'wf_name': wf_name, 'spacegroup': sg_sub} - [tags_trans.update(tags.get(i,{})) + 'surface_plane': spec.get('orient'), 'unique_id': unique_id, + 'wf_name': wf_name, 'spacegroup': sg_sub} + [tags_trans.update(tags.get(i,{})) for i in ['trans','general']] - - # make substrate slab spec + + # make substrate slab spec spec_trans = deepcopy(spec) spec_trans['_queueadapter']={'nodes': sub_nodes} spec_trans['tags'] = tags_trans # optimize substrate slab structure firework copy_vasp_outputs = False if is_bulk_optimized else True - input_fws.append(SubstrateSlabFW(spec=spec_trans, name=name, - structure=struct_sub, slab_params=slab_params, + input_fws.append(SubstrateSlabFW(spec=spec_trans, name=name, + structure=struct_sub, slab_params=slab_params, vasp_input_set=vis_trans, vasp_cmd=bin_3d, db_file=DB_FILE, copy_vasp_outputs=copy_vasp_outputs)) - + # update links dict if not is_bulk_optimized: links_dict['bulk']['links']={input_fws[0]: [input_fws[len(input_fws)-1]]} links_dict.update({'trans':{ 'fws': input_fws[len(input_fws)-1]}}) ############################### - + ###################### ### OPTIMIZE 2D FW ### if not is_2d_optimized: name = "2D Structure Optimization: {}".format(unique_id) - - # center 2d material - struct_2d = center_slab(struct_2d) - + + # center 2d material + struct_2d = center_slab(struct_2d) + # VASP calculation controls uis_2d = uis if not uis_2d else uis_2d - vis_2d = vis_2d or CMDLInterfaceSet(struct_2d, + vis_2d = vis_2d or CMDLInterfaceSet(struct_2d, auto_dipole = dipole, iface=False, vdw = vdw, user_incar_settings = uis_2d) - + # database tags tags_2d = {'composition': struct_2d.composition.reduced_formula, - 'task_label': name, 'wf_name': wf_name, 'spacegroup': sg_film, - 'unique_id': unique_id} + 'task_label': name, 'wf_name': wf_name, 'spacegroup': sg_film, + 'unique_id': unique_id} [tags_2d.update(tags.get(i,{})) for i in ['2d','general']] - + # make 2D spec spec_2d = deepcopy(spec) spec_2d['tags'] = tags_2d - + # optimize 2D structure firework - input_fws.append(HeteroOptimizeFW(spec=spec_2d, structure=struct_2d, + input_fws.append(HeteroOptimizeFW(spec=spec_2d, structure=struct_2d, name=name, vasp_input_set=vis_2d, vasp_cmd=bin_2d, db_file=DB_FILE)) - + # update links dict - links_dict.update({'2d':{'fws':input_fws[len(input_fws)-1]}}) + links_dict.update({'2d':{'fws':input_fws[len(input_fws)-1]}}) ###################### - + ######################## ### OPTIMIZE 3D2D FW ### if not is_3d2d_optimized: name = "3D2D Structure Optimization: {}".format(unique_id) - + # VASP calculation controls uis_3d2d = uis if not uis_3d2d else uis_3d2d - vis_3d2d = vis_3d2d or CMDLInterfaceSet(struct_3d2d, + vis_3d2d = vis_3d2d or CMDLInterfaceSet(struct_3d2d, vdw = vdw, iface = False, user_incar_settings = uis_3d2d) # symmetry data sga = SpacegroupAnalyzer(struct_3d2d) sg_info = {} - [sg_info.update({key:sga.get_symmetry_dataset()[key]}) + [sg_info.update({key:sga.get_symmetry_dataset()[key]}) for key in ['number','hall_number','international','hall','pointgroup']] - + # database tags tags_3d2d={'composition': struct_3d2d.composition.reduced_formula, 'spacegroup': sg_info, - 'task_label': name, 'wf_name': wf_name, 'unique_id': unique_id} + 'task_label': name, 'wf_name': wf_name, 'unique_id': unique_id} [tags_3d2d.update(tags.get(i,[])) for i in ['3d2d','general']] - + # make 3D2D spec spec_3d2d = deepcopy(spec) spec_3d2d['tags'] = tags_3d2d # optimize 3D2D structure firework input_fws.append(HeteroOptimizeFW(spec=spec_3d2d, - structure=struct_3d2d, name=name, + structure=struct_3d2d, name=name, vasp_input_set=vis_3d2d, vasp_cmd=bin_3d, - db_file=DB_FILE)) + db_file=DB_FILE)) # update links dict - links_dict.update({'3d2d':{ - 'fws': input_fws[len(input_fws)-1]}}) + links_dict.update({'3d2d':{'fws': input_fws[len(input_fws)-1]}}) if not is_2d_optimized: links_dict['2d']['links']=input_fws[len(input_fws)-1] elif is_2d_optimized: - links_dict.update({'3d2d':{'fws':input_fws[len(input_fws)-1]}}) + links_dict.update({'3d2d':{'fws':input_fws[len(input_fws)-1]}}) #################################################################################### - + #################################################################################### ### OPTIMIZE Heterostructures WF ### # VASP calculation controls - h_params = {'transformation_params': + h_params = {'transformation_params': heterotransformation_params, 'transformations': ['hetero_interfaces']*len(heterotransformation_params)} uis_iface = uis if not uis_iface else uis_iface - vis_iface = vis_iface or CMDLInterfaceSet(struct_2d, - vdw = vdw, iface = True, auto_dipole = dipole, + vis_iface = vis_iface or CMDLInterfaceSet(struct_2d, + vdw = vdw, iface = True, auto_dipole = dipole, user_incar_settings = uis_iface) # database tags @@ -300,9 +295,8 @@ def get_heterostructures_stabilityWF(struct_2d, struct_sub, struct_3d2d, 'substrate_composition': struct_sub.composition.reduced_formula, 'surface_plane': spec.get('orient'), 'wf_name': wf_name, 'film_spacegroup': sg_film, 'substrate_spacegroup': sg_sub, 'unique_id': unique_id} - [tags_iface.update(tags.get(i,{})) - for i in ['iface','general']] - + [tags_iface.update(tags.get(i,{})) for i in ['iface','general']] + # make iface spec spec_iface = deepcopy(spec) spec_iface['tags'] = tags_iface @@ -310,15 +304,14 @@ def get_heterostructures_stabilityWF(struct_2d, struct_sub, struct_3d2d, # Create heterointerface structure firework # FW name is assigned in the FW below - input_fws.append(GenHeteroStructuresFW(spec=spec_iface, - structure=struct_2d, heterotransformation_params= - h_params, vasp_input_set = vis_iface, - vasp_cmd = bin_3d, db_file = DB_FILE)) + input_fws.append(GenHeteroStructuresFW(spec=spec_iface, structure=struct_2d, + heterotransformation_params=h_params, vasp_input_set = vis_iface, + vasp_cmd = bin_3d, db_file = DB_FILE)) # update links dict links_dict['iface']={'fws': input_fws[len(input_fws)-1]} #################################### - + ########################## ## create fw links dict ## fw_links = {} @@ -326,10 +319,10 @@ def get_heterostructures_stabilityWF(struct_2d, struct_sub, struct_3d2d, if not spec['is_bulk_optimized']: fw_links.update(links_dict['bulk']['links']) # trans, 2d, 3d2d children - child_iface={'is_sub_optimized': 'trans', + child_iface={'is_sub_optimized': 'trans', 'is_2d_optimized': '2d', 'is_3d2d_optimized': '3d2d'} [fw_links.update( - {links_dict[ref]['fws']:[links_dict['iface']['fws']]}) + {links_dict[ref]['fws']:[links_dict['iface']['fws']]}) for is_opt,ref in child_iface.items() if spec[is_opt]==False] # iface children fw_links[links_dict['iface']['fws']]=[] @@ -337,17 +330,239 @@ def get_heterostructures_stabilityWF(struct_2d, struct_sub, struct_3d2d, if not spec['is_3d2d_optimized'] and not spec['is_2d_optimized']: fw_links[links_dict['2d']['fws']].append(links_dict['2d']['links']) ########################## - + ################## FINAL WF NAME ############### name = '{}-on-{}: hkl-[{}]: {}'.format( struct_2d.composition.reduced_formula, struct_sub.composition.reduced_formula, - spec['orient'],spec['unique_id']) + spec['orient'],spec['unique_id']) wf = Workflow(fireworks=input_fws, links_dict=fw_links, name=name) wf = add_namefile(wf) - print('Workflow Name:',name) + print('Workflow Name:', name) print(wf.links) return wf + +def wf_electronic(structure, tags={}, user_additions={}, prev_calc_dir=None, + dos=True, bader=True, cdd=False, user_incar_settings=None, + vasp_input_set=None, vasp_cmd=VASP_CMD, **kwargs): + ''' + Workflow to obtain 1) site-orbital projected DOS, 2) Bader analysis, and 3) + charge density difference for the given structure. If prev_calc_dir is + specified the StaticFW is skipped and a NonSCFFW is performed with reciprocal + k-point density = 100 and NGiF grids 2x the previous values to converge the + charge density for Bader analysis. If cdd is selected the structure spawns 3 + separate fireworks dividing up the input structure. + + Args: + structure (Structure): Input structure. Used to name the workflow if + prev_calc_dir is specified. + tags (dict): A dictionary of tags for the workflow. Can add + individual tags for each structure when doing charge density difference + calculations using 'iso_1' and 'iso_2' keys. + user_additions (dict): A dictionary specifying addtional information and + settings for the workflow. Must provide the 'unique_id'. To override + vasp_input_set defaults for ElectronicFW use 'electronic_set_overrides' + key. Any valid keys for MPNonSCFSet.from_prev_calc() are valid. + A full list keys are: {'grid_density': '0.03 A spacing between points + in NGiF grid', 'unique_id': 1, 'split_idx':{'iso_1':[], 'iso_2':[]}, + 'dedos': 0.05, 'electronic_set_overrides': {}}. + + Other parameters: + prev_calc_dir (str): A path specifying the previous calculation directory. + Retrieves previous calculation output to perform DOS, bader analysis, + and write input files. If None, will create new static calculation + using the provided structure. + dos (bool): If True, peforms high quality site-orbital projected density of + states. Default set to True. + bader (bool): If True, peforms high quality bader analysis for the given + structure. Set 'grid_density' in 'user_additions' to increase the NGiF + grids. Defaults to True with 2x NGiF grid default set by VASP. + cdd (bool): Compute the charge density difference for the given structure. + If True, supply the 'split_idx' key in 'user_additions'. 'split_idx' + can be the dictionary output from hetero2d.manipulate.util.tag_iface or + 2 lists of atom indices under the keys 'iso_1' and 'iso_2'. The indices + produce 2 isolated structures whose charge density is substracted from + the combined structure. Default set to False. + user_incar_settings (dict): INCAR parameters to override StaticFW defaults. + vasp_input_set (VaspInputSet): VASP input set, used to write the input set + for the VASP calculation. Defaults to CMDLInterfaceSet for StaticFW and + CMDLElectronicSet if prev_calc_dir. + vasp_cmd (str): Command to run vasp. + **kwargs (keyword arguments): Other kwargs that are passed to + Firework.__init__ applied to ElectronicFW. + ''' + if not prev_calc_dir and not structure: # ensure struct/prev_calc present + raise ValueError("Must specify structure or previous calculation.") + + wf_name = "{}-Electronic Properties: {}".format( + structure.composition.reduced_formula, + user_additions['unique_id']) # workflow name + + # STATIC INCAR + uis = {"LCHARG": True, "IBRION": -1, "ISMEAR": 0, "SIGMA": 0.05, "NSW": 0} + uis_static = uis.update(user_incar_settings) if user_incar_settings else uis + vis_static = vasp_input_set or CMDLInterfaceSet(structure, + user_incar_settings=uis_static) + electronic_set_overrides = user_additions.get('electronic_set_overrides', {}) + + # create tags + tags.update({'wf_name': wf_name}) # update tags with wf name + tags_iso1, tags_iso2 = [tags.pop(sub_tags) if sub_tags in tags.keys() else {} + for sub_tags in ['iso_1', 'iso_2'] ] + + fws = [] + # COMBINED: STATIC CALCULATION: no prev data + if not prev_calc_dir: # no prev dir; ToDb is removed + static = StaticFW(structure=structure, + name='Static: {}'.format(user_additions['unique_id']), + vasp_input_set=vis_static, + vasp_cmd=vasp_cmd, + prev_calc_loc=None, + prev_calc_dir=prev_calc_dir, + db_file=DB_FILE) + scf = static + fws.append(static) + else: + scf = None + + # NONSCF CALCULATION: from previous dir + # a basic DOS and Bader calc + if not cdd: + electronic = ElectronicFW(name='NSCF: {}'.format(user_additions['unique_id']), + structure=structure, + dedos=user_additions.get('dedos', 0.05), + grid_density=user_additions.get('grid_density', 0.03), + tags=tags, + dos=dos, + bader=bader, + cdd=cdd, + parents=scf, + prev_calc_dir=prev_calc_dir, + vasp_cmd=vasp_cmd, + db_file=DB_FILE, + electronic_set_overrides=electronic_set_overrides, + **kwargs) + fws.append(electronic) + # Charge Density Difference + elif cdd: + split_idx = user_additions.get('split_idx', None) + if split_idx == None: # raise error if cdd dict not present + raise ValueError("You have specified a charge density difference" \ + "calculation wihtout providing a dictionary to split the structure." \ + "Please specify indices using 'split_idx':{} in user_additions.") + + # get iso_1 (2d), iso_2 (substrate) atom indices + if len(split_idx) == 2: + idx_iso1 = split_idx.get('iso_1') + idx_iso2 = split_idx.get('iso_2') + else: + # get a list of atom ids for the 2d material and substrate + idx_iso1 = [el_v for k,v in split_idx.items() if re.search('2d_layer_[\d]', k) + for el_v in v ] + idx_iso2 = list(set(range(0, structure.num_sites)).difference(set(idx_iso1))) + + # generate iso_1 (2d) and iso_2 (substrate) from combined system + struct_iso1 = Structure.from_sites(sites=[structure[i] for i in idx_iso1], + to_unit_cell=True) + struct_iso2 = Structure.from_sites(sites=[structure[i] for i in idx_iso2], + to_unit_cell=True) + + # NONSCF CALCULATION: COMBINED + cdd_combined = ElectronicFW(name='Combined NSCF: {}'.format(user_additions['unique_id']), + structure=structure, + dedos=user_additions.get('dedos', 0.05), + grid_density=user_additions.get('grid_density', 0.03), + tags=tags, + dos=dos, + bader=bader, + cdd=False, + parents=scf, + prev_calc_dir=prev_calc_dir, + vasp_cmd=vasp_cmd, + db_file=DB_FILE, + electronic_set_overrides=electronic_set_overrides, + **kwargs) + + # STATIC CALC: ISO_1 (2d) + tags_iso1.update(tags) + vis_iso1 = deepcopy(vis_static) + vis_iso1.structure = struct_iso1 # 2d + static_iso1 = StaticFW(structure=struct_iso1, + name='ISO 1 Static: {}'.format(user_additions['unique_id']), + vasp_input_set=vis_iso1, + vasp_cmd=vasp_cmd, + prev_calc_loc=None, + prev_calc_dir=None, + parents=None, + db_file=DB_FILE) + + # STATIC CALC: ISO_2 (substrate) + tags_iso2.update(tags) + vis_iso2 = deepcopy(vis_static) + vis_iso2.structure = struct_iso2 # substrate + static_iso2 = StaticFW(structure=struct_iso2, + name='ISO 2 Static: {}'.format(user_additions['unique_id']), + vasp_input_set=vis_iso2, + vasp_cmd=vasp_cmd, + prev_calc_loc=None, + prev_calc_dir=None, + parents=None, + db_file=DB_FILE) + + # NONSCF CALCULATION: ISO_1 ISO_2 + cdd_iso1 = ElectronicFW(name='ISO 1 NSCF: {}'.format(user_additions['unique_id']), + structure=struct_iso1, + dedos=user_additions.get('dedos', 0.05), + grid_density=user_additions.get('grid_density', 0.03), + tags=tags_iso1, + dos=dos, + bader=bader, + cdd=False, + parents=[static_iso1], + prev_calc_dir=None, + vasp_cmd=vasp_cmd, + db_file=DB_FILE, + electronic_set_overrides=electronic_set_overrides, + **kwargs) + cdd_iso2 = ElectronicFW(name='ISO 2 NSCF: {}'.format(user_additions['unique_id']), + structure=struct_iso2, + dedos=user_additions.get('dedos', 0.05), + grid_density=user_additions.get('grid_density', 0.03), + tags=tags_iso2, + dos=dos, + bader=bader, + cdd=False, + parents=[static_iso2], + prev_calc_dir=None, + vasp_cmd=vasp_cmd, + db_file=DB_FILE, + electronic_set_overrides=electronic_set_overrides, + **kwargs) + + cdd_analysis = Firework( + HeteroAnalysisToDb(db_file=DB_FILE, + task_label="Charge Density Difference Analysis", + dos=False, + bader=False, + cdd=cdd, + additional_fields={}), + name="Charge Density Difference Analysis", + spec={"_allow_fizzled_parents": False}, + parents=[cdd_combined, cdd_iso1, cdd_iso2]) + [fws.append(i) for i in [static_iso1, static_iso2, cdd_iso1, cdd_iso2, + cdd_combined, cdd_analysis]] + + # CREATE WORKFLOW + wf = Workflow(fireworks=fws, name=wf_name) + wf = add_namefile(wf) + + # remove VaspToDb for static calculations + wf = change_Tasks(wf, mode='remove', fw_name_constraint=None, + task_name_constraint='VaspToDb') + + print(wf.name) + print(wf.links) + return wf