michellab / BioSimSpace

Code and resources for the EPSRC BioSimSpace project.
https://biosimspace.org
GNU General Public License v3.0
77 stars 19 forks source link

`minimize=True` in free energy protocol #233

Closed kexul closed 3 years ago

kexul commented 3 years ago

Dear BioSimSpace developers,

Recently I found that the default protocol generated by BSS.FreeEnergy.Relative() set minimise=True in somd.cfg. When I passed an already equilibrated system (by Gromacs) to it, the energy would drop a lot during minimization.

I'm concerned that:

  1. Would that makes the system become unequilibrated again?
  2. Does that mean my system is not well minimized? Before passing the system to BSS.FreeEnergy.Relative(), I performed Minimise->NVT equilibration->NPT equilibration by Gromacs. The energy before minimization agrees well with the potential energy printed by Gromacs.
###=======================Minimisation========================###
Running minimisation.
Energy before the minimisation: -408467 kcal mol-1
Tolerance for minimisation: 1
Maximum number of minimisation iterations: 1000
Energy after the minimization: -519065 kcal mol-1
Energy minimization done.
###===========================================================###

Any help is appreciated! 🤗

Here is the full log of somd-freenrg. somd.txt

jmichel80 commented 3 years ago

this is default practice because the equilibrated input is implicitly equilibrated at one lambda value (even though lambda-dependent potentials are not used during equilibration), typically lambda=0.0 (but that depends on how your put together the equilibration-->free energy pipeline). Changing lambda changes interaction with the system, and can introduce steric clashes (depending on the perturbation) that will blow up the MD integrator with high magnitude forces.

Therefore we recommend to energy minimise the input at each lambda value, and discard the beginning of the output trajectories (to allow for requilibration) before free energy estimation. Typically 15-50 ps seems sufficient but that could depend on the exact use case.

kexul commented 3 years ago

Thanks, that makes things much clearer now!