EoS thermal workflow in current version of MPWorks

mpworks/drones/mp_vaspdrone.py

@@ -249,6 +249,9 @@ def process_fw(self, dir_name, d):
             if 'deformed' in d['task_type']:
                 d['deformation_matrix'] = fw_dict['spec']['deformation_matrix']
                 d['original_task_id'] = fw_dict['spec']['original_task_id']
+            if 'modified volume' in d['task_type']:
+                d['strainfactor'] = fw_dict['spec']['strainfactor']
+                d['original_task_id'] = fw_dict['spec']['original_task_id']
             if not self.update_duplicates:
                 if 'optimize structure' in d['task_type'] and 'output' in d:
                     # create a new SNL based on optimized structure

mpworks/firetasks/eos_thermal_tasks.py

@@ -0,0 +1,232 @@
+from monty.os.path import zpath
+import os
+import json
+from pymongo import MongoClient
+import numpy as np
+from decimal import Decimal
+
+__author__ = 'Cormac Toher'
+
+
+from fireworks.utilities.fw_serializers import FWSerializable
+from fireworks.core.firework import FireTaskBase, FWAction
+from pymatgen.io.vasp.inputs import Incar, Poscar, Kpoints
+from pymatgen.analysis.eqn_of_state_thermal.eos_vasp_setup import ModifiedVolumeStructureSet
+from pymatgen.analysis.eqn_of_state_thermal.eos_wf_run import eos_thermal_properties
+from fireworks.core.firework import Firework, Workflow
+from mpworks.firetasks.vasp_io_tasks import VaspWriterTask, VaspToDBTask
+from mpworks.firetasks.custodian_task import get_custodian_task
+from fireworks.utilities.fw_utilities import get_slug
+from pymatgen import Composition
+from pymatgen.matproj.snl import StructureNL
+from mpworks.workflows import snl_to_wf
+from mpworks.firetasks.snl_tasks import AddSNLTask
+from mpworks.snl_utils.mpsnl import get_meta_from_structure, MPStructureNL
+from pymatgen.core.structure import Structure
+from mpworks.workflows.wf_settings import QA_VASP, QA_DB, QA_VASP_SMALL
+
+def update_spec_force_convergence(spec, user_vasp_settings=None):
+    fw_spec = spec
+    update_set = {"ENCUT": 700, "EDIFF": 0.000001, "ALGO":"N", "NPAR":2, "NSW":0}
+    if user_vasp_settings and user_vasp_settings.get("incar"):
+            update_set.update(user_vasp_settings["incar"])
+    fw_spec['vasp']['incar'].update(update_set)
+    old_struct=Poscar.from_dict(fw_spec["vasp"]["poscar"]).structure
+    if user_vasp_settings and user_vasp_settings.get("kpoints"):
+        kpoints_density = user_vasp_settings["kpoints"]["kpoints_density"]
+    else:
+        kpoints_density = 7000
+    k=Kpoints.automatic_density(old_struct, kpoints_density)
+    fw_spec['vasp']['kpoints'] = k.as_dict()
+    return fw_spec
+
+
+class SetupFConvergenceTask(FireTaskBase, FWSerializable):
+    _fw_name = "Setup Force Convergence Task"
+
+    def run_task(self, fw_spec):
+        incar = fw_spec['vasp']['incar']
+        update_set = {"ENCUT": 600, "EDIFF": 0.00005}
+        incar.update(update_set)
+        kpoints = fw_spec['vasp']['kpoints']
+        k = [int(round(2.5*k)) if int(round(2.5*k))%2 else int(round(2.5*k))+1 for k in kpoints['kpoints'][0]]
+        kpoints['kpoints'] = [k]
+        return FWAction()
+
+class SetupEoSThermalTask(FireTaskBase, FWSerializable):
+    _fw_name = "Setup EoS Thermal Task"
+
+    def run_task(self, fw_spec):
+        incar = Incar.from_file(zpath("INCAR"))
+        incar.update({"NSW": 0})
+        incar.update({"LCHARG": True})
+        incar.update({"IBRION": -1})
+        incar.write_file("INCAR")
+        return FWAction()
+
+class SetupModifiedVolumeStructTask(FireTaskBase, FWSerializable):
+    _fw_name = "Setup Modified Volume Struct Task"
+
+    def run_task(self, fw_spec):
+        # Read structure from previous relaxation
+        relaxed_struct = fw_spec['output']['crystal']
+        modified_struct_set = ModifiedVolumeStructureSet(relaxed_struct)
+        poisson_val = fw_spec['poisson_ratio']
+        wf=[]
+        task_id_list=[]
+        for i, mod_struct in enumerate(modified_struct_set.modvol_structs):
+            fws=[]
+            connections={}
+            f = Composition(mod_struct.formula).alphabetical_formula
+            snl = StructureNL(mod_struct, 'Cormac Toher <cormac.toher@duke.edu>',projects=["Thermal"])
+            tasks = [AddSNLTask()]
+            snl_priority = fw_spec.get('priority', 1)
+            spec = {'task_type': 'Add Modified Struct to SNL database', 'snl': snl.as_dict(),
+                    '_queueadapter': QA_DB, '_priority': snl_priority}
+            if 'snlgroup_id' in fw_spec and isinstance(snl, MPStructureNL):
+                spec['force_mpsnl'] = snl.as_dict()
+                spec['force_snlgroup_id'] = fw_spec['snlgroup_id']
+                del spec['snl']
+            fws.append(Firework(tasks, spec, name=get_slug(f + '--' + spec['task_type']), fw_id=-1000+i*10))
+            connections[-1000+i*10] = [-999+i*10]
+            spec = snl_to_wf._snl_to_spec(snl, parameters={'exact_structure':True})
+            spec = update_spec_force_convergence(spec)
+            spec['strainfactor'] = modified_struct_set.strainfactors[i]
+            spec['original_task_id']=fw_spec["task_id"]
+            spec['_priority'] = fw_spec['_priority']*2
+            # Turn off dupefinder for modified structure
+            del spec['_dupefinder']
+            spec['task_type'] = "Calculate static modified volume structure"
+            fws.append(Firework([VaspWriterTask(), SetupEoSThermalTask(), get_custodian_task(spec)],
+                                spec, name=get_slug(f + '--' + fw_spec['task_type']), fw_id=-999+i*10))
+
+            priority = fw_spec['_priority']*3
+            spec = {'task_type': 'VASP db insertion', '_priority': priority,
+                    '_allow_fizzled_parents': True, '_queueadapter': QA_DB, 'poisson_ratio': poisson_val,
+                    'eqn_of_state_thermal':"modified_structure", 'clean_task_doc':True,
+                    'strainfactor':modified_struct_set.strainfactors[i], 'original_task_id':fw_spec["task_id"], 'original_task_id_list':task_id_list}
+            fws.append(Firework([VaspToDBTask(), AddEoSThermalDataToDBTask()], spec, name=get_slug(f + '--' + spec['task_type']), fw_id=-998+i*10))
+            connections[-999+i*10] = [-998+i*10]
+            wf.append(Workflow(fws, connections))
+        return FWAction(additions=wf)        
+
+
+class AddEoSThermalDataToDBTask(FireTaskBase, FWSerializable):
+    _fw_name = "Add EoS Thermal Data to DB Task"
+
+    def run_task(self, fw_spec):
+        db_dir = os.environ['DB_LOC']
+        db_path = os.path.join(db_dir, 'tasks_db.json')
+	poisson_val = fw_spec['poisson_ratio']
+	strainfactor_val = fw_spec['strainfactor']
+        task_id_list = fw_spec['original_task_id_list']
+        i = fw_spec['original_task_id']
+        j = fw_spec['task_id']
+
+        with open(db_path) as f:
+            db_creds = json.load(f)
+        connection = MongoClient(db_creds['host'], db_creds['port'])
+        tdb = connection[db_creds['database']]
+        tdb.authenticate(db_creds['admin_user'], db_creds['admin_password'])
+        tasks = tdb[db_creds['collection']]
+        eos_thermal = tdb['eos_thermal']
+        ndocs = tasks.find({"original_task_id": i, 
+                            "state":"successful"}).count()
+        existing_doc = eos_thermal.find_one({"relaxation_task_id" : i})
+        if existing_doc:
+            print "Updating: " + i
+        else:
+            print "New material: " + i
+        d = {"analysis": {}, "error": [], "warning": []}
+        d["ndocs"] = ndocs
+	o = tasks.find_one({"task_id" : i},
+                           {"pretty_formula" : 1, "spacegroup" : 1,
+                            "snl" : 1, "snl_final" : 1, "run_tags" : 1})
+        if not o:
+            raise ValueError("Cannot find original task id")
+	if o: 
+	    print("o exists")
+        d["strain_tasks"] = {}
+        volume_values = []
+        energy_values = []
+        for k in tasks.find({"original_task_id": i}, {"strainfactor":1, "calculations.output":1, "state":1, "task_id":1}):     
+            kenerg_atom = k['calculations'][0]['output']['final_energy_per_atom']
+            kvol = k['calculations'][0]['output']['crystal']['lattice']['volume']
+            ksites = k['calculations'][0]['output']['crystal']['sites']
+            knatoms = len(ksites)
+            kenerg = kenerg_atom * knatoms
+            energy_values.append(kenerg)
+            volume_values.append(kvol)
+        d["snl"] = o["snl"]
+        if "run_tags" in o.keys():
+            d["run_tags"] = o["run_tags"]
+            for tag in o["run_tags"]:
+                if isinstance(tag, dict):
+                    if "input_id" in tag.keys():
+                        d["input_mp_id"] = tag["input_id"]
+        d["snl_final"] = o["snl_final"]
+        d["pretty_formula"] = o["pretty_formula"]
+
+        # Old input mp-id style
+        if o["snl"]["about"].get("_mp_id"):
+            d["material_id"] = o["snl"]["about"]["_mp_id"]
+
+        # New style
+        elif "input_mp_id" in d:
+            d["material_id"] = d["input_mp_id"]
+        else:
+            d["material_id"] = None
+        d["relaxation_task_id"] = i
+
+        calc_struct = Structure.from_dict(o["snl_final"])
+        d["spacegroup"]=o.get("spacegroup", "Unknown")
+
+
+        if ndocs >= 21:
+            eos_thermal_dict = {}
+            # Perform thermal equation of state fitting and analysis
+            eos_thermal_properties_inst = eos_thermal_properties()
+            eos_thermal_dict = eos_thermal_properties_inst.eos_thermal_run(calc_struct, volume_values, energy_values, poissonratio=poisson_val)
+
+            # Add equation of state results to dict
+            d["temperature"] = eos_thermal_dict["temperature"]
+            d["pressure"] = eos_thermal_dict["pressure"]
+            d["Thermal_conductivity_temp_list"] = eos_thermal_dict["thermal_conductivity"]
+            d["Debye_temperature_temp_list"] = eos_thermal_dict["Debye_temperature"]
+            d["Gruneisen_parameter_temp_list"] = eos_thermal_dict["Gruneisen_parameter"]
+            d["Heat_capacity_Cv_temp_list"] = eos_thermal_dict["Heat_capacity_Cv"]
+            d["Heat_capacity_Cp_temp_list"] = eos_thermal_dict["Heat_capacity_Cp"]
+            d["Volume_temp_list"] = eos_thermal_dict["Volume"]
+            d["Bulk_modulus_temp_list"] = eos_thermal_dict["Bulk_modulus"]
+            d["BM_coeffs"] = eos_thermal_dict["BM_coeffs"].tolist()
+            jtdbest = eos_thermal_dict["Best_fit_temperature"]
+
+            # Add values at specific temperatures (Debye temperature and Gruneisen parameter best fit, values at 300K)
+            d["Debye_temperature"] =  d["Debye_temperature_temp_list"][eos_thermal_dict["Best_fit_temperature"]]
+            d["Gruneisen_parameter"] = d["Gruneisen_parameter_temp_list"][eos_thermal_dict["Best_fit_temperature"]]
+            d["Thermal_conductivity_300K"] = d["Thermal_conductivity_temp_list"][eos_thermal_dict["300K_point"]]
+            d["Heat_capacity_Cv_300K"] = d["Heat_capacity_Cv_temp_list"][eos_thermal_dict["300K_point"]]
+            d["Heat_capacity_Cp_300K"] = d["Heat_capacity_Cp_temp_list"][eos_thermal_dict["300K_point"]]
+            d["Bulk_modulus_300K"] = d["Bulk_modulus_temp_list"][eos_thermal_dict["300K_point"]]
+
+
+
+        else:
+            d['state'] = "Fewer than 21 successful tasks completed"
+            print("Fewer than 21 successful tasks completed")
+            print("ndocs = ", ndocs)
+            print("mp_id = ", i)
+	    return FWAction()
+
+        if o["snl"]["about"].get("_kpoint_density"):
+            d["kpoint_density"]= o["snl"]["about"].get("_kpoint_density")
+
+        if d["error"]:
+            raise ValueError("Thermal equation of state analysis failed: {}".format(d["error"]))
+        elif d['analysis'].get('filter_pass') is None or d["analysis"]["filter_pass"]:
+            d["state"] = "successful"
+        else:
+            d["state"] = "filter_failed"
+        eos_thermal.update({"relaxation_task_id": d["relaxation_task_id"]}, 
+                           d, upsert=True)
+        return FWAction()

mpworks/processors/process_submissions.py

@@ -5,6 +5,7 @@
 from mpworks.submission.submission_mongo import SubmissionMongoAdapter
 from mpworks.workflows.snl_to_wf import snl_to_wf
 from mpworks.workflows.snl_to_wf_elastic import snl_to_wf_elastic
+from mpworks.workflows.snl_to_wf_eos_thermal import snl_to_wf_eos_thermal
 from mpworks.workflows.wf_utils import NO_POTCARS
 from pymatgen.matproj.snl import StructureNL
 
@@ -77,7 +78,9 @@ def submit_new_workflow(self):
 
                     # create a workflow
                     if "Elasticity" in snl.projects:
-                        wf=snl_to_wf_elastic(snl, job['parameters'])
+                        wf = snl_to_wf_elastic(snl, job['parameters'])
+                    if "EosThermal" in snl.projects:
+                        wf = snl_to_wf_eos_thermal(snl, job['parameters'])
                     else:
                         wf = snl_to_wf(snl, job['parameters'])
                     self.launchpad.add_wf(wf)

mpworks/processors/submit_canonical.py

@@ -46,37 +46,61 @@ def clear_env():
     db['band_structure_fs.files'].remove()
 
 
-def submit_tests(names=None, params=None):
+def submit_tests(names=None, params=None, project=None):
     sma = SubmissionMongoAdapter.auto_load()
 
     # note: TiO2 is duplicated twice purposely, duplicate check should catch this
-    compounds = {"Si": 149, "Al": 134, "ZnO": 2133, "FeO": 18905,
-                 "LiCoO2": 601860, "LiFePO4": 585433, "GaAs": 2534, "Ge": 32, "PbTe": 19717,
-                 "YbO": 1216, "SiC": 567551, "Fe3C": 510623, "SiO2": 547211, "Na2O": 2352,
-                 "InSb (unstable)": 10148, "Sb2O5": 1705, "N2O5": 554368, "BaTiO3": 5020,
-                 "Rb2O": 1394, "TiO2": 554278, "TiO2 (2)": 554278, 'BaNbTePO8': 560794,
-                 "AgCl": 22922, "AgCl (2)": 570858, "SiO2 (2)": 555211, "Mg2SiO4": 2895, "CO2": 20066,
-                 "PbSO4": 22298, "SrTiO3": 5532, "FeAl": 2658, "AlFeCo2": 10884, "NaCoO2": 554427,
-                 "ReO3": 547271, "LaH2": 24153, "SiH3I": 28538, "LiBH4": 30209, "H8S5N2": 28143,
-                 "LiOH": 23856, "SrO2": 2697, "Mn": 35, "Hg4Pt": 2312,
-                 "PdF4": 13868, "Gd2WO6": 651333, 'MnO2': 19395, 'VO2': 504800}
+    compounds = {
+	"InSb": 20012, "AgCl": 22922, "InSe": 20485, "CdGeAs2": 4953,
+	"CdGeP2": 3668, "BN": 1639, "GaSb": 1156, "RbBr": 22867, "NaCl": 22862,
+	"GaP": 2490, "GaN": 804, "SnTe": 1883, "AlAs": 2172, "BP": 1479,
+	"ZnGeAs2": 4008, "KBr": 23251, "PbSe": 2201, "CdSe": 2691, "CuGaS2": 5238,
+	"Al2O3": 1143, "LiF": 1138, "NaF": 682, "BeO": 2542, "NaI": 23268,
+	"LiH": 23703, "AgGaS2": 5342, "CuInTe2": 22261, "SiC": 1883, "SrO": 2472,
+	"ZnSiAs2": 3595, "CuGaSe2": 4840, "InN": 22205, "NaBr": 22916, "ZnTe": 2176,
+	"MgO": 1265, "AlN": 661, "CuGaTe2": 3839, "PbTe": 19717, "KCl": 23193,
+	"ZnGeP2": 4524, "InP": 20351, "Sb2Te3": 1201, "Bi2Se3": 541837, "Si": 149,
+	"RbCl": 22867, "CdS": 672, "CdTe": 406, "Cr2O3": 19399, "HgTe": 2730,
+	"AlP": 1550, "C_diamond": 66, "GaAs": 2534, "KF": 463, "PbS": 21276,
+	"IrSb3": 1239, "KI": 22898, "HgSe": 820, "ZnO": 2133, "InAs": 20305,
+	"Bi2Te3": 34202, "ZnS": 10695, "Ge": 32, "CaO": 2605, "CoSb3": 1317,
+	"RbI": 22903, "ZnSe": 1190, "BaO": 1342, "AlSb": 2624, "ZnSb": 753
+    } if project == 'EosThermal' else {
+	"Si": 149, "Al": 134, "ZnO": 2133, "FeO": 18905, "LiCoO2": 601860,
+	"LiFePO4": 585433, "GaAs": 2534, "Ge": 32, "PbTe": 19717, "YbO": 1216,
+	"SiC": 567551, "Fe3C": 510623, "SiO2": 547211, "Na2O": 2352, "InSb (unstable)": 10148,
+	"Sb2O5": 1705, "N2O5": 554368, "BaTiO3": 5020, "Rb2O": 1394, "TiO2": 554278,
+	"TiO2 (2)": 554278, 'BaNbTePO8': 560794, "AgCl": 22922, "AgCl (2)": 570858,
+	"SiO2 (2)": 555211, "Mg2SiO4": 2895, "CO2": 20066, "PbSO4": 22298,
+	"SrTiO3": 5532, "FeAl": 2658, "AlFeCo2": 10884, "NaCoO2": 554427,
+	"ReO3": 547271, "LaH2": 24153, "SiH3I": 28538, "LiBH4": 30209,
+	"H8S5N2": 28143, "LiOH": 23856, "SrO2": 2697, "Mn": 35, "Hg4Pt": 2312,
+	"PdF4": 13868, "Gd2WO6": 651333, 'MnO2': 19395, 'VO2': 504800
+    }
 
     mpr = MPRester()
 
     for name, sid in compounds.iteritems():
         if not names or name in names:
-            sid = mpr.get_materials_id_from_task_id("mp-{}".format(sid))
+	    mpid = "mp-{}".format(sid)
+            sid = mpr.get_materials_id_from_task_id(mpid)
             s = mpr.get_structure_by_material_id(sid, final=False)
+	    query_key = 'elasticity.poisson_ratio'
+	    poisson_ratio = mpr.query(mpid, [query_key])[0][query_key]
 
-            snl = StructureNL(s, 'Anubhav Jain <anubhavster@gmail.com>')
+            snl = StructureNL(s, 'Patrick Huck <phuck@lbl.gov>')
+	    if project is not None:
+		snl.projects.append(project)
 
             parameters = {'priority': 10} if name == 'Si' else {}
+	    if poisson_ratio is not None:
+	        parameters['poisson_ratio'] = poisson_ratio
             if params:
                 parameters.update(params)
-            sma.submit_snl(snl, 'anubhavster@gmail.com', parameters=parameters)
+            sma.submit_snl(snl, 'phuck@lbl.gov', parameters=parameters)
 
 
-def clear_and_submit(clear=False, names=None, params=None):
+def clear_and_submit(clear=False, names=None, params=None, project=None):
     if clear:
         clear_env()
-    submit_tests(names=names, params=params)
+    submit_tests(names=names, params=params, project=project)

mpworks/scripts/submit_canonical_run.py

@@ -27,6 +27,7 @@ def go_testing():
     parser = ArgumentParser(description=m_description)
     parser.add_argument('-c', '--clear', help='clear old databases', action='store_true')
     parser.add_argument('-n', '--names', help='csv of compound names', default=None)
+    parser.add_argument('-p', '--project', help='project name to trigger specific worflow', default=None)
     parser.add_argument('--noboltztrap', help='do NOT run boltztrap', action='store_true')
     parser.add_argument('--exact', help='exact structure', action='store_true')
     args = parser.parse_args()
@@ -38,7 +39,7 @@ def go_testing():
         params['boltztrap'] = False
     if args.exact:
         params['exact_structure'] = True
-    clear_and_submit(args.clear, names, params)
+    clear_and_submit(args.clear, names, params, args.project)
 
 if __name__ == '__main__':
     go_testing()