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.