[gmx-users] Free energy of protonation

[Gmx QA]

Dear list,

I am trying to calculate the free energy for protonating a single monomer
of oleic acid with the martini forcefield. This in preparation for other
work, and I wanted to start by reproducing what is already done.

To do this, I have defined the following in my topology to describe the
change in protonation state going from A (protonated) to state B (charged).
The martini bead type needs to change as well in this process.

[ atoms ]

; idtypeA   resnr   residu  atomcgnrchargeA   massA typeB
chargeB massB

  1 P4  1   OAN COO 1   0   Qa  -1

I also have this very simple mdp file for doing the transformation

free_energy  = yes

init_lambda-state= 1

fep-lambdas  = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0


And I am using fep-lambdas to describe the entire process. With this, I can
run through grompp and mdrun with different values of init-lambda-state,
but the output I get in each case is the same for each different simulation.


I check this by


gmx analyze -f md.xvg -ee


And I always get the same value for the average. Is there something I am
missing or not specifying properly?


Thanks

/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] RDF between larger molecules

[Gmx QA]

Thanks Justin!

Can you also comment on the other issue, the lack of smoothness? I guess
this is because I don't have that many individual coms per bin to
accurately compute the rdf? At least I think that makes sense since when I
increase the value of -bin in gmx rdf I get an increase in "smoothness" and
fewer "jump" back and forth.



Den mån 20 jan. 2020 kl 15:50 skrev Justin Lemkul :

>
>
> On 1/20/20 9:48 AM, Gmx QA wrote:
> > Hi list,
> >
> > I am working on a system that contains lipids, but since there is no
> water
> > present they do not form a typical bilayer.
> >
> > In these systems I also have some smaller drug-like molecules. What would
> > be the best way to get a proper rdf-function between the lipids and the
> > drugs? Since both types of molecules are larger than single atoms, I
> tried
> > to calculate a com-com rdf using
> >
> > -selrpos whole_mol_com -seltype whole_mol_com
> >
> > but the resulting curve is very rugged (i.e. not as smooth as
> > the atom-rdf's typically shown). The minimum distance (the distance
> > below which the pdf is zero) also appears to be quite small, since one
> the
> > molecules is a lipid I would have expected a larger distance than
> > the corresponding value seen in e.g. an oxygen-oxygen rdf for water.
> >
> > Might all of this be consequences of the fact that I have a com-com rdf
> > between a limited number of molecules?
>
> Yes, likely. You could in principle have a COM-COM RDF value very close
> to r=0 if a lipid were to wrap around a drug molecule, so I don't think
> any of your results are unexpected.
>
> -Justin
>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> Assistant Professor
> Office: 301 Fralin Hall
> Lab: 303 Engel Hall
>
> Virginia Tech Department of Biochemistry
> 340 West Campus Dr.
> Blacksburg, VA 24061
>
> jalem...@vt.edu | (540) 231-3129
> http://www.thelemkullab.com
>
> ==
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] RDF between larger molecules

[Gmx QA]

Hi list,

I am working on a system that contains lipids, but since there is no water
present they do not form a typical bilayer.

In these systems I also have some smaller drug-like molecules. What would
be the best way to get a proper rdf-function between the lipids and the
drugs? Since both types of molecules are larger than single atoms, I tried
to calculate a com-com rdf using

-selrpos whole_mol_com -seltype whole_mol_com

but the resulting curve is very rugged (i.e. not as smooth as
the atom-rdf's typically shown). The minimum distance (the distance
below which the pdf is zero) also appears to be quite small, since one the
molecules is a lipid I would have expected a larger distance than
the corresponding value seen in e.g. an oxygen-oxygen rdf for water.

Might all of this be consequences of the fact that I have a com-com rdf
between a limited number of molecules?

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Hydrogen bond autocorrelation

[Gmx QA]

Hi David

Thanks for your email. In this case, I explicitly select only the compound
itself when computing hydrogen bonds, so while I see your point about
water-mediated h-bonds, the results (including the act-function) are the
same regardless of -nomerge.

Could you comment on whether the h-bond computed in gromacs corresponds to
the continuous or the intermittent definition of a h-bond? I think the
latter, but I am not sure. And are there any consequences later that would
depend on which h-bond acf is used, such as for example when fitting an
exponential function for determining the life time of the (single) bond?

Thanks


Den ons 19 juni 2019 kl 08:26 skrev David van der Spoel <
sp...@xray.bmc.uu.se>:

> Den 2019-06-18 kl. 15:58, skrev Gmx QA:
> > Dear list,
> >
> > This has been discussed many times previously on the list, but I still
> have
> > some questions about hydrogen bond autocorrelation functions.
> >
> > I have run a simulation with single molecule of a compound in water which
> > can form exactly one inter-molecular hydrogen bond. From the hbnum.xvg
> > time-series (i.e. the existence function) of this h-bond I then have a
> > python-script to calculate the acf. The script seems to work because the
> > resulting acf-function looks exactly the same as the corresponding
> function
> > computed with xmgrace. But it is very different from what I get from gmx
> > hbond -ac.
> >
> > I read a lot about continuous vs intermittent acfs for h-bonds, the
> latter
> > being called non-continuous in the van der Spoel et al. 2006-paper I
> think.
> > Is the reason for the discrepancy simply that gmx hbond -ac calculates a
> > _different_ acf than you get from the hbnum.xvg file (and xmgrace)?
> >
> > Then, going forward, after I compute the acf and in the most basic case
> fit
> > to y = a* exp(-t/b), does the value of b correspond to the life-time of
> the
> > h-bond? Or it is the integral (from 0 to inf) that gives you the
> lifetime?
> > This seems to be treated differently by different people.
> >
> > Thanks a bunch
> > /PK
> >
> There are a couple of details indeed, in particular for water. Water may
> bind with oxygen or hydrogens. When your compound provides a hydrogen
> bond acceptor, gmx hbond will compute the hydrogen bond between your
> compound and both H on a water molecule as the same H bond, unless you
> tell it not to with an option -nomerge.
>
> --
> David van der Spoel, Ph.D., Professor of Biology
> Head of Department, Cell & Molecular Biology, Uppsala University.
> Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
> http://www.icm.uu.se
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Hydrogen bond autocorrelation

[Gmx QA]

Dear list,

This has been discussed many times previously on the list, but I still have
some questions about hydrogen bond autocorrelation functions.

I have run a simulation with single molecule of a compound in water which
can form exactly one inter-molecular hydrogen bond. From the hbnum.xvg
time-series (i.e. the existence function) of this h-bond I then have a
python-script to calculate the acf. The script seems to work because the
resulting acf-function looks exactly the same as the corresponding function
computed with xmgrace. But it is very different from what I get from gmx
hbond -ac.

I read a lot about continuous vs intermittent acfs for h-bonds, the latter
being called non-continuous in the van der Spoel et al. 2006-paper I think.
Is the reason for the discrepancy simply that gmx hbond -ac calculates a
_different_ acf than you get from the hbnum.xvg file (and xmgrace)?

Then, going forward, after I compute the acf and in the most basic case fit
to y = a* exp(-t/b), does the value of b correspond to the life-time of the
h-bond? Or it is the integral (from 0 to inf) that gives you the lifetime?
This seems to be treated differently by different people.

Thanks a bunch
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Non-symmetric PMF across lipid bilayer

[Gmx QA]

Or maybe perform the initial pull-simulation in the other direction as
well, and combine the results from the two simulations?


Den ons 21 nov. 2018 kl 18:30 skrev Gmx QA :

> Hi Jochen
>
> Thanks for your insights. That might possibly be the case, then? I need to
> visualize my trajectory to find out I guess.
> This is a cyclic peptide being pulled through the membrane, and the
> reaction coordinate is indeed the com of the peptide with respect to the
> membrane center.
>
>  As I said, I (think I) understand the issues with convergence and slow
> degrees of freedom, but I assumed initially that nevertheless my PMF should
> be symmetric (but most likely not fully converged). Could one try to
> restrain the water molecules somehow to keep them out, or is this simply a
> problem of insufficient sampling that cannot be really solved without
> adding a lot more simulation data?
>
> /PK
>
> Den ons 21 nov. 2018 kl 18:22 skrev Jochen Hub :
>
>> Hi PK,
>>
>> I just looked at your PMF profile.xvg - you have a huge offset between
>> the left and right end (140 kJ/mol), suggesting that you have a major
>> problem with convergence. I would not use -cycl here - as such, I
>> disagree here with the previous messages. :)
>>
>> -cycl is acceptable if your PMF looks more or less cyclic (no offset
>> between left and right end). In that case, -cycl makes use of your
>> knowledge that the PMF "should" be cyclic if you would simulate much
>> longer. In other words: Witih -cycl, you get the most plausible PMF in
>> the light of (i) your histograms and (ii) your knowledge that the PMF
>> should be cyclic.
>>
>> I would instead try to find out *why* you have such a massive sampling
>> problem. A reason could be that you have not used a good reaction
>> coordinate (probably the COM of the solute wrt. the membrane?). Maybe
>> you drag a lot of water into the membrane in the xi<0 region, but not in
>> the x>0 region?
>>
>> Good luck,
>> Jochen
>>
>> Am 21.11.18 um 15:49 schrieb Justin Lemkul:
>> >
>> >
>> > On 11/21/18 9:03 AM, Gmx QA wrote:
>> >> Hi Sheryas,
>> >>
>> >> Thanks - so are you saying that the PMF should be asymmetric with
>> respect
>> >> to the leaflets even before I make the PMF symmetric? This seems to
>> >> contradict what Justin just said, but maybe you mean the same thing?
>> >
>> > I think we're saying the same thing, but I must correct an error I
>> made.
>> > I misremembered the WHAM option you need. It is not -sym (which
>> > symmetrizes around zero), you need -cycl, which sets the end points to
>> > be equivalent. This is effectively "symmetry" but with respect to the
>> > end points, as I have been saying, rather than symmetry in the sense of
>> > a "mirror image" at equivalent +/- values of the reaction coordinate.
>> >
>> > -Justin
>> >
>> >> Den ons 21 nov. 2018 kl 14:57 skrev Shreyas Kaptan
>> >> > >>> :
>> >>> I am reasonably sure that if the bilayer composition is asymmetric
>> with
>> >>> respect to leaflets you should see asymmetry in the PMF. In the ideal
>> >>> case,
>> >>> of infinite sampling, you should have zero free energy difference in
>> the
>> >>> bulk solvent for either side.
>> >>>
>> >>> On Wed, Nov 21, 2018 at 2:52 PM Gmx QA 
>> wrote:
>> >>>
>> >>>> Thanks again,
>> >>>>
>> >>>> So then to summarize: Using -sym is appropriate in this case, even
>> >>>> though
>> >>>> the bilayer is asymmetric with respect to lipid composition.  This
>> fact
>> >>>> would show up anyway (in the limit of unlimited sampling?)
>> >>>>
>> >>>>
>> >>>>
>> >>>> Den ons 21 nov. 2018 kl 14:44 skrev Justin Lemkul :
>> >>>>
>> >>>>>
>> >>>>> On 11/21/18 8:39 AM, Per Larsson wrote:
>> >>>>>> Hi,
>> >>>>>>
>> >>>>>> Thanks Justin, but shouldn't the PMF be (more or less) symmetric
>> >>>> anyway,
>> >>>>>> given the inherent bilayer symmetry?
>> >>>>>> In this case I have designed the two leaflets in the bilayer to
>> have
>> >>>>>> non-identical lipid composition, so then I think using -sym 

Re: [gmx-users] Non-symmetric PMF across lipid bilayer

[Gmx QA]

Hi Jochen

Thanks for your insights. That might possibly be the case, then? I need to
visualize my trajectory to find out I guess.
This is a cyclic peptide being pulled through the membrane, and the
reaction coordinate is indeed the com of the peptide with respect to the
membrane center.

 As I said, I (think I) understand the issues with convergence and slow
degrees of freedom, but I assumed initially that nevertheless my PMF should
be symmetric (but most likely not fully converged). Could one try to
restrain the water molecules somehow to keep them out, or is this simply a
problem of insufficient sampling that cannot be really solved without
adding a lot more simulation data?

/PK

Den ons 21 nov. 2018 kl 18:22 skrev Jochen Hub :

> Hi PK,
>
> I just looked at your PMF profile.xvg - you have a huge offset between
> the left and right end (140 kJ/mol), suggesting that you have a major
> problem with convergence. I would not use -cycl here - as such, I
> disagree here with the previous messages. :)
>
> -cycl is acceptable if your PMF looks more or less cyclic (no offset
> between left and right end). In that case, -cycl makes use of your
> knowledge that the PMF "should" be cyclic if you would simulate much
> longer. In other words: Witih -cycl, you get the most plausible PMF in
> the light of (i) your histograms and (ii) your knowledge that the PMF
> should be cyclic.
>
> I would instead try to find out *why* you have such a massive sampling
> problem. A reason could be that you have not used a good reaction
> coordinate (probably the COM of the solute wrt. the membrane?). Maybe
> you drag a lot of water into the membrane in the xi<0 region, but not in
> the x>0 region?
>
> Good luck,
> Jochen
>
> Am 21.11.18 um 15:49 schrieb Justin Lemkul:
> >
> >
> > On 11/21/18 9:03 AM, Gmx QA wrote:
> >> Hi Sheryas,
> >>
> >> Thanks - so are you saying that the PMF should be asymmetric with
> respect
> >> to the leaflets even before I make the PMF symmetric? This seems to
> >> contradict what Justin just said, but maybe you mean the same thing?
> >
> > I think we're saying the same thing, but I must correct an error I made.
> > I misremembered the WHAM option you need. It is not -sym (which
> > symmetrizes around zero), you need -cycl, which sets the end points to
> > be equivalent. This is effectively "symmetry" but with respect to the
> > end points, as I have been saying, rather than symmetry in the sense of
> > a "mirror image" at equivalent +/- values of the reaction coordinate.
> >
> > -Justin
> >
> >> Den ons 21 nov. 2018 kl 14:57 skrev Shreyas Kaptan
> >>  >>> :
> >>> I am reasonably sure that if the bilayer composition is asymmetric with
> >>> respect to leaflets you should see asymmetry in the PMF. In the ideal
> >>> case,
> >>> of infinite sampling, you should have zero free energy difference in
> the
> >>> bulk solvent for either side.
> >>>
> >>> On Wed, Nov 21, 2018 at 2:52 PM Gmx QA  wrote:
> >>>
> >>>> Thanks again,
> >>>>
> >>>> So then to summarize: Using -sym is appropriate in this case, even
> >>>> though
> >>>> the bilayer is asymmetric with respect to lipid composition.  This
> fact
> >>>> would show up anyway (in the limit of unlimited sampling?)
> >>>>
> >>>>
> >>>>
> >>>> Den ons 21 nov. 2018 kl 14:44 skrev Justin Lemkul :
> >>>>
> >>>>>
> >>>>> On 11/21/18 8:39 AM, Per Larsson wrote:
> >>>>>> Hi,
> >>>>>>
> >>>>>> Thanks Justin, but shouldn't the PMF be (more or less) symmetric
> >>>> anyway,
> >>>>>> given the inherent bilayer symmetry?
> >>>>>> In this case I have designed the two leaflets in the bilayer to have
> >>>>>> non-identical lipid composition, so then I think using -sym would
> >>>>>> obliterate any differences between the leaflets, no?
> >>>>> No, because that's completely unknown (and irrelevant) to WHAM. It
> >>>>> sets
> >>>>> the leftmost window to a zero energy and calculates every window's
> >>>>> energy relative to that, so you'll get the steady increase you see.
> If
> >>>>> you tell it that the leftmost and rightmost windows are equal (which
> >>>>> -sym), then the calculation proceeds differently.
> >>>>>
> >>>

Re: [gmx-users] Non-symmetric PMF across lipid bilayer

[Gmx QA]

Thanks again

But thinking about this a bit further - if, as you say, wham sets the
leftmost window to zero and calculates energies for all consecutive windows
relative to that, then shouldn't windows that are equivalent in terms of
the molecules being surrounded only by water, but on opposing sides of the
bilayer, have the same relative energy in the PMFs without any cyclization?
I might be missing the finer details about the wham equations here, but
would appreciate some further insights :-)



Den ons 21 nov. 2018 kl 15:49 skrev Justin Lemkul :

>
>
> On 11/21/18 9:03 AM, Gmx QA wrote:
> > Hi Sheryas,
> >
> > Thanks - so are you saying that the PMF should be asymmetric with respect
> > to the leaflets even before I make the PMF symmetric? This seems to
> > contradict what Justin just said, but maybe you mean the same thing?
>
> I think we're saying the same thing, but I must correct an error I made.
> I misremembered the WHAM option you need. It is not -sym (which
> symmetrizes around zero), you need -cycl, which sets the end points to
> be equivalent. This is effectively "symmetry" but with respect to the
> end points, as I have been saying, rather than symmetry in the sense of
> a "mirror image" at equivalent +/- values of the reaction coordinate.
>
> -Justin
>
> > Den ons 21 nov. 2018 kl 14:57 skrev Shreyas Kaptan <
> shreyaskap...@gmail.com
> >> :
> >> I am reasonably sure that if the bilayer composition is asymmetric with
> >> respect to leaflets you should see asymmetry in the PMF. In the ideal
> case,
> >> of infinite sampling, you should have zero free energy difference in the
> >> bulk solvent for either side.
> >>
> >> On Wed, Nov 21, 2018 at 2:52 PM Gmx QA  wrote:
> >>
> >>> Thanks again,
> >>>
> >>> So then to summarize: Using -sym is appropriate in this case, even
> though
> >>> the bilayer is asymmetric with respect to lipid composition.  This fact
> >>> would show up anyway (in the limit of unlimited sampling?)
> >>>
> >>>
> >>>
> >>> Den ons 21 nov. 2018 kl 14:44 skrev Justin Lemkul :
> >>>
> >>>>
> >>>> On 11/21/18 8:39 AM, Per Larsson wrote:
> >>>>> Hi,
> >>>>>
> >>>>> Thanks Justin, but shouldn't the PMF be (more or less) symmetric
> >>> anyway,
> >>>>> given the inherent bilayer symmetry?
> >>>>> In this case I have designed the two leaflets in the bilayer to have
> >>>>> non-identical lipid composition, so then I think using -sym would
> >>>>> obliterate any differences between the leaflets, no?
> >>>> No, because that's completely unknown (and irrelevant) to WHAM. It
> sets
> >>>> the leftmost window to a zero energy and calculates every window's
> >>>> energy relative to that, so you'll get the steady increase you see. If
> >>>> you tell it that the leftmost and rightmost windows are equal (which
> >>>> -sym), then the calculation proceeds differently.
> >>>>
> >>>> -Justin
> >>>>
> >>>>>
> >>>>>
> >>>>> On Wed, Nov 21, 2018 at 2:24 PM Justin Lemkul 
> >> wrote:
> >>>>>> On 11/21/18 7:22 AM, Gmx QA wrote:
> >>>>>>> Hi all gmx-users
> >>>>>>>
> >>>>>>> I am working on calculating the PMF using umbrella sampling of a
> >>>> (rather
> >>>>>>> large) molecule across a lipid bilayer. I have set up my umbrellas
> >>>> with a
> >>>>>>> 0,2 nm spacing, and run each window for 100 ns.
> >>>>>>>
> >>>>>>> The problem is that the resulting PMF is not symmetric with respect
> >>> to
> >>>>>> the
> >>>>>>> bilayer center. Initially is looks ok, but when the molecule is
> >>> exiting
> >>>>>> the
> >>>>>>> bilayer on the other side again, the PMF does not go back to
> >>> (roughly)
> >>>>>> the
> >>>>>>> same value as before entering the bilayer.
> >>>>>>>
> >>>>>>> I have uploaded the PMF file here:
> >>>>>>> https://files.fm/u/7ec2rshc
> >>>>>> You didn't get a symmetric profile because you didn't ask for one.
> >> Use
> >>>>>> 

Re: [gmx-users] Non-symmetric PMF across lipid bilayer

[Gmx QA]

Hi Sheryas,

Thanks - so are you saying that the PMF should be asymmetric with respect
to the leaflets even before I make the PMF symmetric? This seems to
contradict what Justin just said, but maybe you mean the same thing?

Den ons 21 nov. 2018 kl 14:57 skrev Shreyas Kaptan :

> I am reasonably sure that if the bilayer composition is asymmetric with
> respect to leaflets you should see asymmetry in the PMF. In the ideal case,
> of infinite sampling, you should have zero free energy difference in the
> bulk solvent for either side.
>
> On Wed, Nov 21, 2018 at 2:52 PM Gmx QA  wrote:
>
> > Thanks again,
> >
> > So then to summarize: Using -sym is appropriate in this case, even though
> > the bilayer is asymmetric with respect to lipid composition.  This fact
> > would show up anyway (in the limit of unlimited sampling?)
> >
> >
> >
> > Den ons 21 nov. 2018 kl 14:44 skrev Justin Lemkul :
> >
> > >
> > >
> > > On 11/21/18 8:39 AM, Per Larsson wrote:
> > > > Hi,
> > > >
> > > > Thanks Justin, but shouldn't the PMF be (more or less) symmetric
> > anyway,
> > > > given the inherent bilayer symmetry?
> > > > In this case I have designed the two leaflets in the bilayer to have
> > > > non-identical lipid composition, so then I think using -sym would
> > > > obliterate any differences between the leaflets, no?
> > >
> > > No, because that's completely unknown (and irrelevant) to WHAM. It sets
> > > the leftmost window to a zero energy and calculates every window's
> > > energy relative to that, so you'll get the steady increase you see. If
> > > you tell it that the leftmost and rightmost windows are equal (which
> > > -sym), then the calculation proceeds differently.
> > >
> > > -Justin
> > >
> > > >
> > > >
> > > >
> > > > On Wed, Nov 21, 2018 at 2:24 PM Justin Lemkul 
> wrote:
> > > >
> > > >>
> > > >> On 11/21/18 7:22 AM, Gmx QA wrote:
> > > >>> Hi all gmx-users
> > > >>>
> > > >>> I am working on calculating the PMF using umbrella sampling of a
> > > (rather
> > > >>> large) molecule across a lipid bilayer. I have set up my umbrellas
> > > with a
> > > >>> 0,2 nm spacing, and run each window for 100 ns.
> > > >>>
> > > >>> The problem is that the resulting PMF is not symmetric with respect
> > to
> > > >> the
> > > >>> bilayer center. Initially is looks ok, but when the molecule is
> > exiting
> > > >> the
> > > >>> bilayer on the other side again, the PMF does not go back to
> > (roughly)
> > > >> the
> > > >>> same value as before entering the bilayer.
> > > >>>
> > > >>> I have uploaded the PMF file here:
> > > >>> https://files.fm/u/7ec2rshc
> > > >> You didn't get a symmetric profile because you didn't ask for one.
> Use
> > > >> the -sym option.
> > > >>
> > > >> -Justin
> > > >>
> > > >>> Any comments or suggestions are much appreciated. I understand the
> > > >> problems
> > > >>> and issues about calculations of a converged PMF with larger
> > molecules,
> > > >> but
> > > >>> nevertheless I would have expected my PMF to be symmetric, albeit
> > > perhaps
> > > >>> not converged.
> > > >>>
> > > >>> Thanks
> > > >>> /PK
> > > >> --
> > > >> ==
> > > >>
> > > >> Justin A. Lemkul, Ph.D.
> > > >> Assistant Professor
> > > >> Office: 301 Fralin Hall
> > > >> Lab: 303 Engel Hall
> > > >>
> > > >> Virginia Tech Department of Biochemistry
> > > >> 340 West Campus Dr.
> > > >> Blacksburg, VA 24061
> > > >>
> > > >> jalem...@vt.edu | (540) 231-3129
> > > >> http://www.thelemkullab.com
> > > >>
> > > >> ==
> > > >>
> > > >> --
> > > >> Gromacs Users mailing list
> > > >>
> > > >> * Please search the archive at
> > > >> http://www.gromacs.org/Support/Mailing_Lists/G

[gmx-users] Non-symmetric PMF across lipid bilayer

[Gmx QA]

Hi all gmx-users

I am working on calculating the PMF using umbrella sampling of a (rather
large) molecule across a lipid bilayer. I have set up my umbrellas with a
0,2 nm spacing, and run each window for 100 ns.

The problem is that the resulting PMF is not symmetric with respect to the
bilayer center. Initially is looks ok, but when the molecule is exiting the
bilayer on the other side again, the PMF does not go back to (roughly) the
same value as before entering the bilayer.

I have uploaded the PMF file here:
https://files.fm/u/7ec2rshc

Any comments or suggestions are much appreciated. I understand the problems
and issues about calculations of a converged PMF with larger molecules, but
nevertheless I would have expected my PMF to be symmetric, albeit perhaps
not converged.

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Autocorrelation times from gmx wham and analyze

[Gmx QA]

Hi all,

I am running a series of umbrella sampling simulations of a drug across a
membrane, with the final aim to determine permeability using the local
diffusion model.

I have come across papers where the (necessary) autocorrelation times of
the mean-squared fluctuation of the coordinates are calculated using gmx
wham, and this is also what I have done like this:

gmx wham -if pullf-files.dat -it tpr-files.dat -o -hist -ac -oiact


 and then the ac-times are written to the file iact.xvg.


To understand better the procedure however, I also wanted to do a similar
calculation manually.

I therefore did this:


gmx traj -f ../umbrella0.part0001.xtc -s ../umbrella0.tpr -ox
umbrella0_z_com.xvg -com -nox -noy


 to extract the z-com-coordinates of the drug molecule for a particular us
window.


Then:


gmx analyze -f umbrella0_z_com.xvg -ac -fitfn exp


Followed by:


gmx analyze -f autocorr.xvg  -integrate


to get what I though was going to be the equivalent acf-time tau. However,
for this umbrella window gmx wham gives me tau =  2.25 whereas the second
approach gives a value of 122. What could be the reason for this
difference? There is some info in cmx wham that says that the act
calculation is cut at a value of 0.05, and that is something that I haven't
been able to replicate in my second approach, but any insights are highly
appreciated.


Cheers

/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Gromacs 2018 and GPU PME

[Gmx QA]

Hi Szilard

Thank you for answering. It did indeed show a significant improvement
with, in particular,

$ gmx mdrun -v -deffnm md -pme gpu -nb gpu -ntmpi 8 -ntomp 6 -npme 1
-gputasks 0001

I also now understand better how to control each individual
simulation. Your point on maximizing aggregate performance is well
taken :-)

Thanks again
/PK

2018-02-09 16:25 GMT+01:00, Szilárd Páll <pall.szil...@gmail.com>:
> Hi,
>
> First of all,have you read the docs (admittedly somewhat brief):
> http://manual.gromacs.org/documentation/2018/user-guide/mdrun-performance.html#types-of-gpu-tasks
>
> The current PME GPU was optimized for single-GPU runs. Using multiple GPUs
> with PME offloaded works, but this mode hasn't been an optimization target
> and it will often not give very good performance. Using multiple GPUs
> requires a separate PME rank (as you have realized), only one can be used
> (as we don't support PME decomposition on GPUs) and it comes some inherent
> scaling drawbacks. For this reason, unless you _need_ your single run to be
> as fast as possible, you'll be better off running multiple simulations
> side-by side.
>
> A few tips for tuning the performance of a multi-GPU run with PME offload:
> * expect to get at best 1.5 scaling to 2 GPUs (rarely 3 if the tasks allow)
> * generally it's best to use about the same decomposition that you'd use
> with nonbonded-only offload, e.g. in your case 6-8 ranks
> * map the GPU task alone or at most together with 1 PP rank to a GPU, i.e.
> use the new -gputasks option
> e.g. for your case I'd expect the following to work ~best:
> gmx mdrun -v -deffnm md -pme gpu -nb gpu -ntmpi 8 -ntomp 6 -npme 1
> -gputasks 0001
> or
> gmx mdrun -v -deffnm md -pme gpu -nb gpu -ntmpi 8 -ntomp 6 -npme 1
> -gputasks 0011
>
>
> Let me know if that gave some improvement.
>
> Cheers,
>
> --
> Szilárd
>
> On Fri, Feb 9, 2018 at 8:51 AM, Gmx QA <gmxquesti...@gmail.com> wrote:
>
>> Hi list,
>>
>> I am trying out the new gromacs 2018 (really nice so far), but have a few
>> questions about what command line options I should specify, specifically
>> with the new gnu pme implementation.
>>
>> My computer has two CPUs (with 12 cores each, 24 with hyper threading) and
>> two GPUs, and I currently (with 2018) start simulations like this:
>>
>> $ gmx mdrun -v -deffnm md -pme gpu -nb gpu -ntmpi 2 -npme 1 -ntomp 24
>> -gpu_id 01
>>
>> this works, but gromacs prints the message that 24 omp threads per mpi
>> rank
>> is likely inefficient. However, trying to reduce the number of omp threads
>> I see a reduction in performance. Is this message no longer relevant with
>> gpu pme or am I overlooking something?
>>
>> Thanks
>> /PK
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at http://www.gromacs.org/
>> Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
> mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Gromacs 2018 and GPU PME

[Gmx QA]

Hi list,

I am trying out the new gromacs 2018 (really nice so far), but have a few
questions about what command line options I should specify, specifically
with the new gnu pme implementation.

My computer has two CPUs (with 12 cores each, 24 with hyper threading) and
two GPUs, and I currently (with 2018) start simulations like this:

$ gmx mdrun -v -deffnm md -pme gpu -nb gpu -ntmpi 2 -npme 1 -ntomp 24
-gpu_id 01

this works, but gromacs prints the message that 24 omp threads per mpi rank
is likely inefficient. However, trying to reduce the number of omp threads
I see a reduction in performance. Is this message no longer relevant with
gpu pme or am I overlooking something?

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Cylinder pulling through bilayer

[Gmx QA]

Anyone?

Thanks again
/PK

2017-06-07 21:22 GMT+02:00 Gmx QA <gmxquesti...@gmail.com>:

> Dear list
>
> I am attempting to pull a small molecule though a bilayer using the pull
> geometry cylinder with gromacs v 2016.
>
> This is the relevant portion of my mdp-file:
>
> pull  = yes
> pull-ngroups = 2
> pull-ncoords = 1
> pull-coord1-groups   = 1 2
> pull-group1-name = LIG
> pull-group2-name = MEM
> pull-coord1-type = umbrella
> pull-coord1-geometry = cylinder
> pull-coord1-rate = 0.1
> pull-coord1-vec  = 0 0 1
> pull-coord1-k= 1000
> pull-coord1-start= yes
> pull-coord1-init = 0
> pull-cylinder-r  = 1.5
>
> The pull-rate is very fast because I'm only doing preliminary test. At the
> start, the drug molecule is in -z position compared to the membrane.
>
> When doing grompp:
>
> $ gmx grompp -f umbrella_md_test.mdp -c npt.gro -p topol.top -o
> pull_test.tpr
>
> The output makes no sense:
> Using a fourier grid of 72x72x192, spacing 0.113 0.113 0.111
> Pull group  natoms  pbc atom  distance at start  reference at t=0
>13618
>2 32500 64286-nan nm   -nan nm
> Estimate for the relative computational load of the PME mesh part: 0.44
> This run will generate roughly 14 Mb of data
>
>
> I.e. nan's for distance. If I however switch in the mdp file so that
> pull-group1-name = MEM and pull-group2-name = LIG, the distance gets
> correctly calculated. But this does not seem to be what is prescribed in
> the manual for cylinder pulling, where is says that the cylinder is formed
> from the first group (should be the drug molecule) and through the com of
> the reference group (the membrane in my case),
>
> I think there is something I am missing?
>
> Thanks!
> /PK
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Cylinder pulling through bilayer

[Gmx QA]

Dear list

I am attempting to pull a small molecule though a bilayer using the pull
geometry cylinder with gromacs v 2016.

This is the relevant portion of my mdp-file:

pull  = yes
pull-ngroups = 2
pull-ncoords = 1
pull-coord1-groups   = 1 2
pull-group1-name = LIG
pull-group2-name = MEM
pull-coord1-type = umbrella
pull-coord1-geometry = cylinder
pull-coord1-rate = 0.1
pull-coord1-vec  = 0 0 1
pull-coord1-k= 1000
pull-coord1-start= yes
pull-coord1-init = 0
pull-cylinder-r  = 1.5

The pull-rate is very fast because I'm only doing preliminary test. At the
start, the drug molecule is in -z position compared to the membrane.

When doing grompp:

$ gmx grompp -f umbrella_md_test.mdp -c npt.gro -p topol.top -o
pull_test.tpr

The output makes no sense:
Using a fourier grid of 72x72x192, spacing 0.113 0.113 0.111
Pull group  natoms  pbc atom  distance at start  reference at t=0
   13618
   2 32500 64286-nan nm   -nan nm
Estimate for the relative computational load of the PME mesh part: 0.44
This run will generate roughly 14 Mb of data


I.e. nan's for distance. If I however switch in the mdp file so that
pull-group1-name = MEM and pull-group2-name = LIG, the distance gets
correctly calculated. But this does not seem to be what is prescribed in
the manual for cylinder pulling, where is says that the cylinder is formed
from the first group (should be the drug molecule) and through the com of
the reference group (the membrane in my case),

I think there is something I am missing?

Thanks!
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Gromacs 2016 pull options

[Gmx QA]

Ah, ok, so the rate does not specify the velocity by which the pull-group
moves?

Does it have the same function in the case of constant-force pulling?

2017-02-20 14:20 GMT+01:00 Justin Lemkul <jalem...@vt.edu>:

>
>
> On 2/20/17 8:19 AM, Gmx QA wrote:
>
>> Hi Justin
>>
>> Thank you for your response.
>>
>> I will test that obviously, but just for my understanding then:
>> Is it not enough to specify the pull-coord1-rate, you also need the
>> pull-coord1-k to make it work?
>>
>> They seem to be specifying sort of the same thing really…
>>
>>
> The rate is the velocity at which the biasing spring extends.  The force
> constant calculates how much force is thus applied to the restrained groups
> based on that displacement.  They are independent of one another,
> functionally and mathematically.
>
> -Justin
>
>
> Thanks
>> /PK
>>
>> 2017-02-20 14:13 GMT+01:00 Justin Lemkul <jalem...@vt.edu>:
>>
>>
>>>
>>> On 2/20/17 8:09 AM, Gmx QA wrote:
>>>
>>> Dear list
>>>>
>>>> Could someone please clarify a few issues I have with the 2016 version
>>>> of
>>>> the pull code. Or rather, explain why the below approach is not working.
>>>>
>>>> I am trying to pull a small molecule though a bilayer. Initially, the
>>>> drug
>>>> is "below" the membrane, i e has a smaller z-coordinate, so I want to
>>>> pull
>>>> in the positive z direction through the membrane.
>>>>
>>>> The following is a snippet of my mdp-file, with the relevant parts for
>>>> the
>>>> pull specs.
>>>>
>>>> pull  = yes
>>>> pull-ngroups = 2
>>>> pull-ncoords = 1
>>>> pull-coord1-groups   = 1 2
>>>> pull-group1-name = drug
>>>> pull-group2-name = membrane
>>>> pull-coord1-type = umbrella
>>>> pull-coord1-geometry = direction-periodic
>>>> pull-coord1-vec  = 0 0 1
>>>> pull-coord1-rate = 0.1
>>>>
>>>> This I think should define two pull groups,one for the drug and one for
>>>> the
>>>> membrane, and it should pull the drug molecule with pull-coord1-rate in
>>>> positive z. I am using direction-periodic since I expect to pull more
>>>> than
>>>> half the box z length.
>>>>
>>>> When I test this, running a 1 ns simulation, during which I expect the
>>>> drug
>>>> to be pulled a significant way, nothing really happens with the com of
>>>> the
>>>> drug, it just "sits" out in the water phase.
>>>>
>>>>
>>>> You didn't set a value of pull-coord1-k, so it defaults to zero and
>>> therefore you get no biasing potential.
>>>
>>> http://manual.gromacs.org/documentation/2016.2/user-guide/
>>> mdp-options.html#com-pulling
>>>
>>> -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
>>>
>>> jalem...@outerbanks.umaryland.edu | (410) 706-7441
>>> http://mackerell.umaryland.edu/~jalemkul
>>>
>>> ==
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at http://www.gromacs.org/Support
>>> /Mailing_Lists/GMX-Users_List before posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>> send a mail to gmx-users-requ...@gromacs.org.
>>>
>>>
> --
> ==
>
> 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
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support
> /Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Gromacs 2016 pull options

[Gmx QA]

Hi Justin

Thank you for your response.

I will test that obviously, but just for my understanding then:
Is it not enough to specify the pull-coord1-rate, you also need the
pull-coord1-k to make it work?

They seem to be specifying sort of the same thing really…

Thanks
/PK

2017-02-20 14:13 GMT+01:00 Justin Lemkul <jalem...@vt.edu>:

>
>
> On 2/20/17 8:09 AM, Gmx QA wrote:
>
>> Dear list
>>
>> Could someone please clarify a few issues I have with the 2016 version of
>> the pull code. Or rather, explain why the below approach is not working.
>>
>> I am trying to pull a small molecule though a bilayer. Initially, the drug
>> is "below" the membrane, i e has a smaller z-coordinate, so I want to pull
>> in the positive z direction through the membrane.
>>
>> The following is a snippet of my mdp-file, with the relevant parts for the
>> pull specs.
>>
>> pull  = yes
>> pull-ngroups = 2
>> pull-ncoords = 1
>> pull-coord1-groups   = 1 2
>> pull-group1-name = drug
>> pull-group2-name = membrane
>> pull-coord1-type = umbrella
>> pull-coord1-geometry = direction-periodic
>> pull-coord1-vec  = 0 0 1
>> pull-coord1-rate = 0.1
>>
>> This I think should define two pull groups,one for the drug and one for
>> the
>> membrane, and it should pull the drug molecule with pull-coord1-rate in
>> positive z. I am using direction-periodic since I expect to pull more than
>> half the box z length.
>>
>> When I test this, running a 1 ns simulation, during which I expect the
>> drug
>> to be pulled a significant way, nothing really happens with the com of the
>> drug, it just "sits" out in the water phase.
>>
>>
> You didn't set a value of pull-coord1-k, so it defaults to zero and
> therefore you get no biasing potential.
>
> http://manual.gromacs.org/documentation/2016.2/user-guide/
> mdp-options.html#com-pulling
>
> -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
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support
> /Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Gromacs 2016 pull options

[Gmx QA]

Dear list

Could someone please clarify a few issues I have with the 2016 version of
the pull code. Or rather, explain why the below approach is not working.

I am trying to pull a small molecule though a bilayer. Initially, the drug
is "below" the membrane, i e has a smaller z-coordinate, so I want to pull
in the positive z direction through the membrane.

The following is a snippet of my mdp-file, with the relevant parts for the
pull specs.

pull  = yes
pull-ngroups = 2
pull-ncoords = 1
pull-coord1-groups   = 1 2
pull-group1-name = drug
pull-group2-name = membrane
pull-coord1-type = umbrella
pull-coord1-geometry = direction-periodic
pull-coord1-vec  = 0 0 1
pull-coord1-rate = 0.1

This I think should define two pull groups,one for the drug and one for the
membrane, and it should pull the drug molecule with pull-coord1-rate in
positive z. I am using direction-periodic since I expect to pull more than
half the box z length.

When I test this, running a 1 ns simulation, during which I expect the drug
to be pulled a significant way, nothing really happens with the com of the
drug, it just "sits" out in the water phase.

Please advice!
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Calculating integrals?

[Gmx QA]

Dear users

I have a question that is not explicitly related to Gromacs, but I hope
someone can answer anyway.

More often than not when you read other people's papers, they are showing
one or two equations with integrals. Often this integral runs in some
direction, or over the entire box from some value to some other value.

If I have some data in an array, and I want to integrate this over a range,
how you do that in general? I mean, using for example the trapezoidal rule
will give me the entire area under the curve, but if I am interested in the
value of the integral as a function of (say) distance?

I am not sure this makes sense, but if someone has a script somewhere were
they are doing something simple, and would be willing to share, that would
be much appreciated…

Best
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] martini lipids self-assembly pressure coupling

[Gmx QA]

Hi all

I am running a simulation with martini lipids with the lipids initially
placed randomly together with W solvent in a box. What pressure control do
I apply?

I have tried isotropic, since the system indeed is isotropic initially, but
I don't see formation of a bilayer, rather multiple layers of lipids and/or
cylindrical weird structures. Going anisotropic always deforms the box, it
becomes longer and longer in z until almost infinity. Semi-isotropic is
correct for bilayers obviously, but can I use that before the bilayer is
formed?

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Center on bilayer

[Gmx QA]

I realize this has been discussed at length before, but I can't for the
life of me get it to work.

I have a system where I want to assemble a bilayer from randomly placed
lipids (using charmm36 that is not the main point).

After about 50 ns I have what seems to be two (rough) leaflets being
formed, but of course one leaflet sits at the "top" of the box and the
other at the bottom, with water in between.

I want a trajectory with the bilayer in the center.

Reading some of the previous suggestions has led me to understand that
simply using -center on the entire POPC group will not work, since they in
fact already have their geometrical center in the box.

I then did try and select a single atom from of the tails in an arbitrary
lipid (the C18-atom) and use that as a centering group, but that does not
work either (in the sense that the two leaflets are still separated by a
chunk of water).

Is there any other way of accomplishing this, such as e.g. shifting the
center of the box or so?

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

Thanks again

In my hands at least, -pbc mol -ur compact seems to be necessary, as
without it I get a lot of  instead of actual numbers.

I tried to find a reference for the chemical potential of TIP3P, but
haven't succeeded yet. Related to gromacs, I find one thread in the
archives that claims that for TIP3p as well as SPC/E the agreement with
experiment is reasonable, but with problems for TIP5P.

https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2014-May/089540.html

Also, does it matter if the rerun is performed on a NVT or NPT trajectory?
This I could of course easily test, but it seems that formally at least the
Widow formula is applicable to NVT only?

Anyway, I think I now know enough to proceed. How reliable are TPI results
if I want to measure the chemical potential of a drug molecule? How hard
are they to converge?

Best
/PK


2016-08-16 16:38 GMT+02:00 João M. Damas <jmda...@itqb.unl.pt>:

> I don't recall any necessary post-processing of the trajectory being
> necessary, as it should deal well with PBC conditions.
>
> I really don't know what is happening, but it seems like the fixes are
> going in the right direction.
>
> Try to plot the  along time, to see if it's converging or not.
> Depending on the trend, you may want to extend the simulation time. Other
> reasons for incorrect excess chemical potential may be that the TIP3p model
> may not reproduce well that quantity.
>
> And yeah, it's kJ/mol.
>
> João
>
> On Tue, Aug 16, 2016 at 2:12 PM, Gmx QA <gmxquesti...@gmail.com> wrote:
>
> > Thanks Joao
> >
> > The water box is equilibrated at 298K (I checked density and temperature
> >  with gmx energy, and they show no drift from reasonable averages).
> >
> > However, I am a bit confused as I just tried the approach on the same
> > trajectory, but extended to 2.5 ns and after one pass through trjconv
> with
> > -pbc mol -ur compact options.
> >
> > This seems to work in the sense that I no longer get any infs, but the
> > value of , while having the right magnitude, is still off by about a
> > factor of 2.
> > After 2.5 ns I get  = -45 (I assume the unit is kJ/mol), whereas if i
> > read correctly the chemical potential should be -23.5 or so. Can it be
> that
> > I now need to increase the run time of my water simulation even more?
> >
> > Thanks
> > /PK
> >
> >
> >
> > 2016-08-16 13:47 GMT+02:00 João M. Damas <jmda...@itqb.unl.pt>:
> >
> > > Well, 5 insertions per frame, doesn't look bad. 1 ns should not be
> an
> > > issue.
> > >
> > > Well, actually that value normally starts at inf and should go to the
> > > equilibrium  as you increase sampling. But to be inf after 5
> > > inserts, it does sound weird.
> > >
> > > How did you build your box of water? Is it equilibrated at that
> > temperature
> > > (298 K)?
> > > Please make make that the water molecule to be inserted has the correct
> > > geometrical coordinates and that it is centered at 0,0,0.
> > >
> > > João
> > >
> > > On Tue, Aug 16, 2016 at 12:45 PM, Gmx QA <gmxquesti...@gmail.com>
> wrote:
> > >
> > > > Hi Joao
> > > >
> > > > Thank you for your reply. My mdp-file is pasted below:
> > > >
> > > > integrator  = tpi
> > > > emtol = 1000.0
> > > > emstep  = 0.01
> > > > nsteps= 5
> > > >
> > > > nstlist = 1
> > > > cutoff-scheme   = group
> > > > ns_type = grid
> > > > coulombtype = pme
> > > > rcoulomb= 1.0
> > > > rvdw= 1.0
> > > > pbc = xyz
> > > > nstxtcout = 5000
> > > >
> > > > ; Temperature coupling is on
> > > > tcoupl = V-rescale
> > > > tc-grps  = system
> > > > tau_t  = 0.1
> > > > ref_t= 298
> > > >
> > > > ; Pressure coupling is on
> > > > pcoupl  = no
> > > > pcoupltype  = isotropic
> > > > tau_p   = 2.0
> > > > ref_p   = 1.0
> > > > compressibility = 4.5e-5
> > > >
> > > > ld_seed = -1
> > > >
> > > > I am using the amber99sb forcefield for tip3p, so I believe the 1.0
> nm
> > > > cutoffs to be correct?
> > > >
> > > > I have tried a couple different v

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

Thanks Joao

The water box is equilibrated at 298K (I checked density and temperature
 with gmx energy, and they show no drift from reasonable averages).

However, I am a bit confused as I just tried the approach on the same
trajectory, but extended to 2.5 ns and after one pass through trjconv with
-pbc mol -ur compact options.

This seems to work in the sense that I no longer get any infs, but the
value of , while having the right magnitude, is still off by about a
factor of 2.
After 2.5 ns I get  = -45 (I assume the unit is kJ/mol), whereas if i
read correctly the chemical potential should be -23.5 or so. Can it be that
I now need to increase the run time of my water simulation even more?

Thanks
/PK



2016-08-16 13:47 GMT+02:00 João M. Damas <jmda...@itqb.unl.pt>:

> Well, 5 insertions per frame, doesn't look bad. 1 ns should not be an
> issue.
>
> Well, actually that value normally starts at inf and should go to the
> equilibrium  as you increase sampling. But to be inf after 5
> inserts, it does sound weird.
>
> How did you build your box of water? Is it equilibrated at that temperature
> (298 K)?
> Please make make that the water molecule to be inserted has the correct
> geometrical coordinates and that it is centered at 0,0,0.
>
> João
>
> On Tue, Aug 16, 2016 at 12:45 PM, Gmx QA <gmxquesti...@gmail.com> wrote:
>
> > Hi Joao
> >
> > Thank you for your reply. My mdp-file is pasted below:
> >
> > integrator  = tpi
> > emtol = 1000.0
> > emstep  = 0.01
> > nsteps= 5
> >
> > nstlist = 1
> > cutoff-scheme   = group
> > ns_type = grid
> > coulombtype = pme
> > rcoulomb= 1.0
> > rvdw= 1.0
> > pbc = xyz
> > nstxtcout = 5000
> >
> > ; Temperature coupling is on
> > tcoupl = V-rescale
> > tc-grps  = system
> > tau_t  = 0.1
> > ref_t= 298
> >
> > ; Pressure coupling is on
> > pcoupl  = no
> > pcoupltype  = isotropic
> > tau_p   = 2.0
> > ref_p   = 1.0
> > compressibility = 4.5e-5
> >
> > ld_seed = -1
> >
> > I am using the amber99sb forcefield for tip3p, so I believe the 1.0 nm
> > cutoffs to be correct?
> >
> > I have tried a couple different values for the number of insertions, and
> > for some frames different values sometimes give non-inf results, but for
> > the most part it does not seem to make a difference.
> >
> > 1 ns is short, I agree, but since the value of  (which I take to be
> the
> > excess chemical potential I am after) immediately goes to inf when the
> > first mu=inf appears, extending does not make sense to me? Is the
> chemical
> > potential reported elsewhere as well?
> >
> > Thanks
> > /PK
> >
> > 2016-08-16 12:35 GMT+02:00 João M. Damas <jmda...@itqb.unl.pt>:
> >
> > > Yes, all water atoms can't be with coordinates 0,0,0. I don't know how
> > did
> > > that even work... Weird.
> > >
> > > 1 ns is very short time scale. Also, how many insertions are you trying
> > per
> > > frame? That could also help improve the sampling.
> > >
> > > The inf values are normal as it just means that the water is too
> > structured
> > > and it's hard to insert a new water from the sampling point of view.
> > >
> > > Can you post the .mdp you used for the tpi?
> > >
> > > João
> > >
> > > On Tue, Aug 16, 2016 at 9:52 AM, Gmx QA <gmxquesti...@gmail.com>
> wrote:
> > >
> > > > Please - anyone?
> > > >
> > > > I gathered from the mail-list that maybe the entire molecule (water
> in
> > > this
> > > > case) should be centered on origo (rather than having all atoms at
> > > (0,0,0),
> > > > but this gives me only a lot of inf for the chemical potential.
> > > >
> > > > How can I calculate the excess chemical potential of water using TPI
> in
> > > > gromacs?
> > > >
> > > >
> > > >
> > > > 2016-08-15 13:31 GMT+02:00 Gmx QA <gmxquesti...@gmail.com>:
> > > >
> > > > > Dear list
> > > > >
> > > > > I am trying to teach myself how the test particle method for excess
> > > > > chemical potential calculations work. To this end, I created a
> small
> > > > system
> > > 

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

Hi Joao

Thank you for your reply. My mdp-file is pasted below:

integrator  = tpi
emtol = 1000.0
emstep  = 0.01
nsteps= 5

nstlist = 1
cutoff-scheme   = group
ns_type = grid
coulombtype = pme
rcoulomb= 1.0
rvdw= 1.0
pbc = xyz
nstxtcout = 5000

; Temperature coupling is on
tcoupl = V-rescale
tc-grps  = system
tau_t  = 0.1
ref_t= 298

; Pressure coupling is on
pcoupl  = no
pcoupltype  = isotropic
tau_p   = 2.0
ref_p   = 1.0
compressibility = 4.5e-5

ld_seed = -1

I am using the amber99sb forcefield for tip3p, so I believe the 1.0 nm
cutoffs to be correct?

I have tried a couple different values for the number of insertions, and
for some frames different values sometimes give non-inf results, but for
the most part it does not seem to make a difference.

1 ns is short, I agree, but since the value of  (which I take to be the
excess chemical potential I am after) immediately goes to inf when the
first mu=inf appears, extending does not make sense to me? Is the chemical
potential reported elsewhere as well?

Thanks
/PK

2016-08-16 12:35 GMT+02:00 João M. Damas <jmda...@itqb.unl.pt>:

> Yes, all water atoms can't be with coordinates 0,0,0. I don't know how did
> that even work... Weird.
>
> 1 ns is very short time scale. Also, how many insertions are you trying per
> frame? That could also help improve the sampling.
>
> The inf values are normal as it just means that the water is too structured
> and it's hard to insert a new water from the sampling point of view.
>
> Can you post the .mdp you used for the tpi?
>
> João
>
> On Tue, Aug 16, 2016 at 9:52 AM, Gmx QA <gmxquesti...@gmail.com> wrote:
>
> > Please - anyone?
> >
> > I gathered from the mail-list that maybe the entire molecule (water in
> this
> > case) should be centered on origo (rather than having all atoms at
> (0,0,0),
> > but this gives me only a lot of inf for the chemical potential.
> >
> > How can I calculate the excess chemical potential of water using TPI in
> > gromacs?
> >
> >
> >
> > 2016-08-15 13:31 GMT+02:00 Gmx QA <gmxquesti...@gmail.com>:
> >
> > > Dear list
> > >
> > > I am trying to teach myself how the test particle method for excess
> > > chemical potential calculations work. To this end, I created a small
> > system
> > > of about 900 water molecules (TIP3p), and simulated for 1 ns.
> > >
> > > I then set up files for inserting an extra water molecule to calculate
> > mu,
> > > the chemical potential.
> > >
> > > The end of my tpi_start.gro looks like this:
> > > 
> > >   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
> > >   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
> > >   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
> > >   885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
> > >   885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
> > >   885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
> > >3.01188   3.01188   3.01188
> > >
> > > So I think this last water molecule is the one to be inserted. The
> > > topology was also updated accordingly.
> > >
> > > I then made a rerun like this:
> > >
> > > $ gmx mdrun -v -deffnm tpi -rerun md.xtc
> > >
> > > And the output is a series of lines like this:
> > >
> > > Reading frame 180 time  900.000   mu  8.924e+00   7.240e+00
> > > mu  7.671e+00   7.242e+00
> > > mu  6.987e+00   7.241e+00
> > > mu  7.010e+00   7.239e+00
> > > mu  7.691e+00   7.241e+00
> > > mu  7.439e+00   7.242e+00
> > > mu  6.757e+00   7.240e+00
> > > mu  8.114e+00   7.243e+00
> > > mu  7.173e+00   7.243e+00
> > > mu  8.768e+00   7.249e+00
> > > Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
> > > mu  9.284e+00   7.259e+00
> > > mu  8.103e+00   7.262e+00
> > > mu  7.835e+00   7.265e+00
> > > mu  7.321e+00   7.265e+00
> > > mu  7.576e+00   7.267e+00
> > > mu  9.348e+00   7.274e+00
> > > mu  7.021e+00   7.272e+00
> > > mu  7.162e+00   7.272e+00
> > > mu  7.695e+00   7.274e+00
> > > Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00
> > >
> > > Now, is  here the average computed chemical potential? What is the
> > TPI
&g

Re: [gmx-users] TPI and chemical potential

[Gmx QA]

Please - anyone?

I gathered from the mail-list that maybe the entire molecule (water in this
case) should be centered on origo (rather than having all atoms at (0,0,0),
but this gives me only a lot of inf for the chemical potential.

How can I calculate the excess chemical potential of water using TPI in
gromacs?



2016-08-15 13:31 GMT+02:00 Gmx QA <gmxquesti...@gmail.com>:

> Dear list
>
> I am trying to teach myself how the test particle method for excess
> chemical potential calculations work. To this end, I created a small system
> of about 900 water molecules (TIP3p), and simulated for 1 ns.
>
> I then set up files for inserting an extra water molecule to calculate mu,
> the chemical potential.
>
> The end of my tpi_start.gro looks like this:
> 
>   884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
>   884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
>   884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
>   885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
>   885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
>   885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
>3.01188   3.01188   3.01188
>
> So I think this last water molecule is the one to be inserted. The
> topology was also updated accordingly.
>
> I then made a rerun like this:
>
> $ gmx mdrun -v -deffnm tpi -rerun md.xtc
>
> And the output is a series of lines like this:
>
> Reading frame 180 time  900.000   mu  8.924e+00   7.240e+00
> mu  7.671e+00   7.242e+00
> mu  6.987e+00   7.241e+00
> mu  7.010e+00   7.239e+00
> mu  7.691e+00   7.241e+00
> mu  7.439e+00   7.242e+00
> mu  6.757e+00   7.240e+00
> mu  8.114e+00   7.243e+00
> mu  7.173e+00   7.243e+00
> mu  8.768e+00   7.249e+00
> Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
> mu  9.284e+00   7.259e+00
> mu  8.103e+00   7.262e+00
> mu  7.835e+00   7.265e+00
> mu  7.321e+00   7.265e+00
> mu  7.576e+00   7.267e+00
> mu  9.348e+00   7.274e+00
> mu  7.021e+00   7.272e+00
> mu  7.162e+00   7.272e+00
> mu  7.695e+00   7.274e+00
> Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00
>
> Now, is  here the average computed chemical potential? What is the TPI
> energy distribution being reported in the tpi.xvg-file?
>
> I tried to google for the chemical potential of water and according to
> e.g. [1] (page 13) it should be around -23.5 kJ/mol. I would really
> appreciate if someone could comment on this, as I need to understand it
> before moving on to more relevant (for me) systems….
>
> Thanks
> /PK
>
>
>
>
>
> [1] Chemical potential of liquids and mixtures via Adaptive Resolution ...
> <https://www.google.se/url?sa=t=j==s=web=21=0ahUKEwj-7-2gqMPOAhVFWiwKHdp4DYY4FBAWCBwwAA=https%3A%2F%2Fopus4.kobv.de%2Fopus4-zib%2Ffiles%2F5097%2FZIB-Report_14-25.pdf=AFQjCNECRk9dif8TGxKJeeztHV6T_FmVuw=-vsjUB9HRdlcnt4vyKWcbw>
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] TPI and chemical potential

[Gmx QA]

Dear list

I am trying to teach myself how the test particle method for excess
chemical potential calculations work. To this end, I created a small system
of about 900 water molecules (TIP3p), and simulated for 1 ns.

I then set up files for inserting an extra water molecule to calculate mu,
the chemical potential.

The end of my tpi_start.gro looks like this:

  884SOL OW 2650   2.466   0.788   2.747 -0.5599 -0.0387  0.0570
  884SOLHW1 2651   2.529   0.835   2.802  0.0024 -0.9901  0.2261
  884SOLHW2 2652   2.506   0.703   2.732 -1.4472 -0.3598 -0.5022
  885SOL OW 2653   0.000   0.000   0.000  0.  0.  0.
  885SOLHW1 2654   0.000   0.000   0.000  0.  0.  0.
  885SOLHW2 2655   0.000   0.000   0.000  0.  0.  0.
   3.01188   3.01188   3.01188

So I think this last water molecule is the one to be inserted. The topology
was also updated accordingly.

I then made a rerun like this:

$ gmx mdrun -v -deffnm tpi -rerun md.xtc

And the output is a series of lines like this:

Reading frame 180 time  900.000   mu  8.924e+00   7.240e+00
mu  7.671e+00   7.242e+00
mu  6.987e+00   7.241e+00
mu  7.010e+00   7.239e+00
mu  7.691e+00   7.241e+00
mu  7.439e+00   7.242e+00
mu  6.757e+00   7.240e+00
mu  8.114e+00   7.243e+00
mu  7.173e+00   7.243e+00
mu  8.768e+00   7.249e+00
Reading frame 190 time  950.000   mu  7.799e+00   7.252e+00
mu  9.284e+00   7.259e+00
mu  8.103e+00   7.262e+00
mu  7.835e+00   7.265e+00
mu  7.321e+00   7.265e+00
mu  7.576e+00   7.267e+00
mu  9.348e+00   7.274e+00
mu  7.021e+00   7.272e+00
mu  7.162e+00   7.272e+00
mu  7.695e+00   7.274e+00
Reading frame 200 time 1000.000   mu  7.104e+00   7.273e+00

Now, is  here the average computed chemical potential? What is the TPI
energy distribution being reported in the tpi.xvg-file?

I tried to google for the chemical potential of water and according to e.g.
[1] (page 13) it should be around -23.5 kJ/mol. I would really appreciate
if someone could comment on this, as I need to understand it before moving
on to more relevant (for me) systems….

Thanks
/PK





[1] Chemical potential of liquids and mixtures via Adaptive Resolution ...

-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Reproducing Charmm36 POPC electron density

[Gmx QA]

Hi Justin,

Thanks, that's reassuring then. My reason for worrying was in part because
the reported electron densities in several papers I have read seem to match
very closely. But since as you say the APL values also agree then
everything might be ok after all?

Thanks
/PK


2016-02-24 15:52 GMT+01:00 Justin Lemkul <jalem...@vt.edu>:

>
>
> On 2/24/16 7:25 AM, Gmx QA wrote:
>
>> Dear list
>>
>> I have a simulation of a system of 64 POPC lipids per leaflet using the
>> Charmm36 lipid parameter set, and I am calculating the electron density
>> across the bilayer to compare to literature values. Unfortunately I cannot
>> quite get the same results as eg. in Klauda et al (
>> http://www.sciencedirect.com/science/article/pii/S0005273614002193)
>>
>> I have made the starting structure using the charmm-gui, and gone through
>> all the equilibration steps as suggested in the gromacs tar-ball from that
>> website.
>>
>> My production mdp file is the following:
>>
>> integrator  = md
>> dt  = 0.002
>> nsteps  = 5000
>> nstlog  = 1
>> nstxout = 50
>> nstvout = 50
>> nstfout = 0
>> nstcalcenergy   = 10
>> nstenergy   = 1000
>> nstxtcout   = 50
>> ;
>> 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  = Nose-Hoover
>> tc_grps = POPC SOL
>> tau_t   = 1.0 1.0
>> ref_t   = 303.15  303.15
>> ;nsttcouple = 1
>> ;
>> pcoupl  = Parrinello-Rahman
>> pcoupltype  = semiisotropic
>> tau_p   = 5.0 5.0
>> compressibility = 4.5e-5 4.5e-5
>> ref_p   = 1.0 1.0
>> ;
>> constraints = h-bonds
>> constraint_algorithm= LINCS
>> continuation= yes
>> ;
>> nstcomm = 10
>> comm_mode   = linear
>> comm_grps   = POPC SOL
>> ;
>> refcoord_scaling= com
>>
>> Simulations are run with gromacs 5.0.4, and currently 100 ns long.
>>
>> Before running vmx density, I center my trajectory like this:
>> $ gmx trjconv -f run.trr -s em.tpr -o run_center.trr -center
>>
>> selecting POPC as the group to center on.
>>
>>   Then, my command for computing the electron density is the following
>> (also
>> using gmx density 5.0.4):
>>
>> $ gmx density -f run_center.trr -s em.tpr -o  density_traj_center.xvg
>> -dens
>> electron -ei electrons.dat -center
>>
>> Again, selecting POPC to center the density profile on.
>>
>> The resulting profile can be found here (i hope that works)
>> http://www.filedropper.com/compareelectrondensity
>>
>> As you can see, the depth of the central trough is not as deep as the
>> literature value, nor is the height of the density the same even though
>> the
>> peak-peak distance appears to be similar, and values in the other regions
>> also match _fairly_ closely.
>>
>> The average APL value over the entire trajectory is 64.0 +- 1.33, which is
>> in line with the reported 64.7 from Klauda, but the electron density
>> worries me a bit. Perhaps anyone (Justin?) can shed any light on this?
>>
>>
> Why?  The results match almost perfectly.  The magnitude of error is
> within a few % and if your APL matches then I'm sure everything is fine.
> We rigorously tested the lipid inputs, cutoff setup, etc. and everything
> agrees quite well between lipid simulations in CHARMM, GROMACS, AMBER,
> NAMD, and OpenMM.  There are minor variations, as would be expected in
> comparing any two software programs, but nothing is really aberrant.
>
> http://dx.doi.org/10.1021/acs.jctc.5b00935
>
> -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
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul

[gmx-users] Reproducing Charmm36 POPC electron density

[Gmx QA]

Dear list

I have a simulation of a system of 64 POPC lipids per leaflet using the
Charmm36 lipid parameter set, and I am calculating the electron density
across the bilayer to compare to literature values. Unfortunately I cannot
quite get the same results as eg. in Klauda et al (
http://www.sciencedirect.com/science/article/pii/S0005273614002193)

I have made the starting structure using the charmm-gui, and gone through
all the equilibration steps as suggested in the gromacs tar-ball from that
website.

My production mdp file is the following:

integrator  = md
dt  = 0.002
nsteps  = 5000
nstlog  = 1
nstxout = 50
nstvout = 50
nstfout = 0
nstcalcenergy   = 10
nstenergy   = 1000
nstxtcout   = 50
;
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  = Nose-Hoover
tc_grps = POPC SOL
tau_t   = 1.0 1.0
ref_t   = 303.15  303.15
;nsttcouple = 1
;
pcoupl  = Parrinello-Rahman
pcoupltype  = semiisotropic
tau_p   = 5.0 5.0
compressibility = 4.5e-5 4.5e-5
ref_p   = 1.0 1.0
;
constraints = h-bonds
constraint_algorithm= LINCS
continuation= yes
;
nstcomm = 10
comm_mode   = linear
comm_grps   = POPC SOL
;
refcoord_scaling= com

Simulations are run with gromacs 5.0.4, and currently 100 ns long.

Before running vmx density, I center my trajectory like this:
$ gmx trjconv -f run.trr -s em.tpr -o run_center.trr -center

selecting POPC as the group to center on.

 Then, my command for computing the electron density is the following (also
using gmx density 5.0.4):

$ gmx density -f run_center.trr -s em.tpr -o  density_traj_center.xvg -dens
electron -ei electrons.dat -center

Again, selecting POPC to center the density profile on.

The resulting profile can be found here (i hope that works)
http://www.filedropper.com/compareelectrondensity

As you can see, the depth of the central trough is not as deep as the
literature value, nor is the height of the density the same even though the
peak-peak distance appears to be similar, and values in the other regions
also match _fairly_ closely.

The average APL value over the entire trajectory is 64.0 +- 1.33, which is
in line with the reported 64.7 from Klauda, but the electron density
worries me a bit. Perhaps anyone (Justin?) can shed any light on this?

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] viscosity using green-kubo

[Gmx QA]

Dear gmx-users

I sent mail a couple days ago about viscosity calculations, but no-one
answered so maybe it was too lengthy, so better luck this time I hope.

I want to calculate viscosity using Green-Kubo formalism, and have run a
box with tip3p waters for 100 ps in NVE, saving data every 1 fs.

For green-kubo, I then first use gmx energy to get (one of) the
off-diagonal elements in the pressure tensor (Pxy in this case)

$ gmx energy -f nve_prod.edr -o nve_prod_pressure_tensor_pxy.xvg

Then, I use gmx analyze to calculate the auto-correlation function of the
Pxy pressure component

$ gmx analyze -f nve_prod_pressure_tensor_pxy.xvg -ac
nve_prod_pressure_tensor_pxy_acf.xvg -nonormalize

And then again gmx analyze to integrate it.
$ gmx analyze -f nve_prod_pressure_tensor_pxy_acf.xvg -integrate

That gives me :
Calculating the integral using the trapezium rule
Integral 1  1940.24341  +/-0.0

Which I then plug into the green-kubo formula
eta = V/(kB*T)* int(…) = …=

For my V=1.59e-26 m^3, T=303 K and with kB=1.38e-23 I get eta = 0.0074

Now, I am not sure about units, or about the validity of this answer. I
tried changing the unit of the pressure from bar to pascal by multiplying
with 10, but the results makes equally less sense.

I have also read about averaging over several of the off-diagonal pressure
tensor components, which I could try, but is there anything fundamental
that I am missing above?

Hope for better luck with some answers this time!
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Viscosity of water from gmx tcaf

[Gmx QA]

No one that can comment on this? The very least on what data should be
fitted in visco_k.xvg to get the viscosity.

/PK

2015-12-15 10:14 GMT+01:00 Gmx QA <gmxquesti...@gmail.com>:

> Dear users
>
> I am trying to teach myself how to calculate the viscosity of different
> simulations from gmx tcaf, and to that end I made a test system with a box
> of about 4000 water molecules, so that I can see if I can get a viscosity
> that matches that of water.
>
> I have read some mails on this list, and in particular the Berk Hess paper
> about various methods to calculate viscosity in gromacs, but I have some
> additional questions.
>
> From the help-text for gmx tcaf:
> The k-dependent viscosities in the -ov file should be fitted to eta(k) =
> eta_0
> (1 - a k^2) to obtain the viscosity at infinite wavelength.
>
> This is my command line for running gmx tcaf:
> $ gmx tcaf -f run.trr -s run.tpr -ov
>
> which gives me the visc_k.xvg file.
>
> I then tried to fit the data (per above) using xmgrace to y = a0 * (1 - a1
> * y * y), and get a0 = 0.41*10^-3 kg /m*s, which after unit conversion
> seems to be too low? (wikipedia says viscosity of water is 0.001 kg /m *s.
>
> On the other hand, gmx tcaf also says:
>
> When the box is cubic, one can use the option -oc, which averages the
> TCAFs over all k-vectors with the same length. This results in more
> accurate TCAFs. Both the cubic TCAFs and fits are written to -oc The cubic
> eta estimates are also written to -ov.
>
> Since I have a cubic box, I run
> $ gmx tcaf -f run.trr -s run.tpr -ov -oc
>
> which gives me another visc_k,xvg file with some (fitted?) data in it.
> What is this data?
> The following is also printed to stdout
>
> Averaged over k-vectors:
> k  1.273  tau  1.860  Omega  1.790  eta  0.14901 10^-3 kg/(m s)
> k  1.801  tau  1.133  Omega  2.061  eta  0.14087 10^-3 kg/(m s)
> k  2.206  tau  0.420  Omega  0.799  eta  0.09827 10^-3 kg/(m s)
> k  2.547  tau  0.340  Omega  0.599  eta  0.06810 10^-3 kg/(m s)
>
> Does gromacs somehow already compute the viscosity at "infinite
> wavelength"?, or am i missing something.
>
> I'd really, really appreciate some help here.
>
> Best
> /PK
>
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Viscosity of water from gmx tcaf

[Gmx QA]

Dear users

I am trying to teach myself how to calculate the viscosity of different
simulations from gmx tcaf, and to that end I made a test system with a box
of about 4000 water molecules, so that I can see if I can get a viscosity
that matches that of water.

I have read some mails on this list, and in particular the Berk Hess paper
about various methods to calculate viscosity in gromacs, but I have some
additional questions.

>From the help-text for gmx tcaf:
The k-dependent viscosities in the -ov file should be fitted to eta(k) =
eta_0
(1 - a k^2) to obtain the viscosity at infinite wavelength.

This is my command line for running gmx tcaf:
$ gmx tcaf -f run.trr -s run.tpr -ov

which gives me the visc_k.xvg file.

I then tried to fit the data (per above) using xmgrace to y = a0 * (1 - a1
* y * y), and get a0 = 0.41*10^-3 kg /m*s, which after unit conversion
seems to be too low? (wikipedia says viscosity of water is 0.001 kg /m *s.

On the other hand, gmx tcaf also says:

When the box is cubic, one can use the option -oc, which averages the TCAFs
over all k-vectors with the same length. This results in more accurate
TCAFs. Both the cubic TCAFs and fits are written to -oc The cubic eta
estimates are also written to -ov.

Since I have a cubic box, I run
$ gmx tcaf -f run.trr -s run.tpr -ov -oc

which gives me another visc_k,xvg file with some (fitted?) data in it. What
is this data?
The following is also printed to stdout

Averaged over k-vectors:
k  1.273  tau  1.860  Omega  1.790  eta  0.14901 10^-3 kg/(m s)
k  1.801  tau  1.133  Omega  2.061  eta  0.14087 10^-3 kg/(m s)
k  2.206  tau  0.420  Omega  0.799  eta  0.09827 10^-3 kg/(m s)
k  2.547  tau  0.340  Omega  0.599  eta  0.06810 10^-3 kg/(m s)

Does gromacs somehow already compute the viscosity at "infinite
wavelength"?, or am i missing something.

I'd really, really appreciate some help here.

Best
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Calculation of [constrainttypes] lengths

[Gmx QA]

Hello everyone

I know how to calculate the lengths in the [constrainttypes] section has
been brought up before on this list, but the answer from Erik Marklund was
to use bond lengths, angles and moment of inertia.

Well, I tried, but I don't quite arrive at the same answer as in the
forcefield-files, so I'd appreciate if someone could take the time to go
though my math. My aim is to implement some additional vsites, and I'd like
to be as thorough as possible.

I have tried to re-calculate the this from the oplsa.ff/ffbonded.itp -file

; constraints for 2nd type rigid CH3 (angle *-CT-HC is 110.7)
 MCH3BCT 20.167031

My short python-script is as follows:

#!/usr/bin/python


import math

# opls equilibrium values

b0CH = 0.109
b0CC = 0.1522
angle = 110.7

# opls masses
mH = 1.008
mC = 12.011

# mass of the virtual site is then:
mvs = (3*mH+mC)/2

# Start to compute the moment of inertia of the CH3 group when rotating
around
# the "next" carbon atom
# Use the cosine rule to calculate the distance between each hydrogen and
the
# carbon at the center of rotation.
x2 = (b0CC * b0CC) + (b0CH * b0CH) - (2 * b0CC * b0CH *
math.cos(math.radians(angle)))

b = math.sqrt(x2)

# The final moment of inertia for three hydrogens and one carbon is thus
I = (b*b*mH)*3 + b0CC*b0CC*mC

# When using vsites, the moment of inertia must be the same, so use that to
# compute the new constraint length for two mass centers
c = math.sqrt(I/(2*mvs))
print c

When I run this script, I get c = 0.167072934788 whereas the value above is
0.167031. Ok, the difference is only in the third-fourth decimal place, but
since this is based on equilibrium values I'd have expected a better
agreement, perhaps?

Best
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] nhbdist option in gmx hbond?

[Gmx QA]

Anyone?

Thanks
/PK

2015-11-03 10:16 GMT+01:00 Gmx QA <gmxquesti...@gmail.com>:

> Dear list
>
> I am interested in the (new?) option nhbdist in gmx hbond. The help text
> says that it can be used to compare the results to Raman Spectroscopy, but
> after some time with Google I can't find any other documentation about what
> is being compared and how the comparison should be done.
>
> Could someone please point me in the right direction?
>
> Also, is it possible to calculate a Raman spectra with gromacs from an
> existing trajectory?
>
> Thanks!
> /PK
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] nhbdist option in gmx hbond?

[Gmx QA]

Dear list

I am interested in the (new?) option nhbdist in gmx hbond. The help text
says that it can be used to compare the results to Raman Spectroscopy, but
after some time with Google I can't find any other documentation about what
is being compared and how the comparison should be done.

Could someone please point me in the right direction?

Also, is it possible to calculate a Raman spectra with gromacs from an
existing trajectory?

Thanks!
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Differences between gmx bar and pymbar free energy

[Gmx QA]

Hi gmx-users

I am running a series of free-energy of solvation calculations for small
molecules in water and octanol, to calculate partition coefficients.

For this, I have set up a series of lambda-points in water and in octanol,
and then haves used primarily gmx bar to evaluate the results.

Using gmx bar (v.5.0.4) with the following command line
$ gmx bar -f md*xvg -o -oi -oh -b 100 &>FEP.dat

I get for eg. ethanol in water a Delta-G of 14.33 kJ/mol (which is very
close to the value (14.42 kj/mol) reported by eg. Lundborg and Lindahl in
JPCB, 119, 2015, for the same combination of water model and forcefield, so
that is reassuring).

Then I tried the same analysis using pymbar (cloned today from git) like
this:
$ python alchemical_analysis.py -t 300 -p 1 -q xvg -p md -m BAR -v -s 100
&> FEP_pymbar.dat

Now, I get a Delta-G of 64.235 kJ/mol, so clearly not the same as with gmx
bar (and not simply a kJ vs kcal issue either I think)

I feel like I am missing something fundamental since I get so different
values, and would appreciate any feedback. The simulations have been run
using gmx 5.0.4 with a total of 47 lambda points, and each run is 10 ns in
length.

I can provide mdp-files and such if needed.

Cheers
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Differences between gmx bar and pymbar free energy

[Gmx QA]

Hello again

Scratch that, I had the input files numbered incorrectly for pymbar. Now
all is well.

Sorry for the confusion
/PK

2015-10-08 13:22 GMT+02:00 Gmx QA <gmxquesti...@gmail.com>:

> Hi gmx-users
>
> I am running a series of free-energy of solvation calculations for small
> molecules in water and octanol, to calculate partition coefficients.
>
> For this, I have set up a series of lambda-points in water and in octanol,
> and then haves used primarily gmx bar to evaluate the results.
>
> Using gmx bar (v.5.0.4) with the following command line
> $ gmx bar -f md*xvg -o -oi -oh -b 100 &>FEP.dat
>
> I get for eg. ethanol in water a Delta-G of 14.33 kJ/mol (which is very
> close to the value (14.42 kj/mol) reported by eg. Lundborg and Lindahl in
> JPCB, 119, 2015, for the same combination of water model and forcefield, so
> that is reassuring).
>
> Then I tried the same analysis using pymbar (cloned today from git) like
> this:
> $ python alchemical_analysis.py -t 300 -p 1 -q xvg -p md -m BAR -v -s 100
> &> FEP_pymbar.dat
>
> Now, I get a Delta-G of 64.235 kJ/mol, so clearly not the same as with gmx
> bar (and not simply a kJ vs kcal issue either I think)
>
> I feel like I am missing something fundamental since I get so different
> values, and would appreciate any feedback. The simulations have been run
> using gmx 5.0.4 with a total of 47 lambda points, and each run is 10 ns in
> length.
>
> I can provide mdp-files and such if needed.
>
> Cheers
> /PK
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Free energy in octanol and water

[Gmx QA]

Hi Justin

Thanks, that is a good suggestion, I will try that.
Besides running multiple replicates, is there something else that is
obviously wrong with my method?

Best
/PK




2015-09-29 19:09 GMT+02:00 Justin Lemkul <jalem...@vt.edu>:

>
>
> On 9/29/15 4:51 AM, Gmx QA wrote:
>
>> Hi gmx-users
>>
>> I am trying to calculate the free energy of solvation for a drug-molecule
>> in octanol and water.
>>
>> I am comparing two situations:
>>
>> A: Two systems, one with the drug in water only, and one with the drug in
>> pure octanol.
>>
>> B: drug in water only (same as in A), and then one box where I have added
>> 20% water to my octanol system (since mixtures of octanol and water
>> typically contains some water in the octanol phase).
>>
>> Then I run a series of free energy calculations per Justin's tutorial, 48
>> lambda points in total.
>> For system A I get a Delta-G for pure water of 39.07 kJ/mol, and 92.81
>> kJ/mol in pure octanol.
>>
>> For system B, the Delta-G for pure water is the same, but for the
>> octanol-box with 20% water the Delta-G is 95.99 kJ/mol.
>>
>> All of these results are obtained using gmx bar in gromacs 5.0.4.
>>
>> My main concern is that the Delta-G is higher in the water/octanol mixture
>> than in pure octanol, but I would have expected the reverse to be the
>> case?
>>
>>
> From experience, I can tell you it is very hard to get good convergence on
> such mixtures.  You'll note that in npt.gro you effectively have two bands
> of water formed in the box, neither or which is particularly close to the
> API molecule. You'll probably get some microscopic phase separation.  The
> route I've gone in the past is to generate several starting configurations
> (different random seeds when inserting the water molecules) and averaging
> over multiple runs.  You should also do some error analysis of your results
> (both block averaging for convergence, and to compare the two DG values
> between the systems).  Without error bars, you have no idea if those two
> values are at all meaningfully different.
>
> -Justin
>
> Is there something wrong with my protocol? The forcefield is
>> GAFF/amber99sb,and all production runs are 10 ns in length.
>>
>> Best
>> /PK
>>
>> I have uploaded copies of all input-files here:
>>
>> https://dl.dropboxusercontent.com/u/66785227/md.mdp
>> https://dl.dropboxusercontent.com/u/66785227/npt.gro
>> https://dl.dropboxusercontent.com/u/66785227/octanol_GMX.itp
>> https://dl.dropboxusercontent.com/u/66785227/API_gmx.itp
>> https://dl.dropboxusercontent.com/u/66785227/topol.top
>>
>>
> --
> ==
>
> 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
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Free energy in octanol and water

[Gmx QA]

Hi gmx-users

I am trying to calculate the free energy of solvation for a drug-molecule
in octanol and water.

I am comparing two situations:

A: Two systems, one with the drug in water only, and one with the drug in
pure octanol.

B: drug in water only (same as in A), and then one box where I have added
20% water to my octanol system (since mixtures of octanol and water
typically contains some water in the octanol phase).

Then I run a series of free energy calculations per Justin's tutorial, 48
lambda points in total.
For system A I get a Delta-G for pure water of 39.07 kJ/mol, and 92.81
kJ/mol in pure octanol.

For system B, the Delta-G for pure water is the same, but for the
octanol-box with 20% water the Delta-G is 95.99 kJ/mol.

All of these results are obtained using gmx bar in gromacs 5.0.4.

My main concern is that the Delta-G is higher in the water/octanol mixture
than in pure octanol, but I would have expected the reverse to be the case?

Is there something wrong with my protocol? The forcefield is
GAFF/amber99sb,and all production runs are 10 ns in length.

Best
/PK

I have uploaded copies of all input-files here:

https://dl.dropboxusercontent.com/u/66785227/md.mdp
https://dl.dropboxusercontent.com/u/66785227/npt.gro
https://dl.dropboxusercontent.com/u/66785227/octanol_GMX.itp
https://dl.dropboxusercontent.com/u/66785227/API_gmx.itp
https://dl.dropboxusercontent.com/u/66785227/topol.top
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Free energy, domain decomposition and gmx 465 v 505

[Gmx QA]

Hi all

I am running some free energy tutorial calculations to learn more about how
those are set up in gromacs, but have run into a problem when using
different versions.

I am following the free energy tutorial on the Gromacs website
(fe-tutorial-4.6), so not Justin's tutorial.

With vmx 4.6.5 everything goes according to the tutorial, and I can run
with four thread on my laptop.
Please find what I think is the relevant section of the log file below:

$/usr/local/gmx465/bin/mdrun -v -deffnm topol_465

Initializing Domain Decomposition on 4 nodes

Dynamic load balancing: auto

Will sort the charge groups at every domain (re)decomposition

Initial maximum inter charge-group distances:

two-body bonded interactions: 0.239 nm, LJC-14 q, atoms 4 7

  multi-body bonded interactions: 0.239 nm, Angle, atoms 2 6

Minimum cell size due to bonded interactions: 0.263 nm

Maximum distance for 5 constraints, at 120 deg. angles, all-trans: 0.497 nm

Estimated maximum distance required for P-LINCS: 0.497 nm

This distance will limit the DD cell size, you can override this with -rcon

Using 0 separate PME nodes, as there are too few total

 nodes for efficient splitting

Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25

Optimizing the DD grid for 4 cells with a minimum initial size of 0.622 nm

The maximum allowed number of cells is: X 3 Y 3 Z 2

Domain decomposition grid 2 x 2 x 1, separate PME nodes 0

PME domain decomposition: 2 x 2 x 1

Domain decomposition nodeid 0, coordinates 0 0 0


Using 4 MPI threads



While on the other hand, when I am using v 5.0.5 I cannot start simulations
using more than one thread, otherwise the domain decomposition complains,
even though the excerpt from the log file from that run is (almost)
identical


$/usr/local/gmx505/bin/mdrun -v -deffnm topol_505


Initializing Domain Decomposition on 4 ranks

Dynamic load balancing: auto

Will sort the charge groups at every domain (re)decomposition

Initial maximum inter charge-group distances:

two-body bonded interactions: 0.239 nm, LJC-14 q, atoms 4 7

  multi-body bonded interactions: 0.239 nm, Angle, atoms 2 6

Minimum cell size due to bonded interactions: 0.263 nm

Maximum distance for 5 constraints, at 120 deg. angles, all-trans: 0.497 nm

Estimated maximum distance required for P-LINCS: 0.497 nm

This distance will limit the DD cell size, you can override this with -rcon

Using 0 separate PME ranks, as there are too few total

 ranks for efficient splitting

Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25

Optimizing the DD grid for 4 cells with a minimum initial size of 0.622 nm

The maximum allowed number of cells is: X 3 Y 3 Z 2


---

Program gmx, VERSION 5.0.5

Source code file:
/Users/pk/source/gromacs-5.0.5/src/gromacs/mdlib/domdec.c, line: 6902


Fatal error:

There is no domain decomposition for 4 ranks that is compatible with the
given box and a minimum cell size of 0.62175 nm

Change the number of ranks or mdrun option -rcon or -dds or your LINCS
settings

Look in the log file for details on the domain decomposition

For more information and tips for troubleshooting, please check the GROMACS

website at http://www.gromacs.org/Documentation/Errors

---



Is there something that I have overlooked in v 5 that needs to be taken
into consideration? I am not that familiar with the new parallelization
scheme of v5, maybe there is something there?


Thanks

/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Reference value for berger POPC area per headgroup

[Gmx QA]

Hi

I have a question about which value to aim for as a reference for the
berger lipids and POPC area per head group.

In Justin Lemkul's tutorial, the value 65.8 A^2 is reported, along with a
reference to a 1998 Tieleman et al. paper.
However, I have read that paper now twice and cannot find this value
anywhere.
There are also two popc-systems available for download from the Tieleman
website, but neither of those give 65.8 A^2, at least in my hands.

Could someone please enlighten me about a good reference value for berger
popc apl?

Thanks
/PK
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.