BUGFIX: specify integrator type for custom integrators

This way if the user defined integrator is "brownian", the Simulation class
will not check if the kinetic energy exceeds the threshold. Similarly, user could
define a custom langevin integrator and want the same reporter messages.

examples/customIntegrators/activeBD.py

@@ -1,10 +1,10 @@
-"""
+r"""
 Polymer simulations with ActiveBrownianIntegrator
 -------------------------------------------------
 
 This is a sample python script to run a polychrom simulation with the `ActiveBrownianIntegrator' custom integrator in polychrom.contrib.integrators. This integrator is used to simulate a polymer where each mononmer has a different effective temperature and thus a different diffusion coefficient :math:`D_i = k_B T_i / \xi`. Here, we consider an example where there are just two types of monomers, active (A) and inactive (B), where :math:`D_A > D_B` and the user chooses the ratio :math:`D_A / D_B`.
-
 """
+
 import time
 import numpy as np
 import os, sys
@@ -22,7 +22,7 @@
 #1 is A, 0 is B
 ids[int(N/2):] = 0
 
-def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, ntimesteps=200000, blocksize=100):
+def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, ntimesteps=10, blocksize=100):
     """ Run a single simulation on a GPU of a hetero-polymer with A monomers and B monomers. A monomers
     have a larger diffusion coefficient than B monomers, with an activity ratio of D_A / D_B.
     
@@ -39,7 +39,7 @@ def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, ntimestep
     timestep : int
         timestep to feed the Brownian integrator (in femtoseconds)
     ntimesteps : int
-        number of timesteps to run the simulation for
+        number of blocks to run the simulation for. For a chain of 1000 monomers, need ~100000 blocks of 100 timesteps to equilibrate.
     blocksize : int
         number of time steps in a block
     
@@ -70,7 +70,7 @@ def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, ntimestep
     reporter = HDF5Reporter(folder="active_inactive", max_data_length=100, overwrite=True)
     sim = simulation.Simulation(
         platform="CUDA", 
-        integrator=integrator,
+        integrator=(integrator, "brownian"),
         timestep=timestep,
         temperature=temperature,
         GPU=gpuid,
@@ -116,4 +116,4 @@ def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, ntimestep
 if __name__ == '__main__':
     gpuid = int(sys.argv[1])
     activity_ratio = int(sys.argv[1])
-    run_monomer_diffusion(gpuid, N, ids, activity_ratio)
+    run_monomer_diffusion(gpuid, N, ids, activity_ratio)

polychrom/simulation.py

@@ -174,8 +174,8 @@ def __init__(self, **kwargs):
             Machines with 1 GPU automatically select their GPU.
 
         integrator : "langevin", "variableLangevin", "verlet", "variableVerlet",
-                     "brownian", optional Integrator to use
-                     (see Openmm class reference)
+                     "brownian", or tuple containing Integrator from Openmm class reference and string defining integrator type.
+                      For user-defined integrators, specify type "brownian" to avoid checking if kinetic energy exceeds max_Ek.
                      
         mass : number or np.array
             Particle mass (default 100 amu)
@@ -315,8 +315,9 @@ def __init__(self, **kwargs):
                 )
         else:
             logging.info("Using the provided integrator object")
-            self.integrator = self.integrator_type
-            self.integrator_type = "UserDefined"
+            integrator, integrator_type = self.integrator_type
+            self.integrator = integrator
+            self.integrator_type = integrator_type
             kwargs["integrator"] = "user_defined"
 
         self.N = kwargs["N"]
@@ -766,7 +767,7 @@ def do_block(
 
         if np.isnan(newcoords).any():
             raise IntegrationFailError("Coordinates are NANs")
-        if eK > self.eK_critical and self.integrator_type.lower() != "brownian":
+        if eK > self.eK_critical and self.integrator_type.lower() != "brownian": 
             raise EKExceedsError("Ek={1} exceeds {0}".format(self.eK_critical, eK))
         if (np.isnan(eK)) or (np.isnan(eP)):
             raise IntegrationFailError("Energy is NAN)")