Add files via upload

- Additional implementation of ASE interface for NEB method and use of NNP

Author: ss0832

multioptpy/ase_calculation_tools.py

@@ -4,8 +4,11 @@
 import re
 import numpy as np
 
-from ase import Atoms
-from ase.vibrations import Vibrations
+try:
+    from ase import Atoms
+    from ase.vibrations import Vibrations
+except ImportError:
+    print("ASE is not installed. Please install ASE to use this module.")
 
 from calc_tools import Calculationtools
 from parameter import UnitValueLib, number_element
@@ -80,9 +83,9 @@ def calc_exact_hess(self, atom_obj, positions, element_list):
         vib.clean()  
         exact_hess = exact_hess / self.hartree2eV * (self.bohr2angstroms ** 2)
         if type(element_list[0]) is str:
-            exact_hess = Calculationtools().project_out_hess_tr_and_rot(exact_hess, element_list, positions)
+            exact_hess = Calculationtools().project_out_hess_tr_and_rot_for_coord(exact_hess, element_list, positions)
         else:
-            exact_hess = Calculationtools().project_out_hess_tr_and_rot(exact_hess, [number_element(elem_num) for elem_num in element_list], positions)
+            exact_hess = Calculationtools().project_out_hess_tr_and_rot_for_coord(exact_hess, [number_element(elem_num) for elem_num in element_list], positions)
         self.Model_hess = exact_hess
         return exact_hess
 
@@ -99,16 +102,7 @@ def single_point(self, file_directory, element_list, iter, electric_charge_and_m
         for num, input_file in enumerate(file_list):
             try:
                 if geom_num_list is None:
-                    print("\n",input_file,"\n")
-
-                    with open(input_file,"r") as f:
-                        input_data = f.readlines()
-                    
-                    positions = []
-                  
-                    for word in input_data[2:]:
-                        positions.append(word.split()[1:4])
-                
+                    positions, _, electric_charge_and_multiplicity = xyz2list(input_file, electric_charge_and_multiplicity)
                 else:
                     positions = geom_num_list        
 
@@ -152,7 +146,9 @@ def single_point(self, file_directory, element_list, iter, electric_charge_and_m
                         self.FUNCTIONAL,
                         self.BASIS_SET
                     )
-                    
+                elif self.software_type == "uma-s-1":  # Neural Network Potential
+                    atom_obj = use_FAIRCHEMNNP(atom_obj, electric_charge_and_multiplicity, self.software_path_dict["uma-s-1"])
+
                 elif self.software_type == "mace_mp":  # Neural Network Potential
                     atom_obj = use_MACE_MP(atom_obj, electric_charge_and_multiplicity)
 
@@ -301,6 +297,22 @@ def use_MACE_MP(atom_obj, electric_charge_and_multiplicity):
     atom_obj.calc = macemp
     return atom_obj
 
+def use_FAIRCHEMNNP(atom_obj, electric_charge_and_multiplicity, fairchem_path):
+    try:
+        from fairchem.core import FAIRChemCalculator
+        from fairchem.core.units.mlip_unit import load_predict_unit
+    except ImportError:
+        raise ImportError("FAIRChem.core modules not found")
+    # Load the prediction unit
+    predict_unit = load_predict_unit(path=fairchem_path, device="cpu")
+    
+    # Set up the FAIRChem calculator
+    fairchem_calc = FAIRChemCalculator(predict_unit=predict_unit, task_name="omol")
+    atom_obj.info = {"charge": int(electric_charge_and_multiplicity[0]),
+                     "spin": int(electric_charge_and_multiplicity[1])}
+    atom_obj.calc = fairchem_calc
+    return atom_obj
+
 def use_MACE_OFF(atom_obj, electric_charge_and_multiplicity):
     """Set up and return ASE atoms object with MACE_OFF calculator.
     
@@ -373,25 +385,11 @@ class ASEEngine(CalculationEngine):
     This engine uses the Atomic Simulation Environment (ASE) to interface with various
     quantum chemistry software packages like GAMESSUS, NWChem, Gaussian, ORCA, MACE, and MOPAC.
     """
-    def __init__(self, **kwarg):
+    def __init__(self, **kwargs):
         UVL = UnitValueLib()
 
         self.bohr2angstroms = UVL.bohr2angstroms
         self.hartree2eV = UVL.hartree2eV
-        
-        self.START_FILE = kwarg["START_FILE"]
-        self.N_THREAD = kwarg["N_THREAD"]
-        self.SET_MEMORY = kwarg["SET_MEMORY"]
-        self.FUNCTIONAL = kwarg["FUNCTIONAL"]
-        self.BASIS_SET = kwarg["BASIS_SET"]
-        self.FC_COUNT = kwarg["FC_COUNT"]
-        self.BPA_FOLDER_DIRECTORY = kwarg["BPA_FOLDER_DIRECTORY"]
-        self.Model_hess = kwarg["Model_hess"]
-        self.unrestrict = kwarg["unrestrict"]
-        self.software_type = kwarg["software_type"]
-        self.hessian_flag = False
-        self.software_path_dict = read_software_path()
-
 
     def calculate(self, file_directory, optimize_num, pre_total_velocity, config):
         """Calculate energy and gradients using ASE and the configured software.
@@ -426,7 +424,7 @@ def calculate(self, file_directory, optimize_num, pre_total_velocity, config):
         geometry_list_tmp, element_list, _ = xyz2list(file_list[0], None)
 
         hess_count = 0
-        software_type = config.software_type
+        software_type = config.othersoft
         software_path_dict = read_software_path()
 
         for num, input_file in enumerate(file_list):
@@ -472,25 +470,29 @@ def calculate(self, file_directory, optimize_num, pre_total_velocity, config):
                         config.FUNCTIONAL,
                         config.BASIS_SET
                     )
+                elif software_type == "uma-s-1":  # Neural Network Potential
+                    atom_obj = use_FAIRCHEMNNP(atom_obj, electric_charge_and_multiplicity, software_path_dict["uma-s-1"])
                 elif software_type == "mace_mp":  # Neural Network Potential
                     atom_obj = use_MACE_MP(atom_obj, electric_charge_and_multiplicity)
                 elif software_type == "mace_off":  # Neural Network Potential
                     atom_obj = use_MACE_OFF(atom_obj, electric_charge_and_multiplicity)
+
+
                 elif software_type == "mopac":
                     atom_obj = use_MOPAC(atom_obj, electric_charge_and_multiplicity, input_file)
                 else:
                     print(f"Software type '{software_type}' is not available...")
                     raise ValueError(f"Unsupported software type: {software_type}")
 
                 # Calculate energy and forces
-                e = atom_obj.get_potential_energy(apply_constraint=False) / config.hartree2eV  # eV to hartree
-                g = -1 * atom_obj.get_forces(apply_constraint=False) * config.bohr2angstroms / config.hartree2eV  # eV/ang. to a.u.
+                e = atom_obj.get_potential_energy(apply_constraint=False) / self.hartree2eV  # eV to hartree
+                g = -1 * atom_obj.get_forces(apply_constraint=False) * self.bohr2angstroms / self.hartree2eV  # eV/ang. to a.u.
 
                 # Store results
                 energy_list.append(e)
                 gradient_list.append(g)
                 gradient_norm_list.append(np.sqrt(np.linalg.norm(g)**2/(len(g)*3)))  # RMS
-                geometry_num_list.append(positions / config.bohr2angstroms)  # Convert to Bohr
+                geometry_num_list.append(positions / self.bohr2angstroms)  # Convert to Bohr
                 num_list.append(num)
 
                 # Handle hessian calculation if needed
@@ -505,14 +507,12 @@ def calculate(self, file_directory, optimize_num, pre_total_velocity, config):
                     vib.clean()  
 
                     # Convert hessian units
-                    exact_hess = exact_hess / config.hartree2eV * (config.bohr2angstroms ** 2)
-                    
+                    exact_hess = exact_hess / self.hartree2eV * (self.bohr2angstroms ** 2)
+
                     # Project out translational and rotational modes
                     calc_tools = Calculationtools()
-                    element_number_list = [number_element(elem) for elem in element_list]
-                    exact_hess = calc_tools.project_out_hess_tr_and_rot(
-                        exact_hess, element_number_list, positions / config.bohr2angstroms
-                    )
+                    exact_hess = calc_tools.project_out_hess_tr_and_rot_for_coord(exact_hess, element_list, positions)
+                       
 
                     # Save hessian
                     np.save(os.path.join(config.NEB_FOLDER_DIRECTORY, f"tmp_hessian_{hess_count}.npy"), exact_hess)

multioptpy/ieip.py

@@ -93,10 +93,15 @@ def __init__(self, args):
             self.ECP = ""
 
 
-
+        self.othersoft = args.othersoft
         self.basic_set_and_function = args.functional+"/"+args.basisset
         self.force_data = force_data_parser(args)
-        if args.usextb == "None" and args.usedxtb == "None":
+        if self.othersoft != "None":
+            self.iEIP_FOLDER_DIRECTORY = args.INPUT+"_iEIP_"+args.othersoft+"_"+str(datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S_%f")[:-2])+"/"
+
+        
+        
+        elif args.usextb == "None" and args.usedxtb == "None":
             self.iEIP_FOLDER_DIRECTORY = args.INPUT+"_iEIP_"+self.basic_set_and_function.replace("/","_")+"_"+str(datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S_%f")[:-2])+"/"
         else:
             if args.usedxtb != "None":
@@ -713,7 +718,9 @@ def close_target_force(self, L, Lt, geom_num_list_1, geom_num_list_2):
 
 
     def optimize(self):
-        if self.args.pyscf:
+        if self.othersoft != "None":
+            from ase_calculation_tools import Calculation
+        elif self.args.pyscf:
             from pyscf_calculation_tools import Calculation
         elif self.args.usextb != "None" and self.args.usedxtb == "None":
             from tblite_calculation_tools import Calculation
@@ -762,8 +769,9 @@ def optimize(self):
                          spin_multiplicity = self.spin_multiplicity[i] or electric_charge_and_multiplicity_list[i][1],
                          excited_state = self.excite_state_list[i],
                          dft_grid=self.dft_grid,
-                         ECP = self.ECP))
-            
+                         ECP = self.ECP,
+                         software_type=self.othersoft,))
+
             SP_list[i].cpcm_solv_model = self.cpcm_solv_model
             SP_list[i].alpb_solv_model = self.alpb_solv_model
 

multioptpy/interface.py

@@ -58,7 +58,8 @@ def ieipparser(parser):
     parser.add_argument("-cpcm", "--cpcm_solv_model",  type=str, default=None, help='use CPCM solvent model for xTB (Defalut setting is not using this model.) (ex.) water')
     parser.add_argument("-alpb", "--alpb_solv_model",  type=str, default=None, help='use ALPB solvent model for xTB (Defalut setting is not using this model.) (ex.) water')#ref.: J. Chem. Theory Comput. 2021, 17, 7, 4250–4261 https://doi.org/10.1021/acs.jctc.1c00471
     parser.add_argument("-grid", "--dft_grid", type=int, default=3, help="fineness of grid for DFT calculation (default: 3 (0~9))")
-
+    parser.add_argument("-os", "--othersoft",  type=str, default="None", help='use other QM software. default is not using other QM software. (require python module, ASE (Atomic Simulation Environment)) (ex.) orca, gaussian, gamessus, mace_mp etc.')
+    
     args = parser.parse_args()#model_function_mode
     args.fix_atoms = []
     args.gradient_fix_atoms = []
@@ -197,6 +198,7 @@ def nebparser(parser):
     parser.add_argument("-mem", "--SET_MEMORY",  type=str, default='1GB', help='use mem(ex. 1GB)')
     parser.add_argument("-cineb", "--apply_CI_NEB",  type=int, default='99999', help='apply CI_NEB method')
     parser.add_argument("-sd", "--steepest_descent",  type=int, default='99999', help='apply steepest_descent method')
+
 
     class CGAction(argparse.Action):
         def __call__(self, parser, namespace, values, option_string=None):
@@ -235,6 +237,7 @@ def __call__(self, parser, namespace, values, option_string=None):
 
     parser.add_argument('-modelhess','--use_model_hessian', nargs='?', help="use model hessian. (Default: not using model hessian If you specify only option, Improved Lindh + Grimme's D3 dispersion model hessian is used.) (ex. lindh, gfnff, gfn0xtb, fischer, fischerd3, fischerd4, schlegel, swart, lindh2007, lindh2007d3, lindh2007d4)", action=ModelhessAction, default=None)
     parser.add_argument("-mfc", "--calc_model_hess",  type=int, default=50, help='calculate model hessian per steps (ex.) [steps per one hess calculation]')
+    parser.add_argument("-os", "--othersoft",  type=str, default="None", help='use other QM software. default is not using other QM software. (require python module, ASE (Atomic Simulation Environment)) (ex.) orca, gaussian, gamessus, mace_mp etc.')
 
     parser = parser_for_biasforce(parser)
     args = parser.parse_args()

multioptpy/neb.py

@@ -67,6 +67,7 @@
 from pyscf_calculation_tools import PySCFEngine
 from psi4_calculation_tools import Psi4Engine
 from dxtb_calculation_tools import DXTBEngine
+from ase_calculation_tools import ASEEngine
 
 class NEBConfig:
     """Configuration management class for NEB calculations"""
@@ -131,6 +132,7 @@ def __init__(self, args):
         self.dft_grid = int(args.dft_grid)
         self.spring_constant_k = 0.01
         self.force_const_for_cineb = 0.01
+        self.othersoft = args.othersoft
 
         # FIRE method parameters
         self.FIRE_dt = 0.1
@@ -213,8 +215,11 @@ def make_neb_work_directory(self, input_file):
             tmp_name = input_file 
 
         timestamp = str(datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S_%f")[:-2])
-        
-        if self.usextb == "None" and self.usedxtb == "None":
+        if self.othersoft != "None":
+            return tmp_name + "_NEB_" + self.othersoft + "_" + timestamp + "/"
+
+
+        elif self.usextb == "None" and self.usedxtb == "None":
             return tmp_name + "_NEB_" + self.basic_set_and_function.replace("/", "_") + "_" + timestamp + "/"
         else:
             if self.usextb != "None":
@@ -229,7 +234,9 @@ class CalculationEngineFactory:
     @staticmethod
     def create_engine(config):
         """Create appropriate calculation engine based on configuration"""
-        if config.usextb != "None":
+        if config.othersoft != "None":
+            return ASEEngine()    
+        elif config.usextb != "None":
             return TBLiteEngine()
         elif config.usedxtb != "None":
             return DXTBEngine()

multioptpy/optimization.py

@@ -141,7 +141,7 @@ def _make_init_directory(self, file):
         date = datetime.datetime.now().strftime("%Y_%m_%d")
         base_dir = f"{date}/{self.START_FILE[:-4]}_OPT_"
 
-        if self.args.othersoft != "None":
+        if self.othersoft != "None":
             self.BPA_FOLDER_DIRECTORY = f"{base_dir}ASE_{timestamp}/"
         elif self.args.usextb == "None" and self.args.usedxtb == "None":
             self.BPA_FOLDER_DIRECTORY = f"{base_dir}{self.FUNCTIONAL}_{self.BASIS_SET}_{timestamp}/"