[gmx-users] Freeze groups, energy group exclusions and partially frozen dihedrals

[Mark Abraham]

Jennifer Williams wrote:
> I have two questions about Freeze groups
> 
> I have a large residue made up of SiO2 (like quartz). I also have 
> organic surface groups (imagine a propyl chain) which are anchored to 
> the surface of my large silica structure.
> 
> I have decided to freeze my silica structure (all silicons and oxygens) 
> in all 3 directions. I do this :
> 1. to save computational time,
> 2. because its movement in relation to the other moving parts of the 
> system is not significant
> 3. because although I have obtained some parameters for Si-O bending 
> etc-these were not parametrized for this particular system)
> 
> My topology file currently contains parameters for Si-O bending, O-Si-O 
> stretching etc (i.e those parts of the molecule I now wish to freeze). I 
> thought that if I used the freeze parameter, that these parameters would 
> be ignored and not contribute to the total potential energy. (Indeed 
> someone replied to my last post saying that this was the case and it 
> would be OK to leave them in). This is not the case as I have found 
> large differences in energy if I leave these parameters for the frozen 
> parts in or not. So I now assume I have to remove the parameters which 
> define completely frozen parts of my residue.

Energy differences are not strictly relevant to MD, since only the 
forces are used in the integration scheme. IIRC (but couldn't find on 
mailing list) freeze groups are implemented by zeroing the forces, so 
intra-freeze-group bonded parameters are indeed irrelevant. You save a 
small amount of computational time by omitting these bonded interactions 
- which is easy to do with an #ifdef.

> Question 1.
> 
> I have read that I should use energy group exclusions on the frozen part 
> of my molecule. However, I only achieve a neutral structure if I 
> consider the charge on both the frozen and non-frozen parts of my 
> molecule (I can?t change the charges to make the frozen part neutral). 
> In fact if I try to exclude the frozen group, I get the following warning.
> 
> WARNING 1 [file SIL.top, line 28748]:
>   Can not exclude the lattice Coulomb energy between energy groups
> 
> I assume this is because to ignore the charges on my frozen part would 
> leave me with a large overall charge.

No, this is because you're using PME. The long-range (i.e. lattice) 
component does not permit selective omission of components. Have a 
search on the mailing list for frozen groups - there've been heaps of 
threads. In particular, NPT is a questionable ensemble.

> If I continue with the energy group exclusion (but allowing the 
> coulombic energy to be calculated) I get:
> 
> Steepest Descents converged to Fmax < 200 in 144 steps
> Potential Energy  = -2.5177378e+05
> Maximum force     =  1.9338353e+02 on atom 4084
> Norm of force     =  1.5444931e+01
> 
> Without excluding the frozen energy group I get:
> 
> Steepest Descents converged to Fmax < 200 in 154 steps
> Potential Energy  = -1.4675300e+06
> Maximum force     =  1.7721352e+02 on atom 4084
> Norm of force     =  1.4064426e+01
> 
> Either one looks reasonable to me. I can?t decide which is the most 
> sensible option. Does it make more sense to use energy group exclusions 
> (with the exception of coulombic terms) or to go without excluding 
> groups at all.

Whatever gives you results that seem to correlate with reality. As a 
reviewer, I'd want to know why the non-silica parameter set was likely 
to work with the silica frozen (or unfrozen).

> Question 2
> 
> I am a bit mixed up about what to do when I have flexible surface groups 
> (i.e not frozen) anchored to a frozen surface-specifically about what to 
> do for constraints, bonds, angles and dihedrals which connect a frozen 
> molecule/s to moving one/s.
> 
> e.g If all Silicons and Oxygens are frozen, does it make sense to 
> include the following in my topology file.
> 
> [constraints]
> 
> Si ? C  (where Si is frozen)
> 
> [bonds]
> 
> Si ? C   (where Si is frozen)
> 
> 
> [angles]
> 
> 
> O ? Si - C  (Si and O  atoms frozen)
> 
> Si ? C - C  (Si  atom frozen)
> 
> 
> [dihedrals]
> 
> Si ? O ? Si - C   (3 atoms of dihedral frozen)
> 
> O ? Si ? C - C    (2 atoms of dihedral frozen)
> 
> Si ? C ? C - C     (1 atom of dihedral frozen)
> 
> My feeling is that I still need to include these parameters to define 
> how the flexible surface groups move. Does anyone know if this will give 
> spurious results because part of the molecule is frozen.?

Don't know, but I expect not. Berk's pretty good at making things work :-)

Mark

> I'd really appreciate any advice because this has been confusing me for 
> some time.
> 
> Thanks
> 
> 
> 
> 
>