add flag for poorly behaving RASPA2 inputs

Author: Jake Smith

mofdiff/gcmc/gcmc_wrapper.py

@@ -139,9 +139,6 @@ def calculate_unit_cells(self, forcefield_cutoff):
 
     def write_out(self, output_path):
         with open(output_path, "w") as log_file:
-            # log_file.write("input file: " + self.sorbent_file + "\n")
-            # log_file.write("helium_void_fraction: " + str(self.helium_void_fraction) + "\n")
-            # log_file.write("helium_void_fraction: " + str(self.helium_void_fraction) + "\n")
             log_file.write(self.raspa_output)
 
 # assigns charges to the atoms in the simulation file using the MEPO Qeq charge equilibration method
@@ -192,21 +189,18 @@ def run_gcmc_simulation(
     initialization_cycles=0,
     equilibration_cycles=2000,
     production_cycles=2000,
-    forcefield="UFF-TraPPe",
+    forcefield="UFF",
     forcefield_cutoff=12,
     molecule_definitions="TraPPE",
     unit_cells=[0, 0, 0],
     cleanup=False,
+    rewrite_raspa_input=False,
 ):
     # copy cif file into parent RASPA folder
     shutil.copy(simulation.sorbent_file, raspa_path + "/share/raspa/structures/cif/")
     workdir = simulation.rundir / "raspa_output" / simulation.identifier
     workdir.mkdir(exist_ok=True, parents=True)
 
-    # create and enter temp directory
-    # os.mkdir(simulation.rundir / "raspa_output" / simulation.identifier)
-    # os.chdir(os.getcwd() + "/temp/raspa_output/" + simulation.identifier)
-
     sorbent_file = ".".join(simulation.sorbent_file.split("/")[-1].split(".")[:-1])
 
     # calculate number of unit cells needed if not user defined
@@ -299,8 +293,14 @@ def run_gcmc_simulation(
     # write out raspa input file
     with open(workdir / "simulation.input", "w") as raspa_input:
         raspa_input.write(simulation.raspa_config)
+    
+    # optionally rewrite input file to avoid errors with RASPA reading it
+    if rewrite_raspa_input:
+        raspa_input_path = workdir / "simulation.input"
+        subprocess.run(["mv", raspa_input_path, f"{raspa_input_path}.orig"])
+        command = f"printf '%s\\n' \"$(cat {raspa_input_path}.orig)\" > {raspa_input_path}"
+        subprocess.Popen(command, shell = True)
 
-    # os.chdir(str(workdir))
     # run raspa simulation
     subprocess.run([raspa_sim_path, "simulation.input"], cwd=str(workdir))
 
@@ -313,4 +313,4 @@ def run_gcmc_simulation(
 
     # clear temp directory
     if cleanup:
-        shutil.rmtree(str(workdir))
+        shutil.rmtree(str(workdir))

mofdiff/gcmc/simulation.py

@@ -87,7 +87,7 @@ def extract_raspa_output(raspa_output, has_N2=False):
 
 @timeout(36000)
 def working_capacity_vacuum_swing(cif_file, calc_charges=True,
-                                  rundir='./temp'):
+                                  rundir='./temp', rewrite_raspa_input=False):
     random.seed(4)
     np.random.seed(4)
     # adsorption conditions
@@ -102,7 +102,10 @@ def working_capacity_vacuum_swing(cif_file, calc_charges=True,
 
     if calc_charges:
         gcmc_wrapper.calculate_mepo_qeq_charges(adsorbed)
-    gcmc_wrapper.run_gcmc_simulation(adsorbed)
+    gcmc_wrapper.run_gcmc_simulation(
+        adsorbed,
+        rewrite_raspa_input=rewrite_raspa_input,
+    )
 
     (
         adsorbed_CO2,
@@ -124,7 +127,10 @@ def working_capacity_vacuum_swing(cif_file, calc_charges=True,
 
     if calc_charges:
         gcmc_wrapper.calculate_mepo_qeq_charges(residual)
-    gcmc_wrapper.run_gcmc_simulation(residual)
+    gcmc_wrapper.run_gcmc_simulation(
+        residual,
+        rewrite_raspa_input=rewrite_raspa_input,
+    )
 
     residual_CO2, heat_of_adsorption_CO2_363 = extract_raspa_output(
         residual.raspa_output, has_N2=False
@@ -149,4 +155,4 @@ def run_or_fail(cif_path):
         return working_capacity_vacuum_swing(cif_path)
     except Exception as e:
         print(e)
-        return None
+        return None

mofdiff/scripts/gcmc_screen.py

@@ -10,7 +10,7 @@
 from p_tqdm import p_umap
 
 
-def main(input_dir, ncpu=24):
+def main(input_dir, ncpu=24, rewrite_raspa_input=False):
     rundir = Path(input_dir) / 'gcmc'
     rundir.mkdir(exist_ok=True)
 
@@ -42,7 +42,11 @@ def compute_gcmc(ciffile, max_natom=10000):
                     info = 'atomic overlap'
                 else:
                     adsorption_info = working_capacity_vacuum_swing(
-                        str(ciffile), calc_charges=calc_charges, rundir=rundir)
+                        str(ciffile),
+                        calc_charges=calc_charges,
+                        rundir=rundir,
+                        rewrite_raspa_input=rewrite_raspa_input,
+                    )
                     info = 'success'
         except Exception as e:
             print(f'Error in {ciffile}: {e}')
@@ -58,5 +62,7 @@ def compute_gcmc(ciffile, max_natom=10000):
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
     parser.add_argument('--input', type=str)
+    parser.add_argument('--rewrite_raspa_input', action='store_true')
+    parser.set_defaults(rewrite_raspa_input=False)
     args = parser.parse_args()
-    main(args.input)
+    main(args.input, rewrite_raspa_input=args.rewrite_raspa_input)