Pynta calculate IRC in different mode by user's choice

pynta/main.py

@@ -22,6 +22,11 @@
 import fireworks.fw_config
 import logging
 
+#logger
+logger = logging.getLogger(__name__)
+logging.basicConfig(filename='pynta.log', level=logging.INFO)
+
+
 class Pynta:
     def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fworker_path=None,
         vacuum=8.0,repeats=(3,3,4),slab_path=None,software="Espresso", pbc=(True,True,False),socket=False,queue=False,njobs_queue=0,a=None,
@@ -33,6 +38,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw
                             }, },
         software_kwargs_gas=None,
         TS_opt_software_kwargs=None,
+        irc_mode=None, #choose irc mode: 'skip', 'relaxed', 'fixed'
         lattice_opt_software_kwargs={'kpts': (25,25,25), 'ecutwfc': 70, 'degauss':0.02, 'mixing_mode': 'plain'},
         reset_launchpad=False,queue_adapter_path=None,num_jobs=25,max_num_hfsp_opts=None,#max_num_hfsp_opts is mostly for fast testing
         Eharmtol=3.0,Eharmfiltertol=30.0,Ntsmin=5,frozen_layers=2,fmaxopt=0.05,fmaxirc=0.1,fmaxopthard=0.05):
@@ -59,6 +65,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw
         self.metal = metal
         self.adsorbate_fw_dict = dict()
         self.software_kwargs = software_kwargs
+        self.irc_mode = irc_mode
 
         if software.lower() == 'vasp':
             self.pbc = (True,True,True)
@@ -116,6 +123,8 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw
         self.fmaxirc = fmaxirc
         self.fmaxopthard = fmaxopthard
 
+        logger.info('Pynta class is initiated')
+
     def generate_slab(self,skip_launch=False):
         """
         generates and optimizes a small scale slab that can be scaled to a large slab as needed
@@ -129,6 +138,7 @@ def generate_slab(self,skip_launch=False):
             self.a = a
         else:
             a = self.a
+        logger.info('Construct slab with optimal lattice constant')
         #construct slab with optimial lattice constant
         slab = slab_type(self.metal,self.repeats,a,self.vacuum)
         slab.pbc = self.pbc
@@ -440,24 +450,52 @@ def setup_transition_states(self,adsorbates_finished=False):
         Sets up fireworks to generate and filter a set of TS estimates, optimize each unique TS estimate,
         and run vibrational and IRC calculations on the each unique final transition state
         Note the vibrational and IRC calculations are launched at the same time
+
+        If irc_mode is "fixed", entire slab layers are fixed. 
+        If irc_mode is "relaxed", bottom half of slab layers are fixed and top half of slab layers are relaxed.
+        If irc_mode is not "fixed" nor "relaxed", IRC is not calculated
         """
         if self.software != "XTB":
             opt_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs_TS,
                 "run_kwargs": {"fmax" : self.fmaxopt, "steps" : 70},"constraints": ["freeze up to {}".format(self.freeze_ind)],"sella":True,"order":1,}
         else:
             opt_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs_TS,
                 "run_kwargs": {"fmax" : 0.02, "steps" : 70},"constraints": ["freeze up to "+str(self.nslab)],"sella":True,"order":1,}
+
         vib_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs,
                 "constraints": ["freeze up to "+str(self.nslab)]}
-        IRC_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs,
-                "run_kwargs": {"fmax" : self.fmaxirc, "steps" : 70},"constraints":["freeze up to "+str(self.nslab)]}
+
+        #logging.info
+        print(f"================= IRC mode is: {self.irc_mode} =======================")
+        logger.info(f"================= IRC mode is: {self.irc_mode} =======================")
+        #pass through 
+
         for i,rxn in enumerate(self.rxns_dict):
+            #if irc_mode is "fixed" freeze all slab and conduct MolecularTSEstimate. 
+            if self.irc_mode == "fixed":
+                print("==Entire slab layers are frozen==")
+                logging.info("Entire slab layers are frozen")
+                IRC_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs,
+                    "run_kwargs": {"fmax" : self.fmaxopt, "steps" : 70},"constraints": ["freeze up to "+str(self.nslab)]}
+
+            elif self.irc_mode == "relaxed":
+                print("==Top half of the slab is relaxed==")
+                logger.info("==Top half of the slab is relaxed==")
+                IRC_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs,
+                    "run_kwargs": {"fmax" : self.fmaxopt, "steps" : 70},"constraints": ["freeze up to {}".format(self.freeze_ind)]}
+        # if irc_mode = "skip" : do not conduct IRC
+            else:
+                print("==Skip IRC: IRC is not conducted==")
+                logger.info("==Skip IRC: IRC is not conducted==")
+                IRC_obj_dict = {}
+                pass
+
             ts_path = os.path.join(self.path,"TS"+str(i))
             os.makedirs(ts_path)
             ts_task = MolecularTSEstimate({"rxn": rxn,"ts_path": ts_path,"slab_path": self.slab_path,"adsorbates_path": os.path.join(self.path,"Adsorbates"),
-                "rxns_file": self.rxns_file,"path": self.path,"metal": self.metal,"facet": self.surface_type, "out_path": ts_path,
-                "spawn_jobs": True, "opt_obj_dict": opt_obj_dict, "vib_obj_dict": vib_obj_dict,
-                    "IRC_obj_dict": IRC_obj_dict, "nprocs": 48, "name_to_adjlist_dict": self.name_to_adjlist_dict,
+                    "rxns_file": self.rxns_file,"path": self.path,"metal": self.metal,"facet": self.surface_type, "out_path": ts_path, "irc_mode": self.irc_mode,
+                    "spawn_jobs": True, "opt_obj_dict": opt_obj_dict, "vib_obj_dict": vib_obj_dict, "IRC_obj_dict": IRC_obj_dict,
+                    "nprocs": 48, "name_to_adjlist_dict": self.name_to_adjlist_dict,
                     "gratom_to_molecule_atom_maps":{sm: {str(k):v for k,v in d.items()} for sm,d in self.gratom_to_molecule_atom_maps.items()},
                     "gratom_to_molecule_surface_atom_maps":{sm: {str(k):v for k,v in d.items()} for sm,d in self.gratom_to_molecule_surface_atom_maps.items()},
                     "nslab":self.nslab,"Eharmtol":self.Eharmtol,"Eharmfiltertol":self.Eharmfiltertol,"Ntsmin":self.Ntsmin,

pynta/tasks.py

@@ -468,15 +468,16 @@ def run_task(self, fw_spec):
 class MolecularTSEstimate(FiretaskBase):
     required_params = ["rxn","ts_path","slab_path","adsorbates_path","rxns_file","path","metal","facet",
                         "name_to_adjlist_dict", "gratom_to_molecule_atom_maps",
-                        "gratom_to_molecule_surface_atom_maps","opt_obj_dict",
-                                "vib_obj_dict","IRC_obj_dict","nslab","Eharmtol","Eharmfiltertol","Ntsmin","max_num_hfsp_opts"]
-    optional_params = ["out_path","spawn_jobs","nprocs",]
+                        "gratom_to_molecule_surface_atom_maps","irc_mode",
+                        "vib_obj_dict","opt_obj_dict","nslab","Eharmtol","Eharmfiltertol","Ntsmin","max_num_hfsp_opts"]
+    optional_params = ["out_path","spawn_jobs","nprocs"]
     def run_task(self, fw_spec):
         gratom_to_molecule_atom_maps = {sm: {int(k):v for k,v in d.items()} for sm,d in self["gratom_to_molecule_atom_maps"].items()}
         gratom_to_molecule_surface_atom_maps = {sm: {int(k):v for k,v in d.items()} for sm,d in self["gratom_to_molecule_surface_atom_maps"].items()}
         out_path = self["out_path"] if "out_path" in self.keys() else ts_path
         spawn_jobs = self["spawn_jobs"] if "spawn_jobs" in self.keys() else False
         nprocs = self["nprocs"] if "nprocs" in self.keys() else 1
+        IRC_obj_dict = self["IRC_obj_dict"] if "IRC_obj_dict" in self.keys() else None
 
         ts_path = self["ts_path"]
         rxn = self["rxn"]
@@ -490,6 +491,7 @@ def run_task(self, fw_spec):
         max_num_hfsp_opts = self["max_num_hfsp_opts"]
         slab_path = self["slab_path"]
         slab = read(slab_path)
+        irc_mode = self["irc_mode"]
 
         cas = SlabAdsorptionSites(slab,facet,allow_6fold=False,composition_effect=False,
                             label_sites=True,
@@ -614,17 +616,28 @@ def run_task(self, fw_spec):
                 xyzsout.append(xyzs[Eind])
 
         if spawn_jobs:
-            irc_obj_dict_forward = deepcopy(self["IRC_obj_dict"])
-            irc_obj_dict_forward["forward"] = True
-            irc_obj_dict_reverse = deepcopy(self["IRC_obj_dict"])
-            irc_obj_dict_reverse["forward"] = False
-
-            ctask = MolecularCollect({"xyzs":xyzsout,"check_symm":True,"fw_generators": ["optimize_firework",["vibrations_firework","IRC_firework","IRC_firework"]],
-                "fw_generator_dicts": [self["opt_obj_dict"],[self["vib_obj_dict"],irc_obj_dict_forward,irc_obj_dict_reverse]],
+            print(irc_mode)
+            if irc_mode == "relaxed" or irc_mode == "fixed":
+                irc_obj_dict_forward = deepcopy(self["IRC_obj_dict"])
+                irc_obj_dict_forward["forward"] = True
+                irc_obj_dict_reverse = deepcopy(self["IRC_obj_dict"])
+                irc_obj_dict_reverse["forward"] = False
+
+                ctask = MolecularCollect({"xyzs":xyzsout,"check_symm":True,"fw_generators": ["optimize_firework",["vibrations_firework","IRC_firework","IRC_firework"]],
+                    "fw_generator_dicts": [self["opt_obj_dict"],[self["vib_obj_dict"],irc_obj_dict_forward,irc_obj_dict_reverse]],
                     "out_names": ["opt.xyz",["vib.json","irc_forward.traj","irc_reverse.traj"]],"future_check_symms": [True,False], "label": "TS"+str(index)+"_"+rxn_name})
-            cfw = Firework([ctask],name="TS"+str(index)+"_"+rxn_name+"_collect",spec={"_allow_fizzled_parents": True, "_priority": 5})
-            newwf = Workflow([cfw],name='rxn_'+str(index)+str(rxn_name))
-            return FWAction(detours=newwf) #using detour allows us to inherit children from the original collect to the subsequent collects
+                cfw = Firework([ctask],name="TS"+str(index)+"_"+rxn_name+"_collect",spec={"_allow_fizzled_parents": True, "_priority": 5})
+                newwf = Workflow([cfw],name='rxn_'+str(index)+str(rxn_name))
+                return FWAction(detours=newwf) #using detour allows us to inherit children from the original collect to the subsequent collects
+
+            #if irc_mode == "skip":
+            else:
+                ctask = MolecularCollect({"xyzs":xyzsout,"check_symm":True,"fw_generators": ["optimize_firework",["vibrations_firework"]],
+                        "fw_generator_dicts": [self["opt_obj_dict"],[self["vib_obj_dict"]]],
+                        "out_names": ["opt.xyz",["vib.json"]],"future_check_symms": [True,False], "label": "TS"+str(index)+"_"+rxn_name})
+                cfw = Firework([ctask],name="TS"+str(index)+"_"+rxn_name+"_collect",spec={"_allow_fizzled_parents": True, "_priority": 5})
+                newwf = Workflow([cfw],name='rxn_'+str(index)+str(rxn_name))
+                return FWAction(detours=newwf) #using detour allows us to inherit children from the original collect to the subsequent collects
         else:
             return FWAction()