- HDLC default for QM/MM or ActRegion systems, force_coordsystem to b…

…ypass

- some work on ONIOM but not complete
- MTD: user_cvforce and user_biasvar, work on supporting cbpack, reenabled plumed

ash/interfaces/interface_OpenMM.py

@@ -3464,6 +3464,8 @@ def __init__(self, fragment=None, theory=None, charge=None, mult=None, timestep=
         self.force_file_option=force_file_option #Gradients/forces as a file
         self.energy_file_option=energy_file_option #Energies as a file
         self.atomic_units_force_reporter=atomic_units_force_reporter #Forces in atomic units
+        self.user_cvforce=None # Initializing possibility of user CV object
+        self.user_biasvar=None #Initializing possibility of user biasvariable
 
         #PERIODIC or not
         if self.openmmobject.Periodic is True:
@@ -3788,6 +3790,7 @@ def run(self, simulation_steps=None, simulation_time=None, metadynamics=False, m
             print(f"Workerdir: {workerdir} provided. Entering dir")
             os.chdir(workerdir)
 
+
         #If using Plumed then now we add Plumed-force to system from plumedinput string
         if plumedinput != None:
             import openmmplumed
@@ -3819,9 +3822,9 @@ def run(self, simulation_steps=None, simulation_time=None, metadynamics=False, m
             if metadyn_settings["numCVs"] == 2:
                 #Creating CV biasvariables and forces
                 CV1_bias,cvforce_1 = create_CV_bias(metadyn_settings["CV1_type"],metadyn_settings["CV1_atoms"],metadyn_settings["CV1_biaswidth"],
-                                                    CV_range=metadyn_settings["CV1_range"], reference_pos=reference_pos, reference_particles=metadyn_settings["CV1_atoms"])
+                                                    CV_range=metadyn_settings["CV1_range"], reference_pos=reference_pos, reference_particles=metadyn_settings["CV1_atoms"], user_cvforce=self.user_cvforce, user_biasvar=self.user_biasvar)
                 CV2_bias,cvforce_2 = create_CV_bias(metadyn_settings["CV2_type"],metadyn_settings["CV2_atoms"],metadyn_settings["CV2_biaswidth"],
-                                                    CV_range=metadyn_settings["CV2_range"], reference_pos=reference_pos, reference_particles=metadyn_settings["CV2_atoms"])
+                                                    CV_range=metadyn_settings["CV2_range"], reference_pos=reference_pos, reference_particles=metadyn_settings["CV2_atoms"], user_cvforce=self.user_cvforce, user_biasvar=self.user_biasvar)
 
                 #Gridwidth and min/max values now set. Adding to dict
                 metadyn_settings["CV1_gridwidth"] = CV1_bias.gridWidth
@@ -3844,7 +3847,8 @@ def run(self, simulation_steps=None, simulation_time=None, metadynamics=False, m
             elif metadyn_settings["numCVs"] == 1:
                 #Creating CV biasvariable and force
                 CV1_bias,cvforce_1 = create_CV_bias(metadyn_settings["CV1_type"],metadyn_settings["CV1_atoms"],metadyn_settings["CV1_biaswidth"],
-                                                    CV_range=metadyn_settings["CV1_range"], reference_pos=reference_pos, reference_particles=metadyn_settings["CV1_atoms"])
+                                                    CV_range=metadyn_settings["CV1_range"], reference_pos=reference_pos, reference_particles=metadyn_settings["CV1_atoms"],
+                                                    user_cvforce=self.user_cvforce, user_biasvar=self.user_biasvar)
                 #Gridwidth and min/max values now set. Adding to dict
                 metadyn_settings["CV1_gridwidth"] = CV1_bias.gridWidth
                 metadyn_settings["CV1_minvalue"] = CV1_bias.minValue
@@ -4145,12 +4149,13 @@ def run(self, simulation_steps=None, simulation_time=None, metadynamics=False, m
 
 
                 #OpenMM metadynamics
-                if metadynamics == True:
+                if metadynamics is True:
                     if self.printlevel >= 2:
                         print("Now calling OpenMM native metadynamics and taking 1 step")
                     meta_object.step(self.simulation, 1)
 
                     #getCollectiveVariables
+                    cv1scaling=1
                     if step % metadyn_settings["saveFrequency"]*metadyn_settings["frequency"] == 0:
                         if self.printlevel >= 2:
                             print("MTD: Writing current collective variables to disk")
@@ -4616,6 +4621,7 @@ def OpenMM_metadynamics(fragment=None, theory=None, timestep=0.001, simulation_s
               CV1_atoms=None, CV2_atoms=None, CV1_type=None, CV2_type=None, biasfactor=6,
               height=1,
               CV1_biaswidth=0.5, CV2_biaswidth=0.5, CV1_range=None, CV2_range=None,
+              user_cvforce=None, user_biasvar=None,
               frequency=1, savefrequency=10, printlevel=2,
               biasdir='.', multiplewalkers=False, numcores=1, walkerid=None):
     print_line_with_mainheader("OpenMM metadynamics")
@@ -4667,6 +4673,14 @@ def OpenMM_metadynamics(fragment=None, theory=None, timestep=0.001, simulation_s
                         centerforce_distance=centerforce_distance, centerforce_center=centerforce_center,
                         barostat_frequency=barostat_frequency, specialbox=specialbox, printlevel=printlevel)
 
+    #
+    if user_cvforce is not None:
+        print("User CV-force was given:", user_cvforce)
+        md.user_cvforce=user_cvforce
+    if user_biasvar is not None:
+        print("User Biasvar was given:", user_biasvar)
+        md.user_biasvar=user_biasvar
+
     #Load OpenMM.app
     import openmm
     import openmm.app as openmm_app
@@ -4768,7 +4782,7 @@ def OpenMM_metadynamics(fragment=None, theory=None, timestep=0.001, simulation_s
                                                         restraints=restraints)
     else:
         simulation = md.run(simulation_steps=simulation_steps, simulation_time=simulation_time, metadynamics=native_MTD, metadyn_settings=metadyn_settings,
-                            restraints=restraints)
+                            restraints=restraints, plumedinput=plumedinput)
     print("Metadynamics simulation done")
 
     #Finalizing simulation (writes and updates files)
@@ -4873,7 +4887,7 @@ def Gentle_warm_up_MD(theory=None, fragment=None, time_steps=[0.0005,0.001,0.004
     return
 
 #Function to create CV biases in native OpenMM metadynamics
-def create_CV_bias(CV_type,CV_atoms,biaswidth_cv,CV_range=None, reference_pos=None, reference_particles=None):
+def create_CV_bias(CV_type,CV_atoms,biaswidth_cv,CV_range=None, reference_pos=None, reference_particles=None, user_cvforce=None, user_biasvar=None):
     import openmm
     print("Inside create_CV_bias")
     print("CV_type:", CV_type)
@@ -4912,6 +4926,20 @@ def create_CV_bias(CV_type,CV_atoms,biaswidth_cv,CV_range=None, reference_pos=No
             CV_unit_label="Å"
             biaswidth_cv_unit=openmm.unit.angstroms
             biaswidth_cv_unit_label="Å"
+        elif CV_type.lower() == "cn":
+            CV_min_val=0.0
+            CV_max_val=5.0 #TODO
+            CV_unit=1
+            CV_unit_label="CN"
+            biaswidth_cv_unit=1 # TODO
+            biaswidth_cv_unit_label="CN" #
+        elif CV_type.lower() == "custom":
+            CV_min_val=0.0
+            CV_max_val=5.0 #TODO
+            CV_unit=1
+            CV_unit_label="Custom"
+            biaswidth_cv_unit=1 # TODO
+            biaswidth_cv_unit_label="Custom" #
     else:
         print("CV range given.")
         CV_min_val=CV_range[0]
@@ -4936,6 +4964,16 @@ def create_CV_bias(CV_type,CV_atoms,biaswidth_cv,CV_range=None, reference_pos=No
             CV_unit_label="Å"
             biaswidth_cv_unit=openmm.unit.angstroms
             biaswidth_cv_unit_label="Å"
+        elif CV_type.lower() == "cn":
+            CV_unit=1 #TODO
+            CV_unit_label="CN"
+            biaswidth_cv_unit=1 #TODO
+            biaswidth_cv_unit_label="CN"
+        elif CV_type.lower() == "custom":
+            CV_unit=1 #TODO
+            CV_unit_label="Custom"
+            biaswidth_cv_unit=1 #TODO
+            biaswidth_cv_unit_label="Custom"
     print(f"CV_min_val: {CV_min_val} and CV_max_val: {CV_max_val} {CV_unit_label}")
     print(f"Biaswidth of CV: {biaswidth_cv} {biaswidth_cv_unit_label}")
     # Define collective variables for CV1 and CV2.
@@ -4968,6 +5006,46 @@ def create_CV_bias(CV_type,CV_atoms,biaswidth_cv,CV_range=None, reference_pos=No
         cvforce = openmm.RMSDForce(reference_pos)
         cvforce.setParticles(reference_particles)
         CV_bias = openmm.app.BiasVariable(cvforce, CV_min_val*CV_unit, CV_max_val*CV_unit, biaswidth_cv*biaswidth_cv_unit, periodic=False)
+    elif CV_type.lower() == "cn":
+        print("CV type is CN")
+        try:
+            import cvpack
+        except:
+            print("CV-type CN requires the cvpack package. See https://github.com/RedesignScience/cvpack?tab=readme-ov-file")
+            print("Install like this: mamba install -c mdtools cvpack")
+            ashexit()
+
+        #NH3 group
+        group1=[6]
+        group2=[13,14,15]
+
+        forces = {f.getName(): f for f in self.openmmobject.system.getForces()}
+        nc = cvpack.NumberOfContacts(
+            group1,
+            group2,
+            forces["NonbondedForce"],
+            stepFunction="step(1-x)",
+        )
+
+        cvforce = nc
+
+        #Biasvariable: forceobj, minval, maxval, biaswidth
+        biasvar=openmm.app.BiasVariable(cvforce, 00, 18.0, 2.0, periodic=False)
+
+
+
+
+        ashexit()
+        #cvforce = openmm.CustomTorsionForce('theta')
+        #cvforce.addTorsion(*CV_atoms)
+        #CV_bias = openmm.app.BiasVariable(cvforce, CV_min_val*CV_unit, CV_max_val*CV_unit, biaswidth_cv*biaswidth_cv_unit, periodic=True)
+    elif CV_type.lower() == "custom":
+        #User cv force with atoms already added
+        #cvforce
+        cvforce=user_cvforce
+        CV_bias= user_biasvar
+        print("cvforce:", cvforce)
+        print("CV_bias:", CV_bias)
     else:
         print("unsupported CV_type for native OpenMM metadynamics implementation")
         ashexit()

ash/interfaces/interface_geometric_new.py

@@ -19,7 +19,7 @@
 
 #Wrapper function around GeomeTRICOptimizerClass
 #NOTE: theory and fragment given to Optimizer function but not part of Class initialization. Only passed to run method
-def geomeTRICOptimizer(theory=None, fragment=None, charge=None, mult=None, coordsystem='tric', frozenatoms=None, constraints=None,
+def geomeTRICOptimizer(theory=None, fragment=None, charge=None, mult=None, coordsystem='tric', force_coordsystem=False, frozenatoms=None, constraints=None,
                        constrainvalue=False, maxiter=250, ActiveRegion=False, actatoms=None,
                        convergence_setting=None, conv_criteria=None, print_atoms_list=None, TSOpt=False, hessian=None, partial_hessian_atoms=None,
                        modelhessian=None, subfrctor=1, MM_PDB_traj_write=False, printlevel=2):
@@ -38,7 +38,7 @@ def geomeTRICOptimizer(theory=None, fragment=None, charge=None, mult=None, coord
                         hessian=hessian, partial_hessian_atoms=partial_hessian_atoms,modelhessian=modelhessian,
                         convergence_setting=convergence_setting, conv_criteria=conv_criteria,
                         print_atoms_list=print_atoms_list, subfrctor=subfrctor, MM_PDB_traj_write=MM_PDB_traj_write,
-                        printlevel=printlevel)
+                        printlevel=printlevel, force_coordsystem=force_coordsystem)
 
     #Providing theory and fragment to run method. Also constraints
     result = optimizer.run(theory=theory, fragment=fragment, charge=charge, mult=mult,
@@ -54,7 +54,7 @@ def __init__(self,theory=None, charge=None, mult=None, coordsystem='tric',
                      frozenatoms=None, maxiter=250, ActiveRegion=False, actatoms=None,
                        convergence_setting=None, conv_criteria=None, TSOpt=False, hessian=None,
                        print_atoms_list=None, partial_hessian_atoms=None, modelhessian=None,
-                       subfrctor=1, MM_PDB_traj_write=False, printlevel=2):
+                       subfrctor=1, MM_PDB_traj_write=False, printlevel=2, force_coordsystem=False):
 
             self.printlevel=printlevel
             print_line_with_mainheader("geomeTRICOptimizer initialization")
@@ -72,15 +72,16 @@ def __init__(self,theory=None, charge=None, mult=None, coordsystem='tric',
             if frozenatoms is None:
                 frozenatoms=[]
 
-            if ActiveRegion is True and coordsystem == "tric":
-                print("Warning: ActiveRegion True and coordsystem is TRIC.")
-                print("HDLC coordinate system is usually more robust for large systems.")
-                print("Try coordsystem='hdlc' if you experience problems")
-                #TODO: Look into this more
-                # print("Activeregion true and coordsystem = tric are not compatible")
-                # print("Switching to HDLC")
-                # coordsystem='hdlc'
-                # exit()
+            if ActiveRegion is True and coordsystem.lower() == "tric":
+                print("Warning: ActiveRegion is set but the coordsystem is TRIC. The HDLC coordinate system is usually much more robust for large systems than TRIC.")
+                print("")
+                if force_coordsystem is True:
+                    print("force_coordsystem is True.")
+                    print("Sticking with coordsystem TRIC")
+                else:
+                    print("force_coordsystem is False.")
+                    print("Warning: Now switching to HDLC to avoid likely robustness problems with TRIC. To avoid this behaviour (and force use of TRIC) you can use set the Boolean force_coordsystem to True.")
+                    coordsystem='hdlc'
 
             # Defining some attributes
             self.maxiter=maxiter

ash/modules/module_MM.py

@@ -255,6 +255,16 @@ def update_charges(self,atomlist,charges):
         #print("Charges are now:", charges)
         print("Sum of charges:", sum(charges))
 
+    def update_LJ_epsilons(self,atomlist,epsilons):
+        print("Updating charges.")
+        assert len(atomlist) == len(epsilons)
+        print("Error: update_LJ_epsilons in Nonbondedtheory not ready")
+        ashexit()
+        #for atom,eps in zip(atomlist,epsilons):
+        #    self.atom_charges[atom] = eps
+        #print("Charges are now:", charges)
+        print("Sum of epsilons:", sum(epsilons))
+
     # current_coords is now used for full_coords, charges for full coords
     def run(self, current_coords=None, elems=None, charges=None, connectivity=None, numcores=1, label=None,
             Coulomb=True, Grad=True, qmatoms=None, actatoms=None, frozenatoms=None, charge=None, mult=None):

ash/modules/module_coords.py

@@ -3050,7 +3050,7 @@ def QMPC_fragexpand(theory=None, fragment=None, thresh=5e-4):
 #Now sped up via get_connected_atoms_np. Silly
 def get_boundary_atoms(qmatoms, coords, elems, scale, tol, excludeboundaryatomlist=None, unusualboundary=False):
     timeA = time.time()
-    print("Determining QM-MM boundary")
+    print("Determining QM-MM/HL-LL boundary")
     print("Parameters determing connectivity:")
     print("Scaling factor:", scale)
     print("Tolerance:", tol)