RPMD Contraction Scheme

[Ali Eltareb]

Hi,

I'm trying to implement the ring-polymer contraction scheme (RPC), however I've run into some strange problems. I can't get the speed up of ~6 when using 32 beads per ring-polymer.

I'm using a flexible forcefield model of water (q-TIP4P/F), which contains bonded and nonbonded forces. I first set every force group to its own index, I then set the reciprocal space force group from the nonbonded force to its own index. The below code shows how I was partitioning the force groups.

Code: Select all

k = 0
for force in system.getForces():
    force.setForceGroup(k)
    k =  k + 1

# The second force group is the nonbonded force group that contains the coulomb force and LJ force
system.getForce(2).setReciprocalSpaceForceGroup(5)

integrator = RPMDIntegrator(32, 300*kelvin, 1.0/picosecond, 0.00025*picoseconds, {2:5, 5:1})

This script only speeds up the simulation by ~2 when compared to a uncontracted simulation.

The strange behaviour occurs when I print out the contractions. So If I include a print statement right after the integrator line above, OpenMM gives me a different answer if I was to include a print statement right after simulating for 100 steps. Here is the code I was using and its output

Code: Select all

k = 0
for force in system.getForces():
    print (force,force.setForceGroup(k))
    k =  k + 1

# Check that every force group got its own index
for force in system.getForces():
    print (force,force.getForceGroup())

system.getForce(2).setReciprocalSpaceForceGroup(5)

steps = 100
integrator = RPMDIntegrator(32, 300*kelvin, 1.0/picosecond, 0.00025*picoseconds, {2:3,5:3})
print (integrator.getContractions())

platform = Platform.getPlatformByName('CPU')
simulation = Simulation(pdb.topology, system, integrator, platform)#, properties)
simulation.context.setPositions(pdb.positions)
simulation.context.setVelocitiesToTemperature(T_K)
simulation.step(steps)
print ('!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print (integrator.getContractions())

Code: Select all

<simtk.openmm.openmm.CustomBondForce; proxy of <Swig Object of type 'OpenMM::CustomBondForce *' at 0x7f005dda9ab0> > None
<simtk.openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7f005dda9c00> > None
<simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7f005dda9bd0> > None
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f005dda9b40> > None
<simtk.openmm.openmm.CustomBondForce; proxy of <Swig Object of type 'OpenMM::CustomBondForce *' at 0x7f005dda9c00> > 0
<simtk.openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7f005dda9b10> > 1
<simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7f005dda9ba0> > 2
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f005dda9ae0> > 3
<simtk.openmm.openmm.mapii; proxy of <Swig Object of type 'std::map< int,int > *' at 0x7f005dda9c00> >
!!!!!!!!!!!!!!!!!!!!!!!!!!!
<simtk.openmm.openmm.mapii; proxy of <Swig Object of type 'std::map< int,int > *' at 0x7f005dda9ba0> >

Notice how the contractions are different in both print statements. Why is this occurring? Is there a problem with the code I wrote? Also why is it only showing one mapping when I have specified two force groups to be contracted?

[Peter Eastman]

How much speedup you get from contractions depends on the details of your simulation. The only part you're speeding up is the nonbonded force. Everything else in the simulation (bonded forces, integration, constraints, etc.) still takes just as long. If nonbonded forces only take half the total computation time, then 2x would be the upper limit on how much speedup you could theoretically hope to get from optimizing them.

The relative costs of different parts of the calculation depend on various factors: the size of the system, the Ewald error tolerance, the cutoff distance, what other forces are present, whether you have constraints, etc. You could try running in a profiler to see where the time is being spent.

[Peter Eastman]

Change your print statement to

Code: Select all

print(dict(integrator.getContractions()))

The C++ method returns a map<int, int>. SWIG automatically generates a Python wrapper around it, which behaves exactly like a Python dict. Copying the contents into an actual dict gives you a more useful string representation.