[gmx-users] Harmonic COM Restraints with the Pull Code

[Justin Lemkul]

On 5/9/15 9:52 PM, Eric Dybeck wrote:
> Hello All,
>
> I am looking to run a set of simulations in which molecules are slowly
> restrained to their COM positions. To provide some context, I am going to
> melt a solid, and I would like to restrain the molecules first so that they
> remain in their lattice position as the melting takes places.
>
> I believe that I can achieve this using the pull code, although I've never
> used this feature in gromacs before. As I understand it, the settings in my
> .mdp file should look something like:
>
> ; Pull code
> pull            = umbrella
> pull_geometry   = distance
> pull_dim        = Y Y Y        ; 3D harmonic restraints
> pull_start      = yes           ; define initial COM distance > 0
> pull_ngroups    = 1          ; This may need to be the number of molecules
> pull_group0     =              ; Reference molecule
> pull_group1     =              ; Restrained Molecule
> pull_rate1      = 0          ; No pulling rate, I am just keeping them in a
> fixed position
> pull_k1         = 1000     ; This number I will vary along the alchemical
> path
>
> In most of the pull code examples that I've seen, one atom group in the
> simulation is being pulled away from another atom group also in the
> simulation. In this case, the groups in the simulation will need to get
> pulled toward a fixed 'reference' point that is defined based on an initial
> configuration that I can supply. Is this doable? Obviously I could also
> give an absolute (X,Y,Z) coordinate of the desired COM based on the initial
> .gro file if that is easier to do.
>

The reference coordinates can be constant with grompp -r, which overrides -c for 
the purpose of defining the reference for the restraint.  In your case, though, 
you're effectively trying to restraint a molecule/group to itself in some other 
state.  I doubt the pull code deals elegantly with this situation.  In theory, 
you could define a duplicate for each atom that is just a dummy, with some fixed 
position, that defines a reference for each pulled group.  Effectively, each 
molecule is a "ghost" of itself.  Then you can assign each restraint.  I have no 
idea if that will work; just a guess based on "typical" use of the pull code.

> Also, I will need to restrain all the molecules in the solid, rather than
> just a single molecule. Is it possible to have ~100 groups harmonically
> restrained with the Pull code?
>

Probably, but it's likely to be very inefficient.

If you're trying to melt a solid, wouldn't something like a sheer force using 
acc-grps be more appropriate?

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================