From bfbaeac5283300c9df687d337886a6e75f856ab8 Mon Sep 17 00:00:00 2001 From: Deepti Date: Wed, 4 Jan 2023 18:49:02 -0500 Subject: [PATCH] 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 | 12 ++++++------ polychrom/simulation.py | 11 ++++++----- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/customIntegrators/activeBD.py b/examples/customIntegrators/activeBD.py index 2c0266d..1d333b7 100755 --- a/examples/customIntegrators/activeBD.py +++ b/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) \ No newline at end of file + run_monomer_diffusion(gpuid, N, ids, activity_ratio) diff --git a/polychrom/simulation.py b/polychrom/simulation.py index a7807e9..6756929 100755 --- a/polychrom/simulation.py +++ b/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)")