Add config file

protex/forcefield/config_test.toml

@@ -0,0 +1,36 @@
+name = "Config for IM1H OAC System"
+
+[forcefield]
+
+
+[[states]] # array of tables
+donor = "IM1H"
+donor_name = "H7"
+acceptor = "IM1"
+acceptor_name = "N2"
+
+[[states]]
+donor = "HOAC"
+donor_name = "H"
+acceptor = "OAC"
+acceptor_name = "O2"
+
+[[updates]]
+reaction = ["IM1H", "OAC"]
+r_max = 0.155
+prob = 0.994
+
+[[updates]]
+reaction = ["IM1H", "OAC"]
+r_max = 0.155
+prob = 0.994
+
+[[updates]]
+reaction = ["IM1H", "OAC"]
+r_max = 0.155
+prob = 0.994
+
+[[updates]]
+reaction = ["IM1H", "OAC"]
+r_max = 0.155
+prob = 0.994

protex/system.py

@@ -9,6 +9,11 @@
 import parmed
 import yaml
 
+try:
+    import tomllib
+except ModuleNotFoundError:
+    import toml as tomllib
+
 try:
     import openmm
 except ImportError:
@@ -78,6 +83,35 @@ def __init__(
             [value["r_max"] for value in self.allowed_updates.values()]
         )
 
+    @classmethod
+    def from_config(cls, file: str):
+        config = tomllib.load(file)
+        try:
+            config["states"]
+            config["updates"]
+        except KeyError as e:
+            e.args = (
+                *e.args,
+                "the config file must contain arrays named 'states' and 'updates'",
+            )
+            raise
+
+        state_list = []
+        for state in config["states"]:
+            s = {}
+            s[state["donor"]] = {"donor_name": state["donor_name"]}
+            s[state["acceptor"]] = {"acceptor_name": state["acceptor_name"]}
+            state_list.append(s)
+
+        allowed_updates = {}
+        for update in config["updates"]:
+            allowed_updates[frozenset(update["reaction"])] = {
+                "r_max": update["r_max"],
+                "prob": update["prob"],
+            }
+
+        return cls(state_list, allowed_updates)
+
     def get_update_value_for(self, residue_set: frozenset[str], property: str) -> float:
         """
         returns the value in the allowed updates dictionary

protex/tests/test_simulation.py

@@ -13,6 +13,9 @@
 except ImportError:
     from simtk.openmm.app import DCDReporter, StateDataReporter
 
+# hier ist ein kommentar
+def add(a,b):
+    return 5
 
 @pytest.mark.skipif(
     os.getenv("CI") == "true",

protex/tests/test_system.py

@@ -7,6 +7,10 @@
 from collections import defaultdict
 from sys import stdout
 
+import pytest
+
+import protex
+
 try:  # Syntax changed in OpenMM 7.6
     import openmm as mm
     from openmm import (
@@ -51,13 +55,9 @@
     from simtk.openmm.app import Simulation
     from simtk.unit import angstroms, kelvin, picoseconds
 
-import pytest
-
-import protex
-
 from ..reporter import ChargeReporter, EnergyReporter
 from ..system import IonicLiquidSystem, IonicLiquidTemplates
-from ..testsystems import (
+from ..testsystems import (  # generate_im1h_oac_dummy_system,; generate_im1h_oac_system_clap,
     IM1H_IM1,
     OAC_HOAC,
     generate_im1h_oac_system,
@@ -177,6 +177,11 @@ def test_setup_simulation():
     assert nr_of_particles == 17500 + 500  # +lps for im1 im1h
 
 
+def test_load_config():
+    templates = IonicLiquidTemplates.from_config("protex/forcefield/config_test.toml")
+    print(templates)
+
+
 def test_run_simulation():
 
     simulation = generate_im1h_oac_system()
@@ -615,6 +620,7 @@ def test_drude_forces():
     atom_idxs = defaultdict(list)  # store atom_idxs
     atom_names = defaultdict(list)  # store atom_names
     names = []  # store names
+    pair_12_13_list = ionic_liquid._build_exclusion_list(ionic_liquid.topology)
 
     # iterate over residues, select the first residue for HOAC and OAC and save the individual bonded forces
     for ridx, r in enumerate(topology.residues()):
@@ -636,12 +642,11 @@ def test_drude_forces():
                 if idx1 in atom_idxs[r.name] and idx2 in atom_idxs[r.name]:
                     print(f)
                     force_state[r.name].append(f)
-
             print("thole")
             print(drude_force.getNumScreenedPairs())
             for drude_id in range(drude_force.getNumScreenedPairs()):
                 f = drude_force.getScreenedPairParameters(drude_id)
-                parent1, parent2 = ionic_liquid.pair_12_13_list[drude_id]
+                parent1, parent2 = pair_12_13_list[drude_id]
                 drude1, drude2 = parent1 + 1, parent2 + 1
                 # print(f"thole {idx1=}, {idx2=}")
                 # print(f"{drude_id=}, {f=}")
@@ -672,7 +677,7 @@ def test_drude_forces():
             print(drude_force.getNumScreenedPairs())
             for drude_id in range(drude_force.getNumScreenedPairs()):
                 f = drude_force.getScreenedPairParameters(drude_id)
-                parent1, parent2 = ionic_liquid.pair_12_13_list[drude_id]
+                parent1, parent2 = pair_12_13_list[drude_id]
                 drude1, drude2 = parent1 + 1, parent2 + 1
                 # print(f"thole {idx1=}, {idx2=}")
                 # print(f"{drude_id=}, {f=}")
@@ -968,6 +973,67 @@ def test_write_psf_save_load():
     os.remove("checkpoint.rst")
 
 
+@pytest.mark.skipif(
+    os.getenv("CI") == "true",
+    reason="Files not on remote",
+)
+def test_write_psf_save_load_clap():
+    psf = "/site/raid3/florian/clap/b3lyp/im1h_oac_150_im1_hoac_350.psf"
+    crd = "/site/raid3/florian/clap/b3lyp/im1h_oac_150_im1_hoac_350.crd"
+    PARA_FILES = [
+        "polarizable_flo_dummy.rtf",
+        "polarizable_flo_dummy.prm",
+    ]
+    para_files = [
+        f"/site/raid3/florian/clap/toppar/{para_file}" for para_file in PARA_FILES
+    ]
+
+    simulation = generate_im1h_oac_system_clap(
+        psf_file=psf, crd_file=crd, para_files=para_files, dummy_atom_type="DUM"
+    )
+    # get ionic liquid templates
+    allowed_updates = {}
+    allowed_updates[frozenset(["IM1H", "OAC"])] = {"r_max": 0.145, "prob": 1}
+    allowed_updates[frozenset(["IM1", "HOAC"])] = {"r_max": 0.145, "prob": 1}
+
+    IM1H_IM1 = {"IM1H": {"atom_name": "HN1"}, "IM1": {"atom_name": "NA1"}}
+    OAC_HOAC = {"OAC": {"atom_name": "O1"}, "HOAC": {"atom_name": "H3"}}
+
+    templates = IonicLiquidTemplates([OAC_HOAC, IM1H_IM1], (allowed_updates))
+    # wrap system in IonicLiquidSystem
+    ionic_liquid = IonicLiquidSystem(simulation, templates)
+    # initialize update method
+    update = NaiveMCUpdate(ionic_liquid)
+    # initialize state update class
+    state_update = StateUpdate(update)
+
+    old_psf_file = psf
+    ionic_liquid.write_psf(old_psf_file, "test1.psf")
+
+    # ionic_liquid.simulation.step(50)
+    state_update.update(2)
+
+    ionic_liquid.write_psf(old_psf_file, "test2.psf")
+
+    ionic_liquid.simulation.step(10)
+    state_update.update(2)
+
+    ionic_liquid.write_psf(old_psf_file, "test3.psf")
+
+    ionic_liquid.saveState("state.rst")
+    ionic_liquid.saveCheckpoint("checkpoint.rst")
+
+    ionic_liquid2 = ionic_liquid  # copy.deepcopy(ionic_liquid)
+    ionic_liquid.loadState("state.rst")
+    ionic_liquid2.loadCheckpoint("checkpoint.rst")
+
+    # os.remove("test1.psf")
+    # os.remove("test2.psf")
+    # os.remove("test3.psf")
+    # os.remove("state.rst")
+    # os.remove("checkpoint.rst")
+
+
 @pytest.mark.skipif(
     os.getenv("CI") == "true",
     reason="Will fail sporadicaly.",
@@ -1030,6 +1096,96 @@ def print_force_contrib(simulation):
     # )
 
 
+@pytest.mark.skipif(
+    os.getenv("CI") == "true",
+    reason="Will fail sporadicaly.",
+)
+def test_dummy():
+    def get_time_energy(simulation, print=False):
+        time = simulation.context.getState().getTime()
+        e_pot = simulation.context.getState(getEnergy=True).getPotentialEnergy()
+        if print:
+            print(f"time: {time}, e_pot: {e_pot}")
+        return time, e_pot
+
+    def save_il(ionic_liquid, number):
+        ionic_liquid.write_psf(
+            f"protex/forcefield/dummy/im1h_oac_im1_hoac_1.psf",
+            f"test_{number}.psf",
+        )
+        ionic_liquid.saveCheckpoint(f"test_{number}.rst")
+
+    def load_sim(psf, rst):
+        sim = generate_single_im1h_oac_system(psf_file=psf)
+        sim.loadCheckpoint(rst)
+        return sim
+
+    def load_il(psf, rst, templates):
+        sim = generate_single_im1h_oac_system(psf_file=psf)
+        il = IonicLiquidSystem(sim, templates)
+        il.loadCheckpoint(rst)
+        return il
+
+    def print_force_contrib(simulation):
+        for i, f in enumerate(simulation.system.getForces()):
+            group = f.getForceGroup()
+            state = simulation.context.getState(getEnergy=True, groups={group})
+            print(f.getName(), state.getPotentialEnergy())
+
+    simulation = generate_im1h_oac_dummy_system()
+    nonbonded_force = [
+        f for f in simulation.system.getForces() if isinstance(f, mm.NonbondedForce)
+    ][0]
+    dummy_atoms = []
+    for atom in simulation.topology.atoms():
+        if atom.residue.name == "IM1" and atom.name == "H7":
+            dummy_atoms.append(atom.index)
+            print(atom)
+            print(nonbonded_force.getParticleParameters(atom.index))
+            # nonbonded_force.setParticleParameters(atom.index, 0.0, 0.0, 0.0)
+            # print(nonbonded_force.getParticleParameters(atom.index))
+        if atom.residue.name == "OAC" and atom.name == "H":
+            dummy_atoms.append(atom.index)
+            print(nonbonded_force.getParticleParameters(atom.index))
+            # nonbonded_force.setParticleParameters(atom.index, 0.0, 0.0, 0.0)
+    for exc_id in range(nonbonded_force.getNumExceptions()):
+        f = nonbonded_force.getExceptionParameters(exc_id)
+        idx1 = f[0]
+        idx2 = f[1]
+        chargeProd, sigma, epsilon = f[2:]
+        if idx1 in dummy_atoms or idx2 in dummy_atoms:
+            # nonbonded_force.setExceptionParameters(exc_id, idx1, idx2, 0.0, sigma, 0.0)
+            print(f)
+
+    # get ionic liquid templates
+    allowed_updates = {}
+    allowed_updates[frozenset(["IM1H", "OAC"])] = {"r_max": 0.785, "prob": 1}
+    # allowed_updates[frozenset(["IM1", "HOAC"])] = {"r_max": 0.165, "prob": 1}
+
+    templates = IonicLiquidTemplates([OAC_HOAC, IM1H_IM1], (allowed_updates))
+    # wrap system in IonicLiquidSystem
+    ionic_liquid = IonicLiquidSystem(simulation, templates)
+    # initialize update method
+    update = NaiveMCUpdate(ionic_liquid)
+    # initialize state update class
+    state_update = StateUpdate(update)
+
+    save_il(ionic_liquid, 0)
+    # for i in range(5):
+    #     state_update.update(2)
+    #     ionic_liquid.simulation.step(10)
+    #     print(ionic_liquid.simulation.context.getPlatform().getName())
+
+    state_update.update(2)
+
+    save_il(ionic_liquid, 1)
+
+    # sim2_1 = load_sim("protex/forcefield/single_pairs/im1_hoac_2.psf", "test_2.rst")
+    # sim_2_oldcoord = load_sim(
+    #    "protex/forcefield/single_pairs/im1_hoac_2.psf", "test_1.rst"
+    # )
+
+
 def test_single_harmonic_force(caplog):
     caplog.set_level(logging.DEBUG)