[gmx-users] Fwd: self-assembly and box size

[Justin Lemkul]

On 11/14/16 9:04 PM, Mohsen Ramezanpour wrote:
> Thanks Justin for your comments.
>
> Comments interspersed.
>
> On Mon, Nov 14, 2016 at 6:47 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 11/14/16 5:11 PM, Mohsen Ramezanpour wrote:
>>
>>> Dear Gromacs users,
>>>
>>> I have made an initial system with a number of molecules in a "large" box.
>>>
>>> Large means a box of size 50  50  50 nm?
>>>
>>> I am interested in the self-assembly process of molecules (parametrized in
>>> charmm36 ff) in the box.
>>>
>>> I expected pressure coupling to shrink the box quickly, but it seems it
>>> fails to do so. The simulation box did not change (ca. 0.5 nm in each
>>> direction) which seems strange to me.
>>>
>>> I have used four different P-couple schemes for 4 different type of
>>> self-assemblies:
>>>
>>> 1) isotropic
>>>
>>> pcoupl                  = berendsen
>>> pcoupltype           = isotropic
>>> tau_p                   = 5.0
>>> compressibility     = 4.5e-5
>>> ref_p                    = 1.0
>>>
>>>
>>> 2) semiisotropic
>>>
>>> pcoupl                 = berendsen
>>> pcoupltype           = semiisotropic
>>> tau_p                   = 5.0
>>> compressibility     = 4.5e-5  4.5e-5
>>> ref_p                   = 1.0     1.0
>>>
>>>
>>> 3) anisotropic with *fix shape* of box (*but the size could change*)
>>>
>>> pcoupl                  = berendsen
>>> pcoupltype            = anisotropic
>>> tau_p                    = 5.0
>>> compressibility      = 4.5e-5  4.5e-5  4.5e-5   0       0       0
>>> ref_p                     = 1.0      1.0       1.0       0       0       0
>>>
>>> 4) anisotropic with *both size and shape able to change*
>>>
>>>
>>> pcoupl                  = berendsen
>>> pcoupltype           = anisotropic
>>> tau_p                   = 5.0
>>> compressibility     = 4.5e-5  4.5e-5  4.5e-5  4.5e-5  4.5e-5  4.5e-5
>>> ref_p                   = 1.0     1.0     1.0     0       0       0
>>>
>>>
>>> Here are some other parameters in mdp file which might be useful for
>>> finding the problem:
>>>
>>> integrator              = md
>>> dt                         = 0.002
>>> nsteps                  = 50000
>>> nstlog                   = 1000
>>> nstxout                 = 1000
>>> nstvout                 = 1000
>>> nstfout                  = 1000
>>> nstcalcenergy       = 100
>>> nstenergy             = 1000
>>>
>>> cutoff-scheme      = Verlet
>>> nstlist                  = 20
>>> rlist                     = 1.2
>>> coulombtype        = pme
>>> rcoulomb             = 1.2
>>> vdwtype               = Cut-off
>>> vdw-modifier        = Force-switch
>>> rvdw_switch        = 1.0
>>> rvdw                    = 1.2
>>>
>>> tcoupl                   = berendsen
>>> tc_grps                 = Molecules   water_ions
>>> tau_t                     = 1.0    1.0
>>> ref_t                      = 355    355
>>>
>>> constraints             = h-bonds
>>> constraint_algorithm    = LINCS
>>> continuation            = yes
>>>
>>> nstcomm                 = 100
>>> comm_mode            = linear
>>> comm_grps             = Molecules   WI
>>>
>>>
>> You shouldn't reset COM motion like this.
>
>
> Do you mean deviding it to Molecules and WI? If yes, could you please let
> me know why it is not approperiate?
> How if I use "system" instead?
>

Use System.  Only layered systems should be considered for separate comm_grps. 
Your approach risks the different groups having spurious contributions from COM 
motion removal causing them collide.  I've seen crashes from it before.

>>
>>
>> refcoord_scaling        = com
>>>
>>> gen-vel                 = no
>>>
>>> I can make the initial system smaller, but I think starting from a large
>>> box would be better because it gives all the molecule enough freedom to
>>> move and assemble in better way.
>>>
>>> Is not it right to expect the p-couple shrink the box in nanoseconds? This
>>> was only for 100 ps but I expected more shrinkage of box. What do you
>>> think?
>>>
>>>
>> 100 ps is nothing.
>>
>
> Even if the box is not full with waters and there are a lot of empty spaces
> between molecules?
>

"A lot" is subjective, but in general, I stand by the statement that 100 ps is 
nothing in this context.

>>
>> If you have a box full of water, as your settings above imply,
>
> That is not the case, I have molecules and water, but there are plenty of
> empty spaces.
>

What is the purpose of void space?  You're potentially introducing a lot of 
artificiality.  Perhaps you need to better explain your goals.  Using such a 
huge box is computationally very expensive, for potentially little or no gain, 
plus the headaches associated with bizarre behavior due to voids (which can also 
cause crashes).

>
>> you can't expect much of a change.  Why do you need the box to shrink, or
>> why do you expect it to?
>
>
> Based on the number of lipid molecules and waters I added to the box, and
> since the same system (with the same number of molecules and waters) will
> be in lamellat  or other phases (in equilibrium) of smaller sizes (e.g. 10
>  10  10 nm). i.e. if I make a bilayer out of these, it will be ca. 10*10*10
> nm^3 in size
>
>
>> A box will only shrink a lot if there is substantial void space, e.g.
>> intentional or due to a poor initial approximation of the system that
>> requires equilibration.
>>
> Agreed. Regarding what I described about the system, do you think this is
> because 100 ps is not enough? Even for a small shrink of ~ 2 nm?
>

Impossible to predict.  But again, 100 ps is not nearly enough time to assess 
anything like what you're after.

-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

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