Added barrier height errors and also changed the energy profile figur…

…e to contain the barrier numbers.

src/calcs.py

@@ -1,3 +1,5 @@
+import os
+import torch
 from mace.calculators import MACECalculator
 from newtonnet.utils.ase_interface import MLAseCalculator
 from ase import Atoms
@@ -16,24 +18,31 @@ def calc_mace():
 
 
 def calc():
+    use_cuda = torch.cuda.is_available()
+    device = 'cuda' if use_cuda else 'cpu'
     mlcalculator = MLAseCalculator(
         model_path=[
-                    '/global/u2/k/kumaranu/Documents/NewtonNet/example/predict/training_52/models/best_model_state.tar',
+                    os.path.join(
+                        os.path.expanduser("~"),
+                        'Documents/NewtonNet/example/predict/training_52/models/best_model_state.tar'
+                    ),
                     #'/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_53/models/best_model_state.tar',
                     #'/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_54/models/best_model_state.tar',
                     #'/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_55/models/best_model_state.tar',
                    ],    # path to model file, str or list of str
         settings_path=[
-                       '/global/u2/k/kumaranu/Documents/NewtonNet/example/predict/training_52/run_scripts/config0.yml',
+                       os.path.join(
+                           os.path.expanduser("~"),
+                           'Documents/NewtonNet/example/predict/training_52/run_scripts/config0.yml'
+                       ),
                        #'/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_53/run_scripts/config2.yml',
                        #'/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_54/run_scripts/config1.yml',
                        #'/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_55/run_scripts/config3.yml',
                       ],    # path to configuration file, str or list of str
         hess_method=None,    # method to calculate hessians. 'autograd', 'fwd_diff', 'cnt_diff', or None (default: 'autograd')
         # hess_precision=1e-5,    # hessian gradient calculation precision for 'fwd_diff' and 'cnt_diff', ignored otherwise (default: None)
         disagreement='std',    # method to calculate disagreement among models. 'std', 'std_outlierremoval', 'range':, 'values', or None (default: 'std')
-        device='cuda'   # 'cpu' or list of cuda
-        #device='cpu'   # 'cpu' or list of cuda
+        device=device
     )
     return mlcalculator
 

src/check_graph_iso.py

@@ -5,6 +5,7 @@
 from pymatgen.analysis.graphs import MoleculeGraph
 from pymatgen.analysis.local_env import OpenBabelNN
 from networkx.algorithms.graph_hashing import weisfeiler_lehman_graph_hash
+from plot2 import gen_energy_profile
 
 
 def add_specie_suffix(graph):
@@ -70,7 +71,10 @@ def compare_mols(molecule_1, molecule_2) -> bool:
     return graph_1_hash == graph_2_hash
 
 
-def check_graph_iso(ref_dir):
+def check_graph_iso(
+        ref_dir,
+        threshold=0.05,
+):
     num_directories = len(glob.glob(ref_dir + '/???/'))
     for i in range(num_directories):
         logdir = os.path.join(ref_dir, f'{i:03}')
@@ -94,16 +98,46 @@ def check_graph_iso(ref_dir):
             traj_pymat_mol2 = AseAtomsAdaptor.get_molecule(pdt_traj_atoms[-1])
             bool_product = compare_mols(traj_pymat_mol1, traj_pymat_mol2)
             if bool_reactant and bool_product:
-                continue
+                geodesic_xyz = os.path.join(
+                    ref_dir,
+                    f'{i:03}',
+                    'geodesic_path.xyz',
+                )
+                if not os.path.exists(geodesic_xyz):
+                    print(f"Geodesic's file could not be found for index {i:03}")
+                    continue
+                neb_path_xyz = os.path.join(
+                    ref_dir,
+                    f'{i:03}',
+                    'optimized_path_aseneb_NEBOptimizer_None.xyz'
+                )
+                if not os.path.exists(neb_path_xyz):
+                    print(f"NEB's final path file could not be found for index {i:03}")
+                    continue
+                ts_xyz = os.path.join(
+                    ref_dir,
+                    f'{i:03}',
+                    'TS.xyz',
+                )
+                if not os.path.exists(ts_xyz):
+                    print(f"TS's file could not be found for index {i:03}")
+                    continue
+                barrier_err = gen_energy_profile(
+                    geodesic_xyz,
+                    neb_path_xyz,
+                    ts_xyz,
+                )
+                if abs(barrier_err) > threshold:
+                    print(f'barrier_err: {barrier_err:.3f} eV for index: {i:03}')
             else:
                 print(f'bool_reactant: {bool_reactant}, bool_product: {bool_product} for index {i:03}')
         else:
             print(f'Path could not be found for index {i:03}.')
 
 
 if __name__ == '__main__':
-    # ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_nn_inputs'
-    #ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_nn_inputs1'
-    ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_rgd1_inputs'
+    ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_nn_inputs'
+    # ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_nn_inputs1'
+    # ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_rgd1_inputs'
     # ref_dir = '/home/kumaranu/Downloads/neb_nn_inputs'
     check_graph_iso(ref_dir)

src/plot2.py

@@ -1,8 +1,10 @@
 import os
 import sys
 import numpy as np
-import matplotlib.pyplot as plt
+from calcs import calc
 from ase.io import read
+import matplotlib.pyplot as plt
+
 
 def plot_energy_from_xyz(xyz_file):
     # Read the XYZ file
@@ -20,7 +22,11 @@ def plot_energy_from_xyz(xyz_file):
     # Convert lists to numpy arrays for plotting
     energies = np.array(energies)
     geometry_indices = np.array(geometry_indices)
-
+
+    # Find the index and energy value of the maximum energy point
+    max_energy_idx = np.argmax(energies)
+    max_energy = energies[max_energy_idx]
+
     # Plot energies
     plt.plot(
              geometry_indices,
@@ -29,30 +35,71 @@ def plot_energy_from_xyz(xyz_file):
              markersize=3,
              linewidth=0.5,
              linestyle='-',
-             label=os.path.basename(xyz_file))
+             label=os.path.basename(xyz_file),
+    )
+    return energies[0], max_energy, energies[-1]
+
+
+def gen_energy_profile(
+        geodesic_xyz,
+        neb_path_xyz,
+        ts_xyz,
+):
+    # Plot energies from the first XYZ file
+    _, _, _ = plot_energy_from_xyz(geodesic_xyz)
+
+    # Plot energies from the second XYZ file
+    e_r, e_highest, e_p = plot_energy_from_xyz(neb_path_xyz)
 
-# Check if two arguments are provided
-if len(sys.argv) != 3:
-    print("Usage: python script.py <xyz_file1> <xyz_file2>")
-    sys.exit(1)
+    ts = read(ts_xyz)
+    mlcalculator = calc()
+    mlcalculator.calculate(ts)
+    e_ts = mlcalculator.results['energy']
 
-# Plot energies from the first XYZ file
-plot_energy_from_xyz(sys.argv[1])
+    # Set plot labels and title
+    plt.xlabel('Geometry Index')
+    plt.ylabel('Energy (eV)')
+    plt.title(f'neb barrier f: {e_highest-e_r:.3f} eV, '
+              f'Ref. barrier f: {e_ts-e_r:.3f} eV\n'
+              f'neb barrier r: {e_highest-e_p:.3f} eV, '
+              f'Ref. barrier r: {e_ts-e_p:.3f} eV')
+    plt.legend(['geodesic', 'neb'])
+    plt.grid(True)
 
-# Plot energies from the second XYZ file
-plot_energy_from_xyz(sys.argv[2])
+    # Save the plot in the same directory as the first XYZ file
+    output_file = os.path.splitext(geodesic_xyz)[0] + '_energy_plot.png'
+    plt.savefig(output_file)
 
-# Set plot labels and title
-plt.xlabel('Geometry Index')
-plt.ylabel('Energy (eV)')
-plt.title('Energy vs. Geometry Index')
-plt.legend()
-plt.grid(True)
+    # Show the plot
+    plt.show()
+    return e_highest-e_ts
 
-# Save the plot in the same directory as the first XYZ file
-output_file = os.path.splitext(sys.argv[1])[0] + '_energy_plot.png'
-plt.savefig(output_file)
 
-# Show the plot
-plt.show()
+if __name__ == '__main__':
+    # Check if two arguments are provided
+    # if len(sys.argv) != 3:
+    #     print("Usage: python script.py <xyz_file1> <xyz_file2>")
+    #     sys.exit(1)
+    ref_dir = '/home/kumaranu/Downloads/neb_nn_inputs'
+    index = '001'
+    geodesic_xyz = os.path.join(
+        ref_dir,
+        index,
+        'geodesic_path.xyz',
+    )
+    neb_path_xyz = os.path.join(
+        ref_dir,
+        index,
+        'optimized_path_aseneb_NEBOptimizer_None.xyz'
+    )
+    ts_xyz = os.path.join(
+        ref_dir,
+        index,
+        'TS.xyz',
+    )
 
+    barrier_err = gen_energy_profile(
+        geodesic_xyz,
+        neb_path_xyz,
+        ts_xyz,
+    )