[gmx-users] Freezing Groups

[Mark Abraham]

Reza Salari wrote:
> Hi everybody,
> 
> I want to do a simulation using frozen groups and it seems it is much 
> more difficult than I first thought. I searched the mailing list and I 
> still have problems that I can't figure out. It is much appreciated if 
> help me with these problems or point me in the right direction.
> 
> The scenario is this: I have a protein and a ligand and I want to do 
> free energy calculation. I want to freeze the protein (for a couple of 
> reasons but mainly because it is a big protein and  I want to save 
> computation time so I don't want to use restraints) and I know these 
> about freezing the protein (please correct me if any of them is incorrect):
> -- For minimization, it seems only steep algorithm works.
> -- I should not use pressure coupling and any constraints (bonds,...) 
> for frozen groups.
> -- In my mdp file for production run, I defined these:
>     coulombtype=pme
>     energygrps=protein ligand SOL
>     energygrp_excl=protein protein
>     freezegrps=protein
>     freezedim= Y Y Y
>     comm_mode=no
> 
> -- I know that I should use GROMACS 4.0.2 with fixed freeze code.
>    
> Now when I want to run the grompp for production mdp, I get this warning:
>     " can not exclude the lattice coulomb energy between energy groups"
> 
> I know that this is because I used pme and it seems that I have only two 
> options to solve it: either use other types of long range electrostatics 
> that don't use lattice sum (which I don't want to, because I think PME 
> is much better than others; also I have important electrostatic 
> interactions at protein-ligand interface so I want to use the best 
> method for electrostatics) or use [ exclusions ] part of topology file 
> to exclude interactions of every pair of atoms in the protein (I saw 
> this solution in the mailing this, which is really cumbersome in my case 
> because the protein is big, but I'll do that if I have to).

That looks like a perfect excuse to learn how to write a shell or other 
script. Off the top of my head, I think all you need is a loop over the 
range of protein atom numbers written to a file that you can #include in 
your main topology.

> Here are my questions:
> 1) Isn't there any other solution for this problem?

Not that I know of.

> 2) Is that warning really critical? I mean if I am going to compare the 
> two states of the same protein (same protein with different ligands), 
> don't those "unexcluded coulomb energies" cancel out each other when 
> comparing the total energies of the two states? Or there is no guarantee 
> for this?

True, the mesh contribution of a frozen group should be constant if the 
box stays the same size (OK, it is PME, approximately constant), and 
will be approximately constant if the box does vary a little. Thus the 
error introduced in comparing total energies is small. Since there are 
no forces on the protein atoms, I think that the mesh contribution to 
forces will have no appreciable effect, thus the dynamics would be OK. 
However, the [exclusions] solution is so simple it seems poor economy 
not to use it.

Mark