Page 1 of 2
Flat-bottom potential
Posted: Fri Jul 22, 2022 9:07 am
by mdpoleto
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)])
Re: Flat-bottom potential
Posted: Fri Jul 22, 2022 9:50 am
by peastman
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.
Re: Flat-bottom potential
Posted: Fri Jul 22, 2022 11:56 am
by mdpoleto
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?
Re: Flat-bottom potential
Posted: Fri Jul 22, 2022 12:34 pm
by peastman
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().
Re: Flat-bottom potential
Posted: Wed Jul 27, 2022 11:35 am
by mdpoleto
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.
Re: Flat-bottom potential
Posted: Wed Jul 27, 2022 11:43 am
by peastman
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?
Re: Flat-bottom potential
Posted: Tue Aug 09, 2022 7:55 am
by mdpoleto
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.
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.
Re: Flat-bottom potential
Posted: Tue Aug 09, 2022 8:49 am
by peastman
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?
Re: Flat-bottom potential
Posted: Tue Aug 09, 2022 9:11 am
by mdpoleto
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.
Re: Flat-bottom potential
Posted: Tue Aug 09, 2022 10:29 am
by peastman
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.