[gmx-users] Reproducing results with independent runs

[alex.bjorling]

Dear users,

I am investigating protein crystal packing artifacts by doing equilibrium
simulations starting from a crystal structure. I would like to know if the
relaxations i see are reproducible, in the sense that many simulations with
independent velocities give the same general result.

My plans is to do only one set of (first NVT then NPT) equilibrations with
position restraints. Then, I thought I'd do a shorter NPT run with position
restraints, with more frequent output and using the trr snapshots as
starting points for production runs. 

The only question then is how far apart these snapshots need to be to
guarantee independent velocities. Attached is the velocity autocorrelation
for the Protein group. It seems to me that using snapshots 1ps apart would
do it, since the autocorrelation has decayed by then.

Is this a valid approach?

Best,
Alex

 

--
View this message in context: 
http://gromacs.5086.x6.nabble.com/Reproducing-results-with-independent-runs-tp5012372.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: how is the pulling force measured

[alex.bjorling]

2013/4/30 kim2811 [via GROMACS] 

>
> Hi Alex,
>
> I read the manual but I got confused with the details. My pull_geometry =
> position with pulling direction along y. I am trying to dissociate two
> interacting proteins. So I set the last amino acid of one protein as my
> reference group and the last amino acid of the other protein as the pulled
> group. Now, I don't know exactly how to interpret what reference group and
> pulled group mean, and what is their difference; since in the manual it
> says 'there is no difference in treatment of the reference and pulled group
> (except with the cylinder geometry).' Also, it says that measured force
> output is the force of the pulled group.
>

Ah. I think it would be wiser to use pull_geometry=distance, so that the
two groups (0 and 1, reference and pulled) are connected by a spring. It
might also be better to use groups bigger than one amino acid, but I guess
that depends on what you want to accomplish.

A very useful tutorial in which two interacting peptides are separated can
be found here:
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html


>
> So basically, I am trying to visualize how the system works. Using the
> pull_geometry stated above, are the two groups (reference and pulled)
> attached to two virtual springs? Will this statement be correct?
>

I'm no expert and will have to leave this question unanswered. Sorry.


>
> "Force was applied by moving the clamped ends of the two springs (i = 1,
> 2) in opposite directions with constant velocity v to positions Zi(t) =
> zi(0) + vt, where zi(0) is the initial position of the end of the molecule.
> The forces fi(t) at the two ends were measured at every time step using the
> relationship fi(t) = k (zi(t) - Zi(t))."
>
> I hope my questions make sense. Thank you.
>
> --
>  If you reply to this email, your message will be added to the discussion
> below:
>
> http://gromacs.5086.x6.nabble.com/how-is-the-pulling-force-measured-tp5007769p5007778.html
>  To unsubscribe from how is the pulling force measured, click 
> here
> .
> NAML
>




--
View this message in context: 
http://gromacs.5086.x6.nabble.com/how-is-the-pulling-force-measured-tp5007769p5007781.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: how is the pulling force measured

[alex.bjorling]

kim2811 wrote
> Hi,
> 
> I am not sure if I quite understand what reference group and pulling group
> mean in COM pulling. For the version GROMACS 4.X, I specified in my .mdp
> file the pull_ngroup = 1, and then I have my pull_group0 and pull_group1.
> Am I right in saying that the pull_group0 is fixed and only the
> pull_group1
> is pulled? Or are they simultaneously being pulled by two virtual springs?
> If this is the case, how is the force output measured?

That depends on other options, such as pull_geometry. All of this is really
well explained in the manual, section 7.3.21.

Best,
Alex



--
View this message in context: 
http://gromacs.5086.x6.nabble.com/how-is-the-pulling-force-measured-tp5007769p5007770.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: genbox box size problem

[alex.bjorling]

Justin Lemkul wrote
> I have seen this exact problem posted before.  The coordinates of the
> water box 
> provided by MARTINI have some coordinate values that basically coincide
> across 
> the periodic box.  With genconf, you can apply the -dist option to space
> them 
> out slightly, then re-equilibrate before using the new configuration as
> input 
> into genbox.

I don't think that's true. Stacking boxes with -nbox 3 3 3 (without -dist)
gives a big box which you can minimize and equilibrate with no trouble. Both
with and without -dist 0.21 0.21 0.21 (0.21 is a suitable MARTINI -vdwd for
genbox), you get big water boxes which behave well. The problems start when
you use them as solvents in genbox.

Xavier's suggestion (stack, minimize, equilibrate, manually cut with a
script) works fine, thanks for the suggestion.

Seems to me there might be a bug in genbox, but I'm not the right person to
make such claims.

Alex



--
View this message in context: 
http://gromacs.5086.x6.nabble.com/genbox-box-size-problem-tp5007753p5007768.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: pull=constraint gives zero forces

[alex.bjorling]

Sorry - forgot to mention that before crashing, the run with all other
constraints removed produces a single line of pullf output:

0.  -812.401-4002.84482.04  1951.47 138.953 -1806.55   
-601.0072644.79 447.018 1768.6  -214.64 -199.829-2746.97   
1177.7  476.39  288.535 -559.274123.08  114.493 851.86  550.558

As Thomas Schlesier mentions here,
http://gromacs.5086.n6.nabble.com/pull-constraint-gives-zero-forces-tp5001817.html,
the pullf output apparently contains the forces necessary to enforce the
constraints.

/ Alex


alex.bjorling wrote
> Thanks guys,
> 
> Fixing the bug, recompiling and trying again results in a segfault
> like with version 4.0.7. I interpret this as GROMACS working fine, and
> suppose that there's something wrong with my input. Will continue this
> thread here and would really appreciate any ideas on how to proceed.
> 
> Before the segfault, I get a bunch of LINCS warnings, all concerning
> atoms that have other constraints in the topology. Manually replacing
> these by stiff bonds in the itp gets rid of the LINCS warnings, but
> still produces an immediate segfault.
> 
> The complete mdp follows. (Chris: previously posted via the web forum
> - it seems then you have to click the nabble link to see it).
> 
> Cheers,
> Alex
> 
> 50_constr3.mdp:
> **
> pull= constraint
> pull_geometry   = position
> pull_dim= Y Y Y
> pull_start  = yes
> pull_group0 = BB
> pull_nstxout= 1000
> pull_nstfout= 1000
> pull_ngroups= 7
> pull_constr_tol = 1e-6
> 
> pull_group1 = group1
> pull_init1  = 0.0 0.0 0.0
> pull_rate1  = 0.0 0.0 0.0
> 
> pull_group2 = group2
> pull_init2  = 0.0 0.0 0.0
> pull_rate2  = 0.0 0.0 0.0
> 
> pull_group3 = group3
> pull_init3  = 0.0 0.0 0.0
> pull_rate3  = 0.0 0.0 0.0
> 
> pull_group4 = group4
> pull_init4  = 0.0 0.0 0.0
> pull_rate4  = 0.0 0.0 0.0
> 
> pull_group5 = group5
> pull_init5  = 0.0 0.0 0.0
> pull_rate5  = 0.0 0.0 0.0
> 
> pull_group6 = group6
> pull_init6  = 0.0 0.0 0.0
> pull_rate6  = 0.0 0.0 0.0
> 
> pull_group7 = group7
> pull_init7  = 0.0 0.0 0.0
> pull_rate7  = 0.0 0.0 0.0
> 
> ;
> ; STANDARD MD INPUT OPTIONS FOR MARTINI 2.0 or 2.1
> ;
> 
> ; TIMESTEP IN MARTINI
> ; Most simulations are numerically stable
> ; with dt=40 fs, some (especially rings) require 20-30 fs.
> ; Note that time steps of 40 fs and larger may create local heating or
> ; cooling in your system. Although the use of a heat bath will globally
> ; remove this effect, it is advised to check consistency of
> ; your results for somewhat smaller time steps in the range 20-30 fs.
> ; Time steps exceeding 40 fs should not be used; time steps smaller
> ; than 20 fs are also not required.
> 
> ;define = -DPOSRES
> integrator   = md
> tinit= 0.0
> dt   = 0.02
> nsteps   = 250 ; 50 ns
> nstcomm  = 1
> comm-grps  =
> 
> nstxout  = 0
> nstvout  = 0
> nstfout  = 0
> nstlog   = 1000
> nstenergy= 100
> nstxtcout= 1000
> xtc_precision= 100
> xtc-grps =
> energygrps   = Protein
> 
> ; NEIGHBOURLIST and MARTINI
> ; Due to the use of shifted potentials, the noise generated
> ; from particles leaving/entering the neighbour list is not so large,
> ; even when large time steps are being used. In practice, once every
> ; ten steps works fine with a neighborlist cutoff that is equal to the
> ; non-bonded cutoff (1.2 nm). However, to improve energy conservation
> ; or to avoid local heating/cooling, you may increase the update frequency
> ; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
> ; is computationally less expensive and leads to improved energy
> conservation
> 
> nstlist  = 10
> ns_type  = grid
> pbc  = xyz
> rlist= 1.4
> 
> ; MARTINI and NONBONDED
> ; Standard cut-off schemes are used for the non-bonded interactions
> ; in the Martini model: LJ interactions are shifted to zero in the
> ; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2
> nm.
> ; The treatment of the non-bonded cut-offs is considered to be part of
> ; the force field parameterization, so we recommend not to touch these
> ; values as they will alter the overall balance of the force field.
> ; In principle you can i

[gmx-users] Re: pull=constraint gives zero forces

[alex.bjorling]

Following up on this post. I've tried the same runs using version 4.0.7,
which gave immediate segmentation faults. Not sure if this is a clue or a
trivial consequence of switching versions, but there it is.

Any other ideas why the pullf output just contains zeros?

Cheers,
Alex


alex.bjorling wrote
> I'm using the pull code to maintain the initial structure of a protein
> that otherwise deforms. Using pull=umbrella does what I expect it to,
> but pull=constraint produces zero forces. I'm using version 4.5.5 with
> the MARTINI force field.
> 
> The pull=umbrella mdp contains the following,

> 
> and gives the following pullf output:

> 
> For such runs, the group COM:s stay within 0.1 nm of their initial
> positions, throughout a long trajectory.
> 
> The pull=constraint mdp starts with

> 
> and produces the pullf output

> 
> For these runs, the group COM:s move around freely. I guess I'm doing
> something wrong but can't work out what. I've tried specifying k:s for
> the contraint runs, and tried removing all other constraints in the
> molecule. I've tried to comply with the manual's instructions but to
> no avail.
> 
> Any ideas?
> 
> Cheers,
> Alexander Björling,
> PhD candidate, University of Gothenburg, Sweden





--
View this message in context: 
http://gromacs.5086.n6.nabble.com/pull-constraint-gives-zero-forces-tp5001538p5001793.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] pull=constraint gives zero forces

[alex.bjorling]

Hi,

I'm using the pull code to maintain the initial structure of a protein
that otherwise deforms. Using pull=umbrella does what I expect it to,
but pull=constraint produces zero forces. I'm using version 4.5.5 with
the MARTINI force field.

The pull=umbrella mdp contains the following,


and gives the following pullf output:


For such runs, the group COM:s stay within 0.1 nm of their initial
positions, throughout a long trajectory.

The pull=constraint mdp starts with



and produces the pullf output



For these runs, the group COM:s move around freely. I guess I'm doing
something wrong but can't work out what. I've tried specifying k:s for
the contraint runs, and tried removing all other constraints in the
molecule. I've tried to comply with the manual's instructions but to
no avail.

Any ideas?

Cheers,
Alexander Björling,
PhD candidate, University of Gothenburg, Sweden




--
View this message in context: 
http://gromacs.5086.n6.nabble.com/pull-constraint-gives-zero-forces-tp5001538.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists