Test Langevin Integrator

.github/labeler.yml

@@ -1,11 +1,15 @@
-asfe:
-- transformato/annihilation.py
-restraints:
-- transformato/restraints.py
 charmm:
 - transformato/charmm_factory.py
 analysis:
 - transformato/analysis.py
 endstate_correction:
 - transformato/bin/perform_correction.py
 - transformato/bin/slurm_switching.sh
+openmm files:
+- transformato/bin/*run*.py
+test update:
+- transformato/transformato/test*
+infrastructure update:
+- .github/*
+CI update:
+- devtools/*

.github/workflows/CI.yaml

@@ -7,9 +7,9 @@ on:
   push:
     branches:
       - '**'
-  pull_request:
-    branches:
-      - '**'
+#  pull_request:
+#    branches:
+#      - '**'
   schedule:
     # Nightly tests run on master by default:
     #   Scheduled workflows run on the latest commit on the default or base branch.

transformato/bin/drude_openmm_run.py

@@ -67,9 +67,9 @@
     system = top.createSystem(params, **nboptions)
 
 
-if inputs.vdw == "Force-switch" and fftype != "amber" and env != "vacuum":
-    print(f"Setting the vdw switching function to: Force-switch")
-    system = vfswitch(system, top, inputs)
+# if inputs.vdw == "Force-switch" and fftype != "amber" and env != "vacuum":
+#     print(f"Setting the vdw switching function to: Force-switch")
+#     system = vfswitch(system, top, inputs)
 if hasattr(inputs, "lj_lrc") and inputs.lj_lrc == "yes" and env != "vacuum":
     print(f"We will use LJ Long range correction (LRC)")
     for force in system.getForces():
@@ -82,27 +82,28 @@
             force.setUseLongRangeCorrection(True)
 
 if env != "vacuum":
-    barostat = MonteCarloBarostat(inputs.p_ref * bar, inputs.temp * kelvin)
+    barostat = MonteCarloBarostat(inputs.p_ref * unit.bar, inputs.temp * unit.kelvin)
     system.addForce(barostat)
 
 # integrator = LangevinIntegrator(
 #     inputs.temp * kelvin, 1 / unit.picosecond, inputs.dt * unit.picoseconds
 # )
 
-integrator = DrudeNoseHooverIntegrator(
+integrator = DrudeLangevinIntegrator(
     inputs.temp * kelvin,
-    10 / picosecond,
-    1 * kelvin,
-    200 / picosecond,
-    0.0005 * picoseconds,
+    10 / unit.picosecond,
+    1 * unit.kelvin,
+    200 / unit.picosecond,
+    0.0001 * unit.picoseconds,
 )
 
-# integrator.setMaxDrudeDistance(0.2 * angstroms)
-# if integrator.getMaxDrudeDistance() == 0:
-#     print("No Drude Hard Wall Contraint in use")
-# else:
-#     print("Drude Hard Wall set to {}".format(integrator.getMaxDrudeDistance()))
+integrator.setMaxDrudeDistance(0.2 * unit.angstroms)
+if integrator.getMaxDrudeDistance() == 0:
+    print("No Drude Hard Wall Contraint in use")
+else:
+    print("Drude Hard Wall set to {}".format(integrator.getMaxDrudeDistance()))
 
+# TEST: HELLO 2!
 
 # Set platform
 platform = Platform.getPlatformByName("CUDA")
@@ -151,10 +152,13 @@
 print("\nInitial system energy")
 print(simulation.context.getState(getEnergy=True).getPotentialEnergy())
 
+# Drude VirtualSites
+simulation.context.computeVirtualSites()
+
 # Energy minimization
 if inputs.mini_nstep > 0:
-    print("\nEnergy minimization: %s steps" % inputs.mini_nstep)
-    simulation.minimizeEnergy(maxIterations=inputs.mini_nstep)
+    print("\nEnergy minimization:")
+    simulation.minimizeEnergy()
     print(simulation.context.getState(getEnergy=True).getPotentialEnergy())
 
 # Generate initial velocities
@@ -164,6 +168,29 @@
         simulation.context.setVelocitiesToTemperature(inputs.gen_temp, inputs.gen_seed)
     else:
         simulation.context.setVelocitiesToTemperature(inputs.gen_temp)
+
+## Do some additional pre-equilibration when using Drude particles
+print("Doing a first equilibration run")
+simulation.step(100_000)
+
+print("Doing a second equilibration run")
+simulation.integrator.setStepSize(0.0002 * unit.picoseconds)
+simulation.context.reinitialize(preserveState=True)
+simulation.step(100_000)
+
+print("Doing a third equilibration run")
+simulation.integrator.setStepSize(0.0003 * unit.picoseconds)
+simulation.context.reinitialize(preserveState=True)
+simulation.step(100_000)
+
+print("Doing a fourth equilibration run")
+simulation.integrator.setStepSize(0.0004 * unit.picoseconds)
+simulation.context.reinitialize(preserveState=True)
+simulation.step(100_000)
+
+print("Starting the actual simulation")
+simulation.integrator.setStepSize(inputs.dt * unit.picoseconds)
+simulation.context.reinitialize(preserveState=True)
 
 # Production
 print("\nMD run: %s steps" % inputs.nstep)

transformato/mutate.py

@@ -378,6 +378,11 @@ def __init__(
             self.s1_tlc = s1.tlc
             self.asfe: bool = True
             self.dummy_region_cc1: DummyRegion
+            self.drude: bool = False
+            for atom in self.psf1["waterbox"].view[f":{s1.tlc}"].atoms:
+                if atom.type.startswith("DR"):
+                    self.drude = True
+                    logger.info("Assuming Drude particles are present")
 
     def _check_cgenff_versions(self):
         cgenff_sys1 = self.system["system1"].cgenff_version
@@ -819,29 +824,46 @@ def finish_common_core(
             # all atoms should become dummy atoms in the end
             central_atoms = nx.center(self.graphs["m1"])
 
-            # Assure, that the central atom is no hydrogen
-            for atom in self.psf1["waterbox"][f":{self.s1_tlc}"].atoms:
-                if atom.idx in central_atoms:
-                    if atom.name.startswith("H") == True:
-                        raise RuntimeError(
-                            f"One of the central atoms seems to be a hydrogen atom"
+            if not self.drude:
+
+                # Assure, that the central atom is no hydrogen
+                for atom in self.psf1["waterbox"][f":{self.s1_tlc}"].atoms:
+                    if atom.idx in central_atoms:
+                        if atom.name.startswith("H") == True:
+                            raise RuntimeError(
+                                f"One of the central atoms seems to be a hydrogen atom"
+                            )
+
+                # calculate the ordering or LJ mutations
+                if not odered_connected_dummy_regions_cc1:
+                    odered_connected_dummy_regions_cc1 = (
+                        calculate_order_of_LJ_mutations_asfe(
+                            central_atoms,
+                            self.graphs["m1"].copy(),
                         )
+                    )
 
-            # calculate the ordering or LJ mutations
-            if not odered_connected_dummy_regions_cc1:
-                odered_connected_dummy_regions_cc1 = (
-                    calculate_order_of_LJ_mutations_asfe(
-                        central_atoms,
-                        self.graphs["m1"].copy(),
+                if odered_connected_dummy_regions_cc1:
+                    odered_connected_dummy_regions_cc1 = self._check_for_lp(
+                        odered_connected_dummy_regions_cc1,
+                        self.psf1["waterbox"],
+                        self.s1_tlc,
+                        "m1",
                     )
-                )
+            else:
+                odered_connected_dummy_regions_cc1 = []
 
-            if odered_connected_dummy_regions_cc1:
-                odered_connected_dummy_regions_cc1 = self._check_for_lp(
-                    odered_connected_dummy_regions_cc1,
-                    self.psf1["waterbox"],
-                    self.s1_tlc,
-                    "m1",
+                for atom in self.psf1["vacuum"].view[f":{self.s1_tlc}"]:
+                    if atom.type != "DRUD":
+                        if not atom.type.startswith("L"):
+                            # if not atom.type.startswith("H"):
+                            odered_connected_dummy_regions_cc1.append(atom.idx)
+
+                odered_connected_dummy_regions_cc1 = list(
+                    [odered_connected_dummy_regions_cc1]
+                )
+                logger.info(
+                    f"Using the order as given in the psf file {odered_connected_dummy_regions_cc1}"
                 )
 
             self.dummy_region_cc1 = DummyRegion(

transformato/system.py

@@ -500,7 +500,7 @@ def _generate_rdkit_mol(
                         "PSF to mol conversion did not work! Remove the TLC.sdf file in the tlc folder and try again!"
                     )
 
-            return mol
+        return mol
 
     def generate_atom_tables_from_psf(
         self, psf: pm.charmm.CharmmPsfFile, mol: Chem.rdchem.Mol

transformato/tests/test_drude.py

@@ -19,7 +19,7 @@ def test_create_asfe_system():
     configuration = load_config_yaml(
         config=f"/site/raid5/marvin/drude_test/data/config/acetamide.yaml",
         input_dir=f"/site/raid5/marvin/drude_test/data/",
-        output_dir=f"/site/raid5/marvin/drude_test",
+        output_dir=f"/site/raid5/marvin/drude_test/test_1900/",
     )
 
     s1 = SystemStructure(configuration, "structure1")