[gmx-users] Error in free energy calculation

[Chris Neale]

Dear Sadaf:


You are starting with lambda = 10, which appears to be a state in which VDW
is off, but Coulomb is (partially) on. I’ve not used the FE code in gromacs
for quite some time (since it at least used to be unexpectedly slow
compared to regular simulations), but I suspect that this is your problem —
with q but no VDW you can get two charged atoms getting really close and
generating huge Coulombic forces, leading to system explosion. I believe
that it is normal to have q completely off before you start turning off the
VDW.


Other thoughts, not likely related to your crashes, but I include them
since I looked at your .mdp file already:


why not use nstlist = 10 ? I’m not sure why nstlist = 20 is ever OK with a
2 fs time step, but I realize that gromacs has made some recent strides in
NBlist building, so I might be misinformed here.


I get fine sd integrator stability with lincs-order = 6. Not sure how much
performance a lincs order of 12 is costing you, but you might be OK with 6
here.


Your tau_t is pretty small for velocity Langevin dynamics and you may be in
the over-damped regime. You may get more collective dynamics (and better
overall sampling) with a 1 ps value of tau_t.


Parrinello-Rahman pressure coupling, though giving a better ensemble than
Berendsen, is more fragile and will lead to exploding systems more easily
than Berendsen pressure coupling if your out of pressure/volume
equilibrium. You might at least try using Berendsen pressure coupling and
see if it fixes your issues.


You have gen_vel = no. Suggest to verify that you are in fact loading in
velocities.


Good luck,

Chris.