RPMD Contraction Scheme
Posted: Mon Aug 23, 2021 6:41 pm
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.
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
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?
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> >