Merge: bugfixes and improvements

ash/__init__.py

@@ -113,7 +113,7 @@
 from .interfaces.interface_MNDO import MNDOTheory
 
 from .interfaces.interface_CFour import CFourTheory, run_CFour_HLC_correction, run_CFour_DBOC_correction, convert_CFour_Molden_file
-from .interfaces.interface_xtb import xTBTheory, gxTBTheory
+from .interfaces.interface_xtb import xTBTheory, gxTBTheory,tbliteTheory
 from .interfaces.interface_DFTB import DFTBTheory
 from .interfaces.interface_PyMBE import PyMBETheory
 from .interfaces.interface_MLatom import MLatomTheory

ash/functions/functions_elstructure.py

@@ -2595,3 +2595,14 @@ def density_sensitivity_metric(fragment=None, functional="B3LYP", basis="def2-TZ
     #TODO: Option to plot difference density also
 
     return metrics_dict
+
+def boltzmann_populations(energies, temperature=298.15):
+    print("Inside boltzmann_populations function")
+    beta = 1/(ash.constants.R_gasconst_kcalK*temperature)
+
+    rel_energies=np.array([en-min(energies) for en in energies])*ash.constants.hartokcal
+    print("Relative energies (kcal/mol):", rel_energies)
+    boltzmann_factors = np.exp(-1*rel_energies * beta)
+    populations = boltzmann_factors / np.sum(boltzmann_factors)
+    print("Boltzmann populations at", temperature, "K:", populations)
+    return populations

ash/interfaces/interface_OpenMM.py

@@ -483,10 +483,18 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo
             # Create list of atomnames, used in PDB topology and XML file
             atomnames_full=[j+str(i) for i,j in enumerate(fragment.elems)]
             # Write PDB-file frag.pdb with dummy atomnames
-            write_pdbfile(fragment, outputname="frag", atomnames=atomnames_full)
+            #write_pdbfile(fragment, outputname="frag", atomnames=atomnames_full)
             # Load PDB-file and create topology
-            pdb = openmm.app.PDBFile("frag.pdb")
-            self.topology = pdb.topology
+            #pdb = openmm.app.PDBFile("frag.pdb")
+            #self.topology = pdb.topology
+
+            #Creating new 
+            #fragment.define_topology()
+            #self.topology = fragment.pdb_topology
+            from ash.modules.module_coords import define_dummy_topology
+            self.topology = define_dummy_topology(fragment.elems)
+            print("self.topology:", self.topology)
+            print("self.topology dict:", self.topology.__dict__)
 
             # Create dummy XML file
             xmlfile = write_xmlfile_nonbonded(filename="dummy.xml", resnames=["DUM"], atomnames_per_res=[atomnames_full], atomtypes_per_res=[fragment.elems],
@@ -1497,14 +1505,9 @@ def addexceptions(self, atomlist):
         import itertools
         print("Add exceptions/exclusions. Removing i-j interactions for list:", len(atomlist), "atoms")
 
-        # Has duplicates
-        # [self.nonbonded_force.addException(i,j,0, 0, 0, replace=True) for i in atomlist for j in atomlist]
-        # https://stackoverflow.com/questions/942543/operation-on-every-pair-of-element-in-a-list
-        # [self.nonbonded_force.addException(i,j,0, 0, 0, replace=True) for i,j in itertools.combinations(atomlist, r=2)]
         numexceptions = 0
         numexclusions = 0
         printdebug("self.system.getForces() ", self.system.getForces())
-        # print("self.nonbonded_force:", self.nonbonded_force)
 
         for force in self.system.getForces():
             printdebug("force:", force)
@@ -1847,70 +1850,19 @@ def delete_exceptions(self, atomlist):
         for force in self.system.getForces():
             if isinstance(force, openmm.NonbondedForce):
                 for exc in range(force.getNumExceptions()):
-                    # print(force.getExceptionParameters(exc))
+                    #print(force.getExceptionParameters(exc))
                     # force.getExceptionParameters(exc)
                     p1, p2, chargeprod, sigmaij, epsilonij = force.getExceptionParameters(exc)
                     if p1 in atomlist or p2 in atomlist:
-                        # print("p1: {} and p2: {}".format(p1,p2))
-                        # print("chargeprod:", chargeprod)
-                        # print("sigmaij:", sigmaij)
-                        # print("epsilonij:", epsilonij)
+                        #print("p1: {} and p2: {}".format(p1,p2))
+                        #print("chargeprod:", chargeprod)
+                        #print("sigmaij:", sigmaij)
+                        #print("epsilonij:", epsilonij)
                         chargeprod._value = 0.0
                         force.setExceptionParameters(exc, p1, p2, chargeprod, sigmaij, epsilonij)
-                        # print("New:", force.getExceptionParameters(exc))
+                        #print("New:", force.getExceptionParameters(exc))
         print_time_rel(timeA, modulename="delete_exceptions")
 
-    # # Function to
-    # def zero_nonbondedforce(self, atomlist, zeroCoulomb=True, zeroLJ=True):
-    #     timeA = time.time()
-    #     print("Zero-ing nonbondedforce")
-
-    #     def charge_sigma_epsilon(charge, sigma, epsilon):
-    #         if zeroCoulomb is True:
-    #             newcharge = charge
-    #             newcharge._value = 0.0
-
-    #         else:
-    #             newcharge = charge
-    #         if zeroLJ is True:
-    #             newsigma = sigma
-    #             newsigma._value = 0.0
-    #             newepsilon = epsilon
-    #             newepsilon._value = 0.0
-    #         else:
-    #             newsigma = sigma
-    #             newepsilon = epsilon
-    #         return [newcharge, newsigma, newepsilon]
-
-    #     # Zero all nonbonding interactions for atomlist
-    #     for force in self.system.getForces():
-    #         if isinstance(force, openmm.NonbondedForce):
-    #             # Setting single particle parameters
-    #             for atomindex in atomlist:
-    #                 oldcharge, oldsigma, oldepsilon = force.getParticleParameters(atomindex)
-    #                 newpars = charge_sigma_epsilon(oldcharge, oldsigma, oldepsilon)
-    #                 print(newpars)
-    #                 force.setParticleParameters(atomindex, newpars[0], newpars[1], newpars[2])
-    #             print("force.getNumExceptions() ", force.getNumExceptions())
-    #             print("force.getNumExceptionParameterOffsets() ", force.getNumExceptionParameterOffsets())
-    #             print("force.getNonbondedMethod():", force.getNonbondedMethod())
-    #             print("force.getNumGlobalParameters() ", force.getNumGlobalParameters())
-    #             # Now doing exceptions
-    #             for exc in range(force.getNumExceptions()):
-    #                 print(force.getExceptionParameters(exc))
-    #                 force.getExceptionParameters(exc)
-    #                 p1, p2, chargeprod, sigmaij, epsilonij = force.getExceptionParameters(exc)
-    #                 # chargeprod._value=0.0
-    #                 # sigmaij._value=0.0
-    #                 # epsilonij._value=0.0
-    #                 newpars2 = charge_sigma_epsilon(chargeprod, sigmaij, epsilonij)
-    #                 force.setExceptionParameters(exc, p1, p2, newpars2[0], newpars2[1], newpars2[2])
-    #                 # print("New:", force.getExceptionParameters(exc))
-    #             # force.updateParametersInContext(self.simulation.context)
-    #         elif isinstance(force, openmm.CustomNonbondedForce):
-    #             print("customnonbondedforce not implemented")
-    #             ashexit()
-
     # Updating LJ interactions in OpenMM object. Used to set LJ sites to zero e.g. so that they do not contribute
     # Can be used to get QM-MM LJ interaction energy
     def update_LJ_epsilons(self, atomlist, epsilons):
@@ -3263,9 +3215,9 @@ def write_xmlfile_nonbonded(resnames=None, atomnames_per_res=None, atomtypes_per
     LJforcelines = []
     for resname, atomtypelist, chargelist, sigmalist, epsilonlist in zip(resnames, atomtypes_per_res, charges_per_res,
                                                                          sigmas_per_res, epsilons_per_res):
-        print("atomtypelist:", atomtypelist)
-        print("chargelist.", chargelist)
-        print("sigmalist", sigmalist)
+        #print("atomtypelist:", atomtypelist)
+        #print("chargelist.", chargelist)
+        #print("sigmalist", sigmalist)
         for atype, charge, sigma, epsilon in zip(atomtypelist, chargelist, sigmalist, epsilonlist):
             if charmm == True:
                 #LJ parameters zero here
@@ -3521,22 +3473,28 @@ def __init__(self, fragment=None, theory=None, charge=None, mult=None, timestep=
         #CASE: ONIOMTHeory that might containOpenMMTheory
         elif isinstance(theory, ash.ONIOMTheory):
             print("This is an ONIOMTheory object")
-            print("ONIOMTheory objects are not currently supported")
+            #print("ONIOMTheory objects are not currently supported")
+            self.theory_runtype ="ONIOM"
             #self.QM_MM_object = theory
             self.ONIOM_object = theory
-            self.theory_runtype ="ONIOM"
-
-            for t in theory.theories_N:
-                if isinstance(t,OpenMMTheory):
-                    print("Found OpenMMTheory object inside ONIOMTheory")
-                    self.openmmobject=t
-                    print("Problem: ONIOMTheory containing an OpenMMTheory is currently not supported yet. Complain to developer")
-                    ashexit()
-            #If nothing found then we create:
+            #MMtheory_index = [t.theorytype for t in theory.theories_N].index("MM")
+            #print("MM theory found at index:", MMtheory_index)
+            #self.openmmobject = theory.theories_N[MMtheory_index]
+            #print("self.openmmobject:", self.openmmobject)
+
+            #for t in theory.theories_N:
+            #    if isinstance(t,OpenMMTheory):
+            #        print("Found OpenMMTheory object inside ONIOMTheory")
+            #        self.openmmobject=t
+            #        print("Warnign: ONIOMTheory containing an OpenMMTheory object is currently not officially supported yet. Complain to developer")
+            #        #ashexit()
+
+            #RB NOTE: Creating a new OpenMMTheory object regardless of whether one exists in the ONIOMTheory
             if self.openmmobject is None:
+                print("Creating new OpenMMTheory object to drive simulation")
                 #Creating dummy OpenMMTheory (basic topology, particle masses, no forces except CMMRemoval)
                 self.openmmobject = OpenMMTheory(fragment=fragment, dummysystem=True, platform=platform, printlevel=printlevel,
-                                hydrogenmass=hydrogenmass, constraints=constraints) #NOTE: might add more options here
+                                hydrogenmass=hydrogenmass, constraints=constraints,) #NOTE: might add more options here
             print("Turning on externalforce option.")
             self.openmm_externalforceobject = self.openmmobject.add_custom_external_force()
 
@@ -4620,6 +4578,8 @@ def run(self, simulation_steps=None, simulation_time=None, metadynamics=False, m
 
                 # Run  step to get full system ONIOM gradient.
                 # Updates OpenMM object with ONIOM forces
+                #Note: Unlike QM/MM we don't do any exit_after_customexternalforce_update here because ONIOM object does not update OpenMM object itself
+                # Easier. Drawback that we may have 2 OpenMMTheory objects defined.
                 energy,gradient=self.ONIOM_object.run(current_coords=current_coords, elems=self.fragment.elems, Grad=True, charge=self.charge, mult=self.mult)
                 if self.printlevel >= 2:
                     print("Energy:", energy)

ash/interfaces/interface_crest.py

@@ -259,7 +259,7 @@ def get_crest_conformers(crest_calcdir='crest-calc',conf_file="crest_conformers.
 
     #Getting energies from title lines
     for i in all_titles:
-        en=float(i)
+        en=float(i[0])
         list_xtb_energies.append(en)
 
     for (els,cs,eny) in zip(all_elems,all_coords,list_xtb_energies):

ash/interfaces/interface_fairchem.py

@@ -1,5 +1,4 @@
 import time
-import numpy as np
 import os
 from ash.modules.module_coords import elemstonuccharges
 from ash.functions.functions_general import ashexit, BC,print_time_rel
@@ -8,20 +7,22 @@
 
 # Simple interface to Fairchem
 
-# Use: 
+# Use:
 
 # VIA MODEL NAMES
 # Models available in version 2: uma-s-1p1 (faster,very good), uma-m-1p1 (slower,best)
-# Example: model_name = "uma-s-1p1"  
+# Example: model_name = "uma-s-1p1"
 # Requires hugging-face token activated in shell
 # e.g. export HF_TOKEN=xxxxxxx
 
 #Via model_file
 #model_file="uma-s-1p1.pt"
 
 class FairchemTheory():
-    def __init__(self, model_name=None, model_file=None, task_name=None, device="cuda", seed=41, numcores=1):
+    def __init__(self, model_name=None, model_file=None, task_name=None, device="cuda", seed=41, numcores=1,
+                 printlevel=2):
 
+        module_init_time=time.time()
         # Early exits
         try:
             import fairchem
@@ -43,23 +44,48 @@ def __init__(self, model_name=None, model_file=None, task_name=None, device="cud
 
         print_line_with_mainheader(f"{self.theorynamelabel}Theory initialization")
 
-
+        self.printlevel=printlevel
         self.task_name=task_name
         self.device=device
         self.model_name=model_name
         self.model_file=model_file
         self.seed=seed
         self.numcores=numcores
 
+        # Counter for runcalls
+        self.runcalls=0
+
         if self.device.lower() == 'cpu':
             #Works ??
             os.environ['OMP_NUM_THREADS'] = str(numcores)
 
+        from fairchem.core import pretrained_mlip, FAIRChemCalculator
+        if self.model_name is not None:
+            print("Model set:", self.model_name)
+            predictor = pretrained_mlip.get_predict_unit(self.model_name, device=self.device)
+            self.calc = FAIRChemCalculator(predictor, task_name=self.task_name, seed=self.seed)
+        elif self.model_file is not None:
+            print("Model-file set:", self.model_file)
+            # TODO: can we fix
+            print("Warning: single-atom systems do not work with this approach")
+            self.calc = FAIRChemCalculator.from_model_checkpoint(self.model_file,
+                                                            task_name=self.task_name, device=self.device,
+                                                            seed=self.seed)
+        else:
+            print("Error:Neither model or model_file was set")
+            ashexit()
+        print_time_rel(module_init_time, modulename=f'{self.theorynamelabel} init', moduleindex=2)
+
+    def cleanup(self):
+        pass
+
     def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None, mm_elems=None,
             elems=None, Grad=False, PC=False, numcores=None, restart=False, label=None, Hessian=False,
             charge=None, mult=None):
 
         module_init_time=time.time()
+        if self.printlevel >= 2:
+            print(BC.OKBLUE,BC.BOLD, f"------------RUNNING {self.theorynamelabel} INTERFACE-------------", BC.END)
         # What elemlist to use. If qm_elems provided then QM/MM job, otherwise use elems list
         if qm_elems is None:
             if elems is None:
@@ -68,39 +94,40 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el
             else:
                 qm_elems = elems
 
-        from fairchem.core import pretrained_mlip, FAIRChemCalculator
-
-        if self.model_name is not None:
-            print("Model set:", self.model_name)
-            predictor = pretrained_mlip.get_predict_unit(self.model_name, device=self.device)
-            calc = FAIRChemCalculator(predictor, task_name=self.task_name, seed=self.seed)
-        elif self.model_file is not None:
-            print("Model-file set:", self.model_file)
-            calc = FAIRChemCalculator.from_model_checkpoint(self.model_file, 
-                                                            task_name=self.task_name, device=self.device,
-                                                            seed=self.seed)
+        if self.runcalls == 0:
+            print("First runcall. Creating atoms object")
+            import ase
+            self.atoms = ase.atoms.Atoms(qm_elems,positions=current_coords)
+            self.atoms.info["charge"] = charge
+            self.atoms.info["spin"] = mult
+            # Assigning calculator
+            self.atoms.calc =self.calc
+        elif len(self.atoms.numbers) != len(current_coords):
+            print("Number-of-atoms mismatch (new molecule?). Creating new atoms object")
+            import ase
+            self.atoms = ase.atoms.Atoms(qm_elems,positions=current_coords)
+            self.atoms.info["charge"] = charge
+            self.atoms.info["spin"] = mult
+            # Assigning calculator
+            self.atoms.calc =self.calc
         else:
-            print("Error:Neither model or model_file was set")
-            ashexit()
-
-        import ase
-        atoms = ase.atoms.Atoms(qm_elems,positions=current_coords)
-        # Setting charge/mult
-        atoms.info["charge"] = charge
-        atoms.info["spin"] = mult
-
-        # Assigning calculator
-        atoms.calc = calc
+            print("Updating coordinates in atoms object")
+            self.atoms.set_positions(current_coords)
 
         # Energy
-        en = atoms.get_potential_energy()
+        en = self.atoms.get_potential_energy()
         self.energy = float(en*ash.constants.evtohar)
-
+        if self.printlevel >= 2:
+            print(f"Single-point {self.theorynamelabel} energy:", self.energy)
         if Grad:
-            forces = atoms.get_forces()
+            forces = self.atoms.get_forces()
             self.gradient = forces/-51.422067090480645
 
+        self.runcalls+=1
+        if self.printlevel >= 2:
+            print(BC.OKBLUE,BC.BOLD,f"------------ENDING {self.theorynamelabel}-INTERFACE-------------", BC.END)
         print_time_rel(module_init_time, modulename=f'{self.theorynamelabel} run', moduleindex=2)
+
         if Grad:
             return self.energy, self.gradient
         else: