MTSIntegrator

[gitkol]

Just a clarification, the PLUMED plugin injects PLUMED forces at every time step, correct? Does this mean the outer time step when using the MTSIntegrator? Thanks, Istvan

[peastman]

PlumedForce works just like any other force. You can assign it to any force group you want (call setForceGroup() on it). MTSIntegrator will then evaluate it at whatever point you tell it to use that force group.

[gitkol]

Thanks, Peter. That makes it crystal clear,

[gitkol]

Sorry, Peter. Just a follow-up question. I am running a series of simulations within a single master simulation. The different sub-simulations differ only in using different plumed forces, but the system remains the same. I use the following workflow to switch from one plumed force to another:

# Subtract plumed force and replace it
system.removeForce(plumed1_fi)
plumed2_fi = system.addForce(PlumedForce(plumed2))
# Reinitialize simulation
simulation.context.reinitialize(preserveState=True)
simulation.step(steps)

My question regarding MTSIntegrator is whether I have to redefine it every time I switch forces. Initially, MTS is defined as follows:

# Integrator
nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
nonbonded.setReciprocalSpaceForceGroup(1)
plumed = [f for f in system.getForces() if isinstance(f, Force)][0]
plumed.setForceGroup(2)
integrator = MTSIntegrator(dt, [(2,1), (1,1), (0,3)])
integrator.setConstraintTolerance(constraintTolerance)

When I remove the plumed1 force, group 2 becomes undefined. When I add the plumed2 force group 2 becomes relevant again, but it is not clear to me whether I have to redefine the integrator (and consequently the whole simulation).

Could you clarify the correct workflow?

Thanks a million,

Istvan

[peastman]

That ought to work. A force group is really just a number. When you ask the integrator to compute the forces or energy for a particular group, it goes to every Force object in the system and asks each one to compute any contributions it has for that group. If none of them contributes anything, the force and energy will come out as 0.

When you call reinitialize() on a Context, it throws out all of its internal data structures and rebuilds them from scratch. It's effectively the same as creating a completely new Context with a new (but identical) Integrator.

[gitkol]

That sounds great! I can just define then the integrator once in the header of the simulation script. That makes it so much more elegant. Thanks, Peter!

[gitkol]

Really sorry to have to reopen again, but it looks like there is another complication. Printing system.getForces() results in the following list.

[<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7fca1b11e360> >, <simtk.openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7fca1b11e4e0> >, <simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7fca1b11e660> >, <simtk.openmm.openmm.PeriodicTorsionForce; proxy of <Swig Object of type 'OpenMM::PeriodicTorsionForce *' at 0x7fca1b11e600> >, <simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7fca1b11e6f0> >, <simtk.openmm.openmm.AndersenThermostat; proxy of <Swig Object of type 'OpenMM::AndersenThermostat *' at 0x7fca1b11e690> >, <simtk.openmm.openmm.MonteCarloBarostat; proxy of <Swig Object of type 'OpenMM::MonteCarloBarostat *' at 0x7fca1b11e780> >, <simtk.openmm.openmm.Force; proxy of <Swig Object of type 'OpenMM::Force *' at 0x7fca1b11e540> >]

According to this, the plumed force that I injected just before printing the force list, is the last item in the list called Force. However, if I use the isinstance() filter above using the keyword Force, it'll incorrectly pick the first item that has 'Force' in its name, which happens to be <simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7fca1b11e5a0> > and not the plumed force. Of course, I can do a workaround saying plumed = [f for f in system.getForces()][-1] but that won't be a general solution. It looks like plumed force is injected as the general simtk.openmm.openmm.Force() object. I wonder if Force should be called by some unique name instead.

I apologize if I overlooked something.

Thanks,

Istvan

[peastman]

Force is the superclass of every force type. You want to be checking for a PlumedForce.

[gitkol]

Yes, but PlumedForce does not appear in the system.getForces() list even though I injected it via

drr = """
DISTANCE ATOMS=4,114 LABEL=dist
DRR ARG=dist FULLSAMPLES=500 GRID_MIN=1.7 GRID_MAX=3.8 GRID_SPACING=0.01 TEMP=310 FRICTION=8.0 TAU=0.5 OUTPUTFREQ=10000 HISTORYFREQ=250000 LABEL=drr
PRINT ARG=drr.bias,drr.dist_biasforce,drr.dist_fict,dist STRIDE=1000 FILE=plumed-elf-test.log"""
drr_fi = system.addForce(PlumedForce(drr))

[gitkol]

If I compare the output of system.getForces() before and after injecting the plumed force, the difference is <simtk.openmm.openmm.Force; proxy of <Swig Object of type 'OpenMM::Force *' at 0x7f5858d31570> > added to the list.

[peastman]

Ok, I see. That's an artifact of how the Python wrapping works.

Why do you need to look up the PlumedForce at all? You already know its index in the list, because it was returned by addForce() and you saved it in a variable. You can just write

plumed = system.getForce(plumed2_fi)

[gitkol]

Oh, I see. Then, I can set the group via plumed.setForceGroup(2).

[peastman]

Exactly!

[gitkol]

Yes, this is working. Thanks!