Flat-bottom potential

[Marcelo D Poleto]

Hi,

I am trying to apply a flat-bottom potential between a small set of atoms but I could not find an example of how to set that up. I am currently using the code below, but my distances are often higher than the cutoffs (2A and 2.7A).

Am I doing something wrong here?

Code: Select all

# 0-based list
# LIG-O4  = 13944
# LIG-CX2 = 13940
# E226-HN  = 4481
# G132-HN  = 2236
# S225-OG  = 4474

flat_bottom_force = CustomBondForce(
	'step(r-r0) * (k/2) * (r-r0)^2')
flat_bottom_force.addPerBondParameter('r0')
flat_bottom_force.addPerBondParameter('k')
system.addForce(flat_bottom_force)

flat_bottom_force.addBond(13944, 4481, [2.0*angstrom, flat_constant*kilojoules_per_mole/(nanometer**2)])
flat_bottom_force.addBond(13944, 2236, [2.0*angstrom, flat_constant*kilojoules_per_mole/(nanometer**2)])
flat_bottom_force.addBond(13940, 4474, [2.7*angstrom, flat_constant*kilojoules_per_mole/(nanometer**2)])

[Peter Eastman]

There's an example of a flat bottomed potential at http://docs.openmm.org/latest/userguide ... -container. It's a CustomExternalForce rather than a CustomBondForce, but the principle is the same.

You say the distances are often higher than the cutoffs. Is that a problem? In your force, r0 isn't an upper limit on the bond length. It's a lower limit on when the force is applied. The energy is zero for r<r0. The force doesn't even begin to have an effect until r>r0.

[Marcelo D Poleto]

Hi Peter,

I get it. Thank you! I think I was just surprised to see distances going 1A above r0 even when applying a constant of 1000kJ/mol/nm^2.

Is there any way to retrieve time series of the restraining forces being applied?

[Peter Eastman]

With that force constant, the energy of a 1 A displacement is 1000*0.1*0.1/2 = 5 kJ/mol. At 300K, kT is about 2.5 kJ/mol. So the energy is only 2 kT, well within the range of normal fluctuations.

Is there any way to retrieve time series of the restraining forces being applied?

If you query the configurations, you can compute the distance and restraining force yourself. Or if you put the CustomBondForce in its own force group, you can query the energy and forces from that group using the "groups" argument to Context.getState().

[Marcelo D Poleto]

Hi Peter,

Thanks for the input. I am currently running the code below to check the restraining energies:

Code: Select all

[creates system]

# Including flat-bottom restraint forces
flat_bottom_force = CustomBondForce('step(r-r0) * (k/2) * (r-r0)^2')
flat_bottom_force.addPerBondParameter('r0')
flat_bottom_force.addPerBondParameter('k')
flat_bottom_force.setName("Flat-bottom-force")
flat_bottom_force.setUsesPeriodicBoundaryConditions(True)

# 0-based list
# LIG-O4  = 13944
# LIG-CX2 = 13940
# E226-HN  = 4481
# G132-HN  = 2236
# S225-OG  = 4474
flat_bottom_force.addBond(13944, 4481, [1.0*angstrom, flat_constant*kilojoules_per_mole/(nanometer**2)])
flat_bottom_force.addBond(13944, 2236, [1.0*angstrom, flat_constant*kilojoules_per_mole/(nanometer**2)])
flat_bottom_force.addBond(13940, 4474, [1.7*angstrom, flat_constant*kilojoules_per_mole/(nanometer**2)])
system.addForce(flat_bottom_force)

fbout = open(jobname + "_restraints.dat", "w")

# Setting system
simulation = Simulation(psf.topology, system, integrator)
simulation.context.setPositions(crd.positions)

# Query restraint energies
for i, f in enumerate(system.getForces()):
	f.setForceGroup(i)
	if f.getName() == "Flat-bottom-force":
		fb_index = i
fb_energy = simulation.context.getState(getEnergy=True, groups={fb_index})
print(i, f.getName(), state.getPotentialEnergy()) # already giving 0.0 in return

[setup reporters]

! Iterate every step and recover restraint energies
for i in range(nsteps):
	simulation.step(1)
	fb_energy = simulation.context.getState(getEnergy=True, groups={fb_index})
	line = str(i).ljust(20) + str(fb_energy.getPotentialEnergy()) + "\n"
	fbout.write(line)

fbout.close()

However, fb_energy always returns 0.0. I tried decreasing r0 to a value that is definitely below the initial distance (half of the distance) so fb_energy > 0 but I still get a 0.0 energy in return. I double-checked and distances are being sampled above r0, so definitely fb_energy should be > 0.

I am not sure what I am doing wrong here, but I am fairly positive I am not getting the right result.

[Peter Eastman]

Can you provide a complete, self-contained test case that demonstrates the problem? Without that, I can't reproduce what's you're doing so I can only make guesses about what the problem might be.

Here are a couple of things I noticed from looking at the code.

What is the value of flat_constant? Make sure it isn't zero.

I notice you tell it to apply periodic boundary conditions. Make sure r > r0 after the boundary conditions are applied. Is it possible you're specifying two positions that are far apart, but periodic copies of the particles are closer?

[Marcelo D Poleto]

Hi Peter,

Sorry for the delay to reply on this. I am attaching below a tar ball with a toy model and the code I am using to run it. I cleaned up the code as much as I could. Please note the README file.

flat-bottom_restraint_example.zip

(729.11 KiB) Downloaded 7 times

Basically, the toy model is a alanine dipeptide system, in which I am trying to apply a flat-bottom restraint whenever the distance between HN-HA and HA-O are above 1 Angstrom. I know this is chemically inaccurate, but I want to make sure the restraint is being applied. Whenever I run this code/system, the flat-bottom energy printed is zero throughout the simulation.

I hope that can help illustrating the problem. Thank you very much for your help.

[Peter Eastman]

When I run your script it prints the message

- Initial flat-bottom energy: 32.670849030555445 kJ/mol

Is that not what you get? Which platform are you using?

[Marcelo D Poleto]

That is odd. When I run (with OpenMM 7.6 and a freshly installed 7.7), I get:

- Initial flat-bottom energy: 0.0 kJ/mol

The platform I am using is CUDA.

[Peter Eastman]

I see the problem. You call setForceGroup() on the forces after you create the Simulation. Any change to a System or the Forces it contains must be made before the Simulation/Context is created. Changes made after that are ignored unless you call reinitialize() on the Context. See https://github.com/openmm/openmm/wiki/F ... ons#system for more information about this.