[gmx-users] Re: Wierd results from Umbrella sampling, (Justin A. Lemkul)

[Justin A. Lemkul]

On 5/22/12 6:53 PM, Thomas Schlesier wrote:
> I never worked with the MARTINI (or other coarse-grained) force field, but this
> in the umbrella.mdp
>
> title = Umbrella pulling simulation
> integrator = md
> dt = 0.019
>
> looks suspicious. The dt is about an order of magnitude greater than one uses in
> normal (bond-)constrainted md-simulations. I know that in coarse-grained
> simulations the potentials are smoother than in atomistic force fields and one
> therefore could use a higher timestep, but i don't know if you can go so high.
> Anyway you should look for the umbrella potential. If the force constant and the
> timestep are both too high you could get problems:
> Assume a particle in a harmonic well. If the force constant is high and the
> timestep too high, you wont sample the harmonic well, but 'jump' each time from
> one side to the other (high force + great timestep -> long movement).
> So i would test it first with a timestep of 0.002 (same as you used for pulling
> sim).
>

Testing this would be good, though the MARTINI developers routinely report 
timesteps of 20-40 fs as "normal" use.  I have never obtained stable 
trajectories above 20 fs, but I also do not do much coarse graining.  It could 
very well be that the pull code is not stable with such an integration step.

-Justin

>
> Am 22.05.2012 18:37, schrieb gmx-users-request at gromacs.org:
>> Thank you for your reply, Thomas. Under your explanation, I am well understood
>> of the terms: pull_k and pull_rate.
>> Here I upload both md_umbrella.mdp and md_pulling.mdp.
>> I have to mention that it is coarse gained system with Martini force field.
>>
>> At the same time, i am going to run a simulation without pull code but with
>> LINCS constraint, and also run another system with a huge pull_k but without
>> LINCS. Hope I could get some interesting information.
>>
>> With Regards,
>>
>> Jiangfeng.
>>
>> Message: 3
>> Date: Tue, 22 May 2012 16:36:47 +0200
>> From: Thomas Schlesier<schlesi at uni-mainz.de>
>> Subject: [gmx-users] Re: Wierd results from Umbrella, sampling (Justin
>> A. Lemkul)
>> To:<gmx-users at gromacs.org>
>> Message-ID:<4FBBA47F.5010503 at uni-mainz.de>
>> Content-Type: text/plain; charset="ISO-8859-1"; format=flowed
>>
>> Think it would be best to show the .mdp file, else we can only guess
>> what the parameters are.
>> From the histogram it looks like that the force constant of the
>> restraining potential is too low, since the histograms are really wide,
>> but pull_k1=1000 is a 'normal' value.
>> On thing which confueses me is you said that fluctuations from g_dist
>> are about 0.4nm, but in the histogram i looks like the distance
>> fluctuates about 1nm or even more. for this be sure, that you compare
>> the right distances -> if you pull only in z-direction, the only look
>> into the z-direction from the output from g_dist. Since the x, and
>> y-directions would be unaffected by the umbrella potential.
>>
>> for the error message with the constraints:
>> what happens if you run the system with constraints but without the
>> pull-code?
>>
>> for pull_k1>0 and pull_rate1=0:
>> if you're pulling with an umbrella potential pull_k1 defines the force
>> constant of the potential (hook'sches law). Let's assume you put a
>> spring to your molecule
>> with pull_rate1=0, it means that the other end of the spring doesn't
>> move, and the spring restrains the position of the molecule (molecule
>> cn't diffuse away). in umbrella sampling you want to restrian of
>> molecule (or part of it) relative to another one / part -> so pull_rate1=0
>> with pull_rate1>0, it's like you pull the spring away for the molecule,
>> spring gets strechted -> wants to relax -> molecule follows spring (like
>> that what happens in afm pulling)
>>
>>
>> Am 22.05.2012 15:40, schriebgmx-users-request at gromacs.org:
>>> > Hi Justin, As for the maximum of 0.4 nm fluctuation of my pulled
>>> > protein, I used g_dist to calculate the distance between my pulled
>>> > protein and stable part in a window, where the distance is treated as
>>> > the fluctuation of the protein along z-axis. Well, from pullx.xvg, the
>>> > position moved a lot (3nm for instance.) As for the windows simulation,
>>> > I didn't apply constraint but only the internal constraint in the itp
>>> > file. I still don't understand why it have to do constraint? why not
>>> > give a fully freedom to run the simulation? In the md_umbrella.mdp, I
>>> > set pull_k1=1000; pull_rate1= 0.0, but apparently, I am confused with
>>> > pull_k1>0 combined with pull_rate1=0. In my mdp, i set "none" to LINCS,
>>> > because if I use "all-bonds", an error of "1099 constraints but degrees
>>> > of freedom is only1074" occurs. Actually, there is no any window with a
>>> > designed distance. Here I attach the histo.xvg, profile.xvg. Thank you
>>> > with regards, Jiangfeng. On 5/22/12 9:36 AM, Du Jiangfeng (BIOCH) wrote:
>>>>> >> > Dear Justin,
>>>>> >> >
>>>>> >> > Based on your questions to my simulation, I posted here yesterday
>>>>> hopefully
>>>>> >> > it was the correct way to reply in this forum.
>>>>> >> >
>>> > You've still replied to the entire digest message (which I've cut out); please
>>> > make sure to keep replies free of superfluous posts in the future. The archive
>>> > is already pretty hopeless, but let's not make it worse:)
>>> >
>>>>> >> > In this morning I got a list of new windows of umbrella sampling, the
>>>>> overlap
>>>>> >> > is sufficient enough, but I saw another problem: In the histogram figure,
>>>>> >> > the base of peak covers the distance of 2 nm instead of 0.2 nm., that's
>>>>> >> > horrible! However, when I checked back to the simulation results of each
>>>>> >> > window, the fluctuation of my pulled protein is only 0.4nm in maximum.
>>>>> So the
>>>>> >> > base of peak shouldn't cover such long distance, right?
>>>>> >> >
>>> > If the peaks aren't corresponding to the desired restraint distances, then
>>> there
>>> > are several potential problems:
>>> >
>>> > 1. Your restraints aren't set up the way you think they are (check grompp
>>> output
>>> > and .tpr file contents to be sure)
>>> > 2. Your restraints are ineffectual (in which case you may need to revisit the
>>> > force constant)
>>> >
>>> > I can't determine from your description what's going on. What do you mean by a
>>> > maximum of 0.4 nm fluctuation? In what quantity? What do the contents of
>>> > pullx.xvg show you for the problematic window, and for that matter, the
>>> others?
>>> > Are any of them producing the desired restraint distances?
>>> >
>>> > Again I will ask you to share an image of the histogram and PMF profile; these
>>> > would be very helpful to see.
>>> >
>>> > -Justin
>>> >
>>> > --
>>> > ========================================
>

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

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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