CustomCVForce crashes simulation with very high kinetic energy.

[Martin Floor]

I am trying to implement a fraction-of-native-contacts (Qf) force to apply it in umbrella sampling calculations. First, I am using a CustomBondForce to add each native contact contribution as:

Code: Select all

energy_function = '(1/N)*(1/(1+exp(beta*(r-lambda*r0))))'
Qf = CustomBondForce(energy_function)

Where N is the total number of native contacts and beta and lambda are global parameters for adjusting each contribution.

Then, I am passing this object as a collective variable (scalar) to a CustomCVForce using an Umbrella potential to restraint the Qf value as:

Code: Select all

QfCV = CustomCVForce('0.5*k*(Qf-Qf0)^2')

The remaining potential energy of the system is an implementation of a coarse-grained Go-like (structure-based) energy function that, otherwise, works fine. The problem is that the simulation is crashing with a "Potential energy is NaN" error, even with meagre constant forces or timesteps (within practical reason).

Just before crashing, I see the kinetic energy of the system getting very high and, in the trajectory, atoms are flying away in specific segments. I am not sure how the CustomCVForce object calculates the vectorial component of the force, mainly because it uses a scalar number as the driving force. But I am guessing that the problem is that the force applied is not restorative of the fraction of native contacts; kinetic energy is getting accumulated while the configuration stays far from the restrained Qf0 value.

Does anyone know how to deal with such a scenario, or maybe, has an alternative on how to implement this specific force?

[Peter Eastman]

I am not sure how the CustomCVForce object calculates the vectorial component of the force, mainly because it uses a scalar number as the driving force.

It just takes the gradient. The energy of the CustomCVForce depends on Qf, and that depends on the positions of individual atoms. Applying the chain rule gives the derivative of the energy with respect to atom positions.

What is the value of beta? Your custom bond force adds a contribution for each contact that varies from 0 to 1/N. The distance over which that transition happens is determined by beta. If it happens over a very short distance, that could lead to a very large force.

You also can directly query the forces to see what's happening. If you put the CustomCVForce into its own force group, you can query just that force on its own.

Code: Select all

customcv.setForceGroup(1) # Do this BEFORE creating the simulation
...
forces = simulation.context.getState(getForces=True, groups={1}).getForces()