added more files from the analysis scripts on nersc

src/analyze_all_data.py

@@ -0,0 +1,14 @@
+import numpy as np
+
+f = open('all_data.txt', 'r')
+INF = f.readlines()
+for line in INF:
+    # print(line.split())
+    tmp_data_list = line.split()
+    index = tmp_data_list[0].strip(':')
+    # print(index)
+    ts_e_err = float(tmp_data_list[4])
+    if abs(ts_e_err) > 0.15:
+        print(index, ts_e_err)
+f.close()
+

src/bb.py

@@ -0,0 +1,112 @@
+import os
+import sys
+import json
+import logging
+import argparse
+import shutil, errno
+from ase.atoms import Atoms
+from tabulate import tabulate
+from ase.io import read, write
+from ase.io.jsonio import decode
+from maggma.stores import MongoStore
+import numpy as np
+
+
+def get_docs(launchpad_file, query, collections_name='fireworks'):
+    tasks_store = MongoStore.from_launchpad_file(launchpad_file, collections_name)
+    tasks_store.connect()
+    docs = list(tasks_store.query(query))
+    return docs
+
+
+def atomic_number_to_symbol(atomic_number):
+    periodic_table = {
+        1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B',
+        6: 'C', 7: 'N', 8: 'O', 9: 'F', 10: 'Ne',
+    }
+    return periodic_table.get(atomic_number, 'X')  # 'X' for unknown element
+
+
+def atoms_from_json(atoms_json):
+    atoms_dict = json.loads(atoms_json)
+
+    atomic_numbers = atoms_dict['numbers']['__ndarray__'][2]
+    positions = atoms_dict['positions']['__ndarray__'][2]
+
+    natoms = int(len(positions)/3)
+    coords = np.reshape(positions, (natoms, 3))
+
+    atomic_symbols = [atomic_number_to_symbol(i) for i in atomic_numbers]
+    return Atoms(symbols=atomic_symbols, positions=coords)
+
+
+def download_data(docs, base_dir):
+    for doc in docs:
+        tag = doc['metadata']['tag']
+        index = tag.split('-')[-1].split('_')[0]
+        # print('tag:', tag)
+        if tag.startswith('TS0'):
+            name = 'TS'
+            structure = atoms_from_json(doc['output']['atoms']['atoms_json'])
+            '''
+            print(
+                  name,
+                  structure,
+                  structure.get_positions()
+                 )
+            '''
+        elif tag.startswith('irc'):
+            name = tag.split('-')[1]
+            structure = atoms_from_json(doc['output']['atoms']['atoms_json'])
+            '''
+            print(
+                  name,
+                  structure,
+                  structure.get_positions()
+                 )
+            '''
+        else:
+            continue
+
+        output_dir = os.path.join(base_dir, f"{index}")
+
+        os.makedirs(output_dir, exist_ok=True)
+
+        output_file = os.path.join(output_dir, f'{name}.xyz')
+        write(output_file, structure, format='xyz')
+
+
+def main():
+    tag = 'nov15'
+    noise_level = '00'
+
+    query0 = {
+              "metadata.tags.class": tag,
+               "metadata.tag": {
+                        '$regex': f'.*0.*_noise{noise_level}'
+                       }
+             }
+
+    launchPad_file = os.path.join(os.environ["HOME"],"fw_config/my_launchpad.yaml")
+    collections = 'quacc_results0'
+
+    docs0 = get_docs(launchPad_file, query0, collections_name=collections)
+    download_data(docs0, "/global/scratch/users/kumaranu/saved_files/neb_nn_inputs/")
+
+    '''
+    for doc in docs0:
+        # print(doc.keys())
+        print(doc['metadata']['tag'])
+        #a = 'irc-reverse0-021_noise00'
+        index = a.split('-')[-1].split('_')[0]
+        if a.split('-')[1] == 'reverse0':
+            print('reverse')
+
+        #os.system(f'cp -rL {src} {dst}')
+    print(len(docs0))
+    '''
+
+
+if __name__ == "__main__":
+    main()
+

src/calc_properties.py

@@ -0,0 +1,37 @@
+from ase.io import read, write
+from calcs import calc
+
+# Read molecular configurations from XYZ file
+atoms_list = read('/global/scratch/users/kumaranu/saved_files/neb_nn_inputs/263/output.xyz', index=':')
+
+# Create a calculator instance (e.g., EMT)
+mlcalculator = calc()
+
+# Iterate over each molecular configuration
+for atoms in atoms_list:
+    mlcalculator.calculate(atoms)
+
+    # Calculate energies and forces
+    energy = mlcalculator.results['energy']
+    forces = mlcalculator.results['forces']
+
+    # Store energies and forces in the atoms object
+    atoms.info['energy'] = energy
+    atoms.arrays['forces'] = forces
+
+# Write the updated atoms object with energies and forces to a new XYZ file
+write('/global/scratch/users/kumaranu/saved_files/neb_nn_inputs/263/output1.xyz', atoms_list)
+
+
+'''
+h2 = Atoms(numbers=[1, 1], positions=[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+mlcalculator.calculate(h2)
+
+print(mlcalculator.results['energy'])    # mean of calculated molecular energies, shape (1,)
+print(mlcalculator.results['forces'])    # mean of calculated atomic forces, shape (n_atom, 3)
+print(mlcalculator.results['hessian'])    # mean of calculated atomic Hessian, shape (n_atom, 3, n_atom, 3)
+print(mlcalculator.results['energy_disagreement'])    # disagreement among calculated molecular 
+'''
+
+
+

src/calcs.py

@@ -18,13 +18,13 @@ def calc_mace():
 def calc():
     mlcalculator = MLAseCalculator(
         model_path=[
-                    '/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_52/models/best_model_state.tar',
+                    '/global/u2/k/kumaranu/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/home/users/kumaranu/Documents/NewtonNet/example/predict/training_52/run_scripts/config0.yml',
+                       '/global/u2/k/kumaranu/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',

src/calcs1.py

@@ -0,0 +1,43 @@
+from mace.calculators import MACECalculator
+from newtonnet.utils.ase_interface import MLAseCalculator
+from ase import Atoms
+
+
+def calc_mace():
+
+    mlcalculator = MACECalculator(
+        model_paths='/global/home/users/kumaranu/Documents/gpu_jobs/MACE_model.model',
+        #model_paths='/global/home/users/kumaranu/Documents/mep_tests/src/MACE_model_cpu.model',
+        device='cuda',
+        #device='cpu',
+        default_dtype="float64",
+    )
+    return mlcalculator
+
+
+def calc():
+    mlcalculator = MLAseCalculator(
+        model_path=[
+                    '/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_52/models/best_model_state.tar',
+                   ],    # path to model file, str or list of str
+        settings_path=[
+                       '/global/home/users/kumaranu/Documents/NewtonNet/example/predict/training_52/run_scripts/config0.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
+    )
+    return mlcalculator
+
+'''
+h2 = Atoms(numbers=[1, 1], positions=[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+mlcalculator.calculate(h2)
+
+print(mlcalculator.results['energy'])    # mean of calculated molecular energies, shape (1,)
+print(mlcalculator.results['forces'])    # mean of calculated atomic forces, shape (n_atom, 3)
+print(mlcalculator.results['hessian'])    # mean of calculated atomic Hessian, shape (n_atom, 3, n_atom, 3)
+print(mlcalculator.results['energy_disagreement'])    # disagreement among calculated molecular 
+'''
+

src/check.py

@@ -0,0 +1,18 @@
+import glob, os, sys
+
+logdir = '/global/scratch/users/kumaranu/saved_files/neb_nn_inputs'
+
+os.chdir(logdir)
+for i in range(265):
+    if os.path.exists(f'{i:03}'):
+        os.chdir(f'{i:03}')
+        if os.path.exists('neb_band_aseneb_NEBOptimizer_None.txt'):
+            INF = open('neb_band_aseneb_NEBOptimizer_None.txt', 'r').readlines()
+            #if int(INF[-1].split()[1]) > 200:
+            print(f'{i:03}: {INF[-1].split()[1]}')
+        else:
+            print(f'{i:03}: No NEB file.')
+        os.chdir('../')
+    else:
+        print(f'{i:03} does not exist.')
+

src/check_graph_iso.py

@@ -0,0 +1,18 @@
+import glob, os
+from ase.io import read
+
+ref_dir = '/global/cfs/cdirs/m2834/kumaranu/neb_nn_inputs'
+
+for i in range(265):
+    logdir = os.path.join(ref_dir, f'{i:03}')
+    if os.path.exists(logdir):
+        reactant_traj_file = os.path.join(logdir, 'reactant_opt.traj')
+        if os.path.exists(reactant_traj_file):
+            traj_atoms = read(reactant_traj_file, index=':', format="traj")
+        else:
+            print(f"Reactant's traj file could not be found for index {i:03}")
+        #for atoms in traj_atoms:
+        #    print(i, atoms)
+    else:
+        print(f'Path could not be found for index {i:03}.')
+

src/check_neb_results.py

@@ -0,0 +1,82 @@
+import glob, os, sys
+from ase.io import read
+from calcs1 import calc
+import numpy as np
+from ase import Atoms
+from multiprocessing import Pool
+
+
+def align_geom(refgeom, geom):
+    """Find translation/rotation that moves a given geometry to maximally overlap
+    with a reference geometry. Implemented with Kabsch algorithm.
+
+    Args:
+        refgeom:    The reference geometry to be rotated to
+        geom:       The geometry to be rotated and shifted
+
+    Returns:
+        RMSD:       Root-mean-squared difference between the rotated geometry
+                    and the reference
+        new_geom:   The rotated geometry that maximumally overal with the reference
+    """
+    center = np.mean(refgeom, axis=0)   # Find the geometric center
+    ref2 = refgeom - center
+    geom2 = geom - np.mean(geom, axis=0)
+    cov = np.dot(geom2.T, ref2)
+    v, sv, w = np.linalg.svd(cov)
+
+    if np.linalg.det(v) * np.linalg.det(w) < 0:
+        sv[-1] = -sv[-1]
+        v[:, -1] = -v[:, -1]
+    u = np.dot(v, w)
+    new_geom = np.dot(geom2, u) + center
+    rmsd = np.sqrt(np.mean((new_geom - refgeom) ** 2))
+    return rmsd, new_geom
+
+
+def process_file(i):
+    try:
+        if os.path.exists(f'{i:03}'):
+            os.chdir(f'{i:03}')
+            if os.path.exists('neb_band_aseneb_NEBOptimizer_None.txt'):
+                ase_traj = read('optimized_path_aseneb_NEBOptimizer_None.xyz', ':')
+                energies = [i.get_potential_energy() for i in ase_traj]
+                ts_index = np.argmax(energies)
+
+                mlcalculator = calc()
+                ts_atoms = read('TS.xyz')
+                mlcalculator.calculate(ts_atoms)
+                energy_ts = mlcalculator.results['energy']
+                rmsd, new_geom = align_geom(ase_traj[ts_index].get_positions(), ts_atoms.get_positions())
+
+                fwd_barrier_neb = energies[ts_index] - energies[0]
+                rev_barrier_neb = energies[ts_index] - energies[-1]
+                ref_fwd_barrier = energy_ts - energies[0]
+                ref_rev_barrier = energy_ts - energies[-1]
+
+                ts_e_err = energy_ts - energies[ts_index]
+                fwd_barrier_err = fwd_barrier_neb - ref_fwd_barrier
+                rev_barrier_err = rev_barrier_neb - ref_rev_barrier
+                os.chdir('../')
+                return f'{i:03}: TS energy err: {ts_e_err:.6f} fwd barrier err: {fwd_barrier_err:.6f} Rev barrier err: {rev_barrier_err:.6f} Geometry RMSD: {rmsd:.6f}'
+            else:
+                os.chdir('../')
+    except Exception as e:
+        return f'Error processing file {i}: {str(e)}'
+
+if __name__ == '__main__':
+
+    logdir = '/global/scratch/users/kumaranu/saved_files/neb_nn_inputs'
+    os.chdir(logdir)
+
+    with Pool(5) as pool:
+        results = pool.map(process_file, range(265))
+
+    missing_indices = [i for i, result in enumerate(results) if result is None]
+    if missing_indices:
+        print(f'Missing results for indices: {missing_indices}')
+
+    for result in results:
+        if result is not None:
+            print(result)
+