[gmx-users] System equilibration and interaction energy

[Mark Abraham]

aloks at mbu.iisc.ernet.in wrote:
> Dear Gromacs users,
> 
> I have two queries :
> 1. I am doing simulation for a protein (156 residues ) with a sugar. Out
> of 156 protein resiudes i have restrained 130 residues and rest part  i
> have kept free. My doubt is that when we can say that system has
> equilibrated?

You will never be able to say the system is equilibrated - only that it 
may be equilibrated subject to the restraints.

>    if on the basis of rmsd then i dont feel it would be advisable to look the
> rmsd of protein because major part i have restrained, if i consider  sugar
> rmsd then i can see a kind of periodicity in the rmsd. if i look the rmsd
> of the system then it seems to be stablized at 4.5nm.
>   if i look the potential energy of the system then from the begining
> itself it seems to be equilibrated with a fluctuation of ~700kj/mol. I
> guess that may be because of very high contribution from the solvent.
> 
>    For nearly 26000 atoms i have done the simulation upto 4 ns. So i would
> just like to know that what is the best way to define system
> equilibration when a major part is restrained.

You have the same methods available to guide your guesses as for 
unrestrained MD - check that observables like total energy, density, 
position RMSD (for systems expected to be in a PE minimum) are doing 
sensible things. Since the information content of the restrained atoms 
is low, you may like to exclude them from the maths.

> 2. my second question is about caluclation of the interaction energy
> betwen the protein and the ligand. one possible way seems to be calculate
> all the energy pairs given by gromacs with ligand and protein (ljsr, ljlr,
> lj14 and for coulombs) and sum it or do i need to do two more simulations,
> one with protein alone in the same box size with same number of water
> molecules and other with the ligand alone and then substracting the sum of
> potential energies of independent simulations from the protein-ligand
> energy?

Your second method would be totally weird. If you are using an NVT 
ensemble, you will now have a different density, so your solvent will be 
doing different things interacting with itself and the remaining solute 
such that there's no reason to expect energy additivity. If NPT, the 
system will contract and there's still no such reason. I believe the 
first method is the accepted practice. The interaction energy is purely 
through-space in the absence of bonding, and so only the LJ and Coulomb 
terms contribute... thus they are all you need to consider.

Mark