Re: [gmx-users] RE: Gibbs Energy Calculation and charges

[Michael Shirts]

I likely won't have much time to look at it tonight, but you can see
exactly what the option is doing to the topology.  run gmxdump on the
tpr.  All of the stuff that couple-intramol does is in grompp, so the
results will show up in the detailed listings of the interactions, and
which ones have which values set for the A and B states.

On Wed, Oct 30, 2013 at 5:36 PM, Dallas Warren  wrote:
> Michael, thanks for taking the time to comment and have a look.
>
> The real issue I am having is a bit deeper into the topic than that, my last 
> reply was just an observation on something else.  Will summarise what I have 
> been doing etc.
>
> I have a molecule that are calculating the Gibbs energy of hydration and 
> solvation (octanol).  In a second topology the only difference is that the 
> atomic charges have been doubled.  Considering that charges are scaled 
> linearly with lambda, the normal charge values of dH/dl from lambda 0 to 1 
> obtained should reproduce that of the double charged molecule from lambda 0.5 
> to 1.0.  Is that a correct interpretation?  Since the only difference should 
> be that charge of the atoms and over that range the charge will be identical.
>
> I was using couple-intramol = no and the following are the results from those 
> simulations.
>
> For the OE atom within the molecule, I have plotted the following graphs of 
> dH/dl versus charge of that atom for both of the topologies.
> octanol - http://ozreef.org/stuff/octanol.gif
> water - http://ozreef.org/stuff/water.gif
> mdp file - http://ozreef.org/stuff/gromacs/mdout.mdp
>
> The mismatch between the two topologies is the real issue that I am having.  
> I was hoping to get the two to overlap.
>
> My conclusion based on this is that there is actually something else being 
> changed with the topology by GROMACS when the simulations are being run.  The 
> comments in the manual allude to that, but not entirely sure what is going on.
>
>> From the manual:
>>
>>couple-intramol:
>>
>>no
>> All intra-molecular non-bonded interactions for moleculetype
>>couple-moltype are replaced by exclusions and explicit pair
>>interactions. In this manner the decoupled state of the molecule
>>corresponds to the proper vacuum state without periodicity effects.
>>yes
>> The intra-molecular Van der Waals and Coulomb interactions are
>>also turned on/off. This can be useful for partitioning free-energies
>>of relatively large molecules, where the intra-molecular non-bonded
>>interactions might lead to kinetically trapped vacuum conformations.
>>The 1-4 pair interactions are not turned off.
>
> Chris Neale commented:
>>Ah, I see. I guess that you are using couple-intramol = no (the default in 
>>v4.6.3 at least). That means
>>that the intramolecular charge-charge interactions are always at 
>>full-strength (and therefore different).
>>I would expect that normal at lambda=0 should be the same as double at 
>>lambda=0.5 only for
>>couple-intramol = yes
>>If you were using couple-intramol = yes already, then I am as confused as you 
>>are.
>
> >From Chris' comment, I then went back and repeated the calculation using 
> >couple-intramol = yes with the results of this in the below image (plus the 
> >previous when set to no for comparison)
>
> http://ozreef.org/stuff/gromacs/couple-intramol.png
>
> My rather garbled comment that you saw was my attempt to make sense of the 
> fact that this setting made things worse (and try to add something to get 
> this issue to have another pass on the list, hoping someone like yourself 
> will comment), and the value of dH/dl for c-i=yes at lambda=1 matches with 
> c-i=no at lambad=0.  This you can see with the previously linked to graphs.
>
> Hopefully that helps understand what I have been doing and the issues.
>
> Catch ya,
>
> Dr. Dallas Warren
> Drug Delivery, Disposition and Dynamics
> Monash Institute of Pharmaceutical Sciences, Monash University
> 381 Royal Parade, Parkville VIC 3052
> dallas.war...@monash.edu
> +61 3 9903 9304
> -
> When the only tool you own is a hammer, every problem begins to resemble a 
> nail.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] RE: Gibbs Energy Calculation and charges

[Michael Shirts]

I think the grammar got a little garbled there, so I'm not sure quite
what you are claiming.

One important thing to remember; 1-4 interactions are treated as
bonded interactions right now FOR COUPLE intramol (not for lambda
dependence of the potential energy function), so whether
couple-intramol is set to yes or no does not affect these interactions
at all.  It only affects the nonbondeds with distances greater than
1-5.  At least to me, this is nonintuitive (and we're coming up with a
better scheme for 5.0), but might that explain what you are getting?

On Tue, Oct 29, 2013 at 9:44 PM, Dallas Warren  wrote:
> Just want this to make another pass, just in case those in the know missed it.
>
> Using couple-intrmol = yes the resulting dH/dl plot actually looks like that 
> at lamba = 1 it is actually equal to couple-intramol = no with lambda = 0.
>
> Should that be the case?
>
> Catch ya,
>
> Dr. Dallas Warren
> Drug Delivery, Disposition and Dynamics
> Monash Institute of Pharmaceutical Sciences, Monash University
> 381 Royal Parade, Parkville VIC 3052
> dallas.war...@monash.edu
> +61 3 9903 9304
> -
> When the only tool you own is a hammer, every problem begins to resemble a 
> nail.
>
>> -Original Message-
>> From: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> boun...@gromacs.org] On Behalf Of Dallas Warren
>> Sent: Tuesday, 22 October 2013 2:49 PM
>> To: Discussion list for GROMACS users
>> Subject: [gmx-users] RE: Gibbs Energy Calculation and charges
>>
>> I have completed the simulations using couple-intramol = yes
>>
>> Reminder, have a molecule that are calculating the Gibbs energy of
>> hydration and solvation (octanol).  In one topology the only difference
>> is the atomic charges have been doubled.  Considering that charges are
>> scaled linearly with lambda, the normal charge values of dH/dl from
>> lambda 0 to 1 obtained should reproduce that of the double charged
>> molecule from lambda 0.5 to 1.0.
>>
>> Here is the mdout.mdp file from one of the simulations using couple-
>> intramol = yes.  The difference with the simulations run using couple-
>> intramol = no is only that setting, checked using the diff command.
>> Copy of the mdout.mdp file can be found:
>>
>> http://ozreef.org/stuff/gromacs/mdout.mdp
>>
>> Changing to couple-intrmol = yes makes things significantly worse.
>>
>> Here are the plots of dH/dl obtained as a function of the charge of one
>> of the atoms, OE (calculated based on fact that at lambda = 0 fully
>> charged, lambda = 1 no charged)
>>
>> http://ozreef.org/stuff/gromacs/couple-intramol.png
>>
>> Anyone with some insight on what is going on here?
>>
>> >From the manual:
>> 
>> couple-intramol:
>>
>> no
>> All intra-molecular non-bonded interactions for moleculetype
>> couple-moltype are replaced by exclusions and explicit pair
>> interactions. In this manner the decoupled state of the molecule
>> corresponds to the proper vacuum state without periodicity effects.
>> yes
>> The intra-molecular Van der Waals and Coulomb interactions are
>> also turned on/off. This can be useful for partitioning free-energies
>> of relatively large molecules, where the intra-molecular non-bonded
>> interactions might lead to kinetically trapped vacuum conformations.
>> The 1-4 pair interactions are not turned off.
>> 
>>
>> Catch ya,
>>
>> Dr. Dallas Warren
>> Drug Delivery, Disposition and Dynamics
>> Monash Institute of Pharmaceutical Sciences, Monash University
>> 381 Royal Parade, Parkville VIC 3052
>> dallas.war...@monash.edu
>> +61 3 9903 9304
>> -
>> When the only tool you own is a hammer, every problem begins to
>> resemble a nail.
>>
>>
>> > -Original Message-
>> > From: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> > boun...@gromacs.org] On Behalf Of Dallas Warren
>> > Sent: Thursday, 17 October 2013 3:32 PM
>> > To: Discussion list for GROMACS users
>> > Subject: [gmx-users] RE: Gibbs Energy Calculation and charges
>> >
>> > Chris,
>> >
>> > Thank you, that appears to be the issue then.
>> >
>> > Running them again now with couple-intramol = yes
>> >
>> > Will report back once that is completed with the results.
>> >
>> > Catch ya,
>> >
>> > Dr. Dallas Warren
>> > Drug Delivery, Disposition and Dynamics
>> > Monash Institute of Pharmaceutical Sciences, Monash University
>> > 381 Royal Parade, Parkville VIC 3052
>> > dallas.war...@monash.edu
>> > +61 3 9903 9304
>> > -
>> > When the only tool you own is a hammer, every problem begins to
>> > resemble a nail.
>> >
>> >
>> > > -Original Message-
>> > > From: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> > > boun...@gromacs.org] On Behalf Of Christopher Neale
>> > > Sent: Thursday, 17 October 2013 2:15 PM
>> > > To: gmx-users@gromacs.org
>> > > Subject: [gmx-users] Gibbs Energy Calculation and charges
>> > >
>> > > Ah, I see. I 

Re: [gmx-users] lmc-stats

[Michael Shirts]

Hi, Andrew-

The choice of lmc-stats only affects the calculation of the weights,
it does not affect the calculation of the transition matrix.

There are two possibilities for the transition matrix; the estimated
one (which is just called 'Transition Matrix' -- we should probably
have a better name), and the empirical one.  The empirical transition
matrix is just transition counts. The first one is an ensemble average
of P(k|x); for neighbor moves, this is just the average of the barker
transition probabilities (plus the choice of going up or down),

Let me know if this is useful.  Not that many people have used this
code, thus there are likely ways that it can be improved for better
utility.  I generally haven't tried to calculate free energy
differences from the transition matrix, though it should give
consistent results with the other methods.

Best,
~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shi...@virginia.edu
(434)-243-1821

On Mon, Oct 28, 2013 at 11:13 PM, Andrew S. Paluch  wrote:
> When performing free energy calculations using the expanded ensemble method,
> there is a barker and metropolis option for lmc-stats. Are the corresponding
> transition probabilities computed with or without the weighting factors?
> That is, are these probabilities biased or not?
>
> Thank you,
>
> Andrew
>
> --
>
> Andrew S. Paluch, PhD
> Department of Chemical, Paper, and Biomedical Engineering
> Miami University
> paluc...@miamioh.edu
> (513) 529-0784
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Replica Exchange with Solute Tempering

[Michael Shirts]

Hi, all-

Rest essentially scales the solute-solvent interactions, but maintains
the solute-solute interactions. This can be done solely with
Hamiltonian replica exchange, which is in 4.6.  It's a bit tricky,
though.  We plan on having something that does this automatically in
5.0 or 5.1, but it's not there yet.

What do you intend to do?  It could be that there's a better way to do
what you want to do with Hamiltonian replica exchange, as well.



On Sat, Oct 26, 2013 at 3:52 PM, David Osguthorpe
 wrote:
>
> Hi
>
> Ive been looking for a sample input for a REST simulation but
> so far have not been able to find one.
>
> Does anybody have an example to share?
>
> Ive seen the tutorial by Mark which mentions REST but I dont
> see this included in the archived examples.
>
> Am I missing something?
>
> Thanks for any help
>
> David
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] energy drift - comparison of double and single precision

[Michael Shirts]

Hi, all-

At this point, any fixes are going to be in the 5.0 version, where the
integrators will be a bit different.  If you upload your system files
to redmine.gromacs.org (not just the .mdp), then I will make sure this
gets tested there.

On Fri, Oct 25, 2013 at 10:14 AM, Guillaume Chevrot
 wrote:
> Dear Xavier,
>
> 2013/10/12 XAvier Periole 
>
>>
>> Could you try to reduce the nstcalcenergy flag from 100 to 10 and then one?
>>
>>
> FYI, I tried 10 and 1 and the energy drift is exactly the same.
>
>
>
>
>>  Similar flags apply to temperature and pressure and I believe might
>> seriously affect energy conservation.
>>
>> XAvier.
>>
>> > On Oct 12, 2013, at 0:50, Mark Abraham  wrote:
>> >
>> > Didn't see any problem in the .mdp. -4500 kJ/mol in 10ns over (guessing)
>> > 30K atoms is 0.015 kJ/mol/ns/atom. k_B T is at least 100 times larger
>> than
>> > that. I bet the rest of the lysozyme model physics is not accurate to
>> less
>> > than 1% ;-) There are some comparative numbers at
>> > http://dx.doi.org/10.1016/j.cpc.2013.06.003 - the two systems are rather
>> > different but they share the use of SETTLE.
>> >
>> > Note that using md-vv guarantees the 2007 paper is inapplicable, because
>> > GROMACS did not have a velocity Verlet integrator back then. Sharing the
>> > .log files might be informative.
>> >
>> > Mark
>> >
>> >
>> > On Fri, Oct 11, 2013 at 11:38 PM, Guillaume Chevrot <
>> > guillaume.chev...@gmail.com> wrote:
>> >
>> >> Hi,
>> >>
>> >> sorry for my last post! I re-write my e-mail (with some additional
>> >> information) and I provide the links to my files ;-)
>> >>
>> >> I compared the total energy of 2 simulations:
>> >> lysozyme in water / NVE ensemble / single precision / Gromacs 4.6.3
>> >> lysozyme in water / NVE ensemble / double precision / Gromacs 4.6.3
>> >>
>> >> ... and what I found was quite ... disturbing (see the plots of the
>> total
>> >> energy: http://dx.doi.org/10.6084/m9.figshare.820153). I observe a
>> >> constant
>> >> drift in energy in the case of the single precision simulation.
>> >>
>> >> Did I do something wrong*? Any remarks are welcomed! Here is the link to
>> >> the ‘mdout.mdp’ file (http://dx.doi.org/10.6084/m9.figshare.820154) so
>> you
>> >> can check what mdp options I used.
>> >>
>> >> My second question is: if I did not do something wrong, what are the
>> >> consequences on the simulation? Can I trust the results of single
>> precision
>> >> simulations?
>> >>
>> >> Regards,
>> >>
>> >> Guillaume
>> >>
>> >> *PS: I am not the only one encountering this behavior. In the
>> literature,
>> >> this problem has already been mentioned:
>> >> http://jcp.aip.org/resource/1/jcpsa6/v126/i4/p046101_s1
>> >>
>> >>
>> >>
>> >>
>> >> 2013/10/11 Mark Abraham 
>> >>
>> >>> On Oct 11, 2013 7:59 PM, "Guillaume Chevrot" <
>> >> guillaume.chev...@gmail.com>
>> >>> wrote:
>> 
>>  Hi all,
>> 
>>  I recently compared the total energy of 2 simulations:
>>  lysozyme in water / NVE ensemble / single precision
>>  lysozyme in water / NVE ensemble / double precision
>> 
>>  ... and what I found was quite ... disturbing (see the attached figure
>> >> -
>>  plots of the total energy). I observe a constant drift in energy in
>> the
>>  case of the single precision simulation.
>> 
>>  Did I do something wrong*? Any remarks are welcomed! I join the
>> >>> ‘mdout.mdp’
>>  file so you can check what mdp options I used.
>> >>>
>> >>> Maybe. Unfortunately we cannot configure the mailing list to allow
>> people
>> >>> to send attachments to thousands of people, so you will need to do
>> >>> something like provide links to files on a sharing service.
>> >>>
>> 
>>  My second question is: if I did not do something wrong, what are the
>>  consequences on the simulation? Can I trust the results of single
>> >>> precision
>>  simulations?
>> >>>
>> >>> Yes, as you have no doubt read in the papers published by the GROMACS
>> >> team.
>> >>>
>>  Regards,
>> 
>>  Guillaume
>> 
>>  *PS: I am not the only one encountering this behavior. In the
>> >> literature,
>>  this problem has already been mentioned:
>>  http://jcp.aip.org/resource/1/jcpsa6/v126/i4/p046101_s1
>> >>>
>> >>> ... which is six years old, examining the properties of code seven
>> years
>> >>> old. Life has moved on! :-) Even if you have found a problem, it is a
>> big
>> >>> assumption that this is (still) the cause.
>> >>>
>> >>> Mark
>> >>>
>>  --
>>  gmx-users mailing listgmx-users@gromacs.org
>>  http://lists.gromacs.org/mailman/listinfo/gmx-users
>>  * Please search the archive at
>> >>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>  * Please don't post (un)subscribe requests to the list. Use the
>>  www interface or send it to gmx-users-requ...@gromacs.org.
>>  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >>> --
>> >>> gmx-users mailing listgmx-users@gr

Re: [gmx-users] a new GROMACS simulation tool

[Michael Shirts]

Is there a link to the documentation?  It's a little difficult to know
exactly what this supposed to be doing.  Is it a GUI interface to
gromacs?

In general, it would be great to get these sort of extensions
coordinated with the main gromacs development tree, since otherwise
they would tend to get out of sync, lowering utility to everyone.

On Tue, Oct 22, 2013 at 10:34 AM, Kevin Chen  wrote:
> Hi Everyone,
>
> I'm writing to let you guys know that we have developed a web-based tool MD
> simulation tool for GROMACS.  It is a software package primarily developed
> for biological MD and offers a huge amount of possible options and settings
> for tailoring the simulations. Seamlessly integrated with newly developed
> GUI interfaces, the tool provides comprehensive setup, simulation, analysis
> and job submission tools. Most importantly, unlike other GROMACS GUI
> applications, user can actually run really simulations using the dedicated
> HPC resources. That been said, there's no proposal and installation
> required.  This tool could be a great fit for both teaching and research
> projects. Users inexperienced in MD can work along prepared workflows, while
> experts may enjoy a significant relief from the tedium of typing and
> scripting. As for now, we'd like to invite people to participate in user
> testing on this newly developed tool. Let me know if you'd like to try it
> out. We will set up an account for you.
>
> Best Regards,
>
> Kevin Chen, Ph.D.
> Information Technology at Purdue (ITaP)
> West Lafayette, IN 47907-2108
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] energy drift - comparison of double and single precision

[Michael Shirts]

Hi, Guillaume-

No one can tell you if you did anything wrong if you didn't tell us
what you did!  There are literally thousands of combinations of
options in running an NVE simulation, a substantial fraction of which
are guaranteed not to conserve energy.

If you post the files (inputs and relevant output files) to
redmine.gromacs.org, then people would be able to see if it was an
error in how you ran the simulation, or an actual bug.



On Fri, Oct 11, 2013 at 1:58 PM, Guillaume Chevrot
 wrote:
> Hi all,
>
> I recently compared the total energy of 2 simulations:
> lysozyme in water / NVE ensemble / single precision
> lysozyme in water / NVE ensemble / double precision
>
> ... and what I found was quite ... disturbing (see the attached figure -
> plots of the total energy). I observe a constant drift in energy in the
> case of the single precision simulation.
>
> Did I do something wrong*? Any remarks are welcomed! I join the ‘mdout.mdp’
> file so you can check what mdp options I used.
>
> My second question is: if I did not do something wrong, what are the
> consequences on the simulation? Can I trust the results of single precision
> simulations?
>
> Regards,
>
> Guillaume
>
> *PS: I am not the only one encountering this behavior. In the literature,
> this problem has already been mentioned:
> http://jcp.aip.org/resource/1/jcpsa6/v126/i4/p046101_s1
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] question about OPLS-AA force field -required bond constraints

[Michael Shirts]

OPLS-AA was generally derived with Monte Carlo, which means that all
bonds were exactly constrained.  But read the papers!

On Thu, Oct 10, 2013 at 12:27 PM, Justin Lemkul  wrote:
>
>
> On 10/10/13 12:17 PM, Martin, Erik W wrote:
>>
>>
>> Hi, I'm new to both Gromacs and OPLS.  I have always used CHARMM or AMBER
>> before so am a bit confused about some idiosyncrasies.
>>
>> When I go to a production run on my system (a protein in a truncated
>> octahedron box of Tip4 water) I have removed all constraints.  When I run
>> grompp, I get an error telling me that a bond between a HG proton and OG
>> oxygen of protein has an oscillation period of 9e-3 ps and is incompatible
>> with my 1 fs time step.  I am used to being able to run a completely
>> unconstrained simulation with a 1 fs time step… so I'm assuming this is a
>> property of OPLS-AA?
>>
>
> No, it's not.  It's a property of the bonds.  Nearly all force fields (if
> not all of them) constrain X-H bonds at minimum, to remove the fastest
> motions in the system and allow for a longer time step.  Using dt on the
> same order as the fastest frequencies leads to instability.  See manual
> section 6.7.
>
>
>> Could someone please let me know what the bond constraint requirements are
>> for OPLS-AA and GROMOS forcefields (I'm going to need to use that next).
>> I
>> can't seem to find this information by searching but imagine its on a
>> website
>> somewhere?
>>
>
> Information like this should be in the primary literature for the parameter
> sets.  I don't know of any sort of database online that tracks such things,
> and moreover there are no absolute rules AFAIK.  I generally follow the
> methods of the paper that derived the parameter set.  Unfortunately,
> sometimes they don't say.  For Gromos96 53A5/53A6, validation simulations
> were carried out with all bonds constrained.
>
> -Justin
>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
>
> ==
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Problem running free energy simulations

[Michael Shirts]

Hmm. This really isn't quite enough information to go on.  Can you
file a redmine issue, and include the files used to generate the run
that crashed (.mdp, .gro,. .top), as well as the files that show it
failing (.log)?

http://redmine.gromacs.org/

On Thu, Oct 3, 2013 at 4:32 AM, Jernej Zidar  wrote:
> Dear Michael.
>   The simulations at each lambda point starts from the same structure
> that I equilibrated (NPT ensemble) for 20 nanoseconds. The system has
> ~7500 atoms in a box the size 5 nm x 5 nm x 5 nm.  The molecule of
> interest is located in the center of the unit cell.
>
> Thanks,
> JErnej
>>
>> Sounds like the simulation is blowing up.  How soon does it start crashing.
>>
>> Also, what configurations are you using to start your free energy
>> simulations at each lambda?
>
>
>
>
> --
> Windows: Re-Boot, Linux: Be-Root.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Problem running free energy simulations

[Michael Shirts]

Sounds like the simulation is blowing up.  How soon does it start crashing.

Also, what configurations are you using to start your free energy
simulations at each lambda?

On Wed, Oct 2, 2013 at 10:31 PM, Jernej Zidar  wrote:
> Hi all,
>   I'm trying to determine the free energy of solvation for a molecule
> in n-octanol. I'm separately turning off the Coulomb and Lennard-Jones
> interactions as instructed in the free energy tutorial. The
> Lennard-Jones simulations keep on crashing for most values of lambda
> with the message:
> Program mdrun, VERSION 4.6.3
> Source code file: /home/zidar/gromacs-4.6.3/src/mdlib/domdec_top.c, line: 393
>
> Fatal error:
> 1 of the 79007 bonded interactions could not be calculated because
> some atoms involved moved further apart than the multi-body cut-off
> distance (1 nm) or the two-body cut-off distance (1 nm), see option
> -rdd, for pairs and tabulated bonds also see option -ddcheck
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
>
> - - -
>
>   I tried using "-noddcheck", "-pd" (separate simulations) but the
> simulations still fail.
>
>   The molecule surveyed is not charged (branched hydrocarbon) and I
> can simulate it using any ensemble without a problem. I'm using CgenFF
> without adding any new atom types so all the interactions should be
> properly accounted for.
>
>   Any advice on how to solve this will be greatly appreciated. I was
> able to perform free energy simulations of a similar molecule in
> octanol without any major problem.
>
> Thanks,
> Jernej
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] "Illegal instruction" error from alchemical-gromacs.py

[Michael Shirts]

Hi, Chris-

The best place to file this issue is the SimTK pymbar page, rather
than alchemistry.org, since it's a pymbar problem.

We have collaborators that may have updated the pymbar.py recently.
I'll try to get this stabilized in the very near future.

Testing quickly, my best guess is that it's not alchemical-gromacs.py
(latest version seems to work fine for me), but the pymbar
installation that is causing problems.  Give the illegal instruction
calculations, the helper C++ code is probably causing a problem;
perhaps it compiled using headers from the wrong version of numpy.

In most cases, the best way is to explicitly force it to use the pure
python. If you just delete the compiled _pymbar.so, then it will
default to pure Python.   It may be slower (5-10x) but it should still
finish in a minute or two.

I'm working on getting instructions for g_bar up there in
alchemistry.org as well, but there are some bugs in g_bar that I am
dealing with.

On Thu, Sep 26, 2013 at 10:35 PM, Christopher Neale
 wrote:
> Dear Users:
>
> I'm having difficulty running MBAR after some free energy calculations (MBAR 
> via alchemical-gromacs.py obtained from alchemistry.org).
> The input options to alchemical-gromacs.py have obviously changed since the 
> site at
> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy
> was updated, but I've taken a look at the python source and it seems as if my 
> invocation below should work.
> Nevertheless, I get an "Illegal instruction" error (See below).
>
> gpc-logindm02-$ ls ANALYSIS/
> ethanol.0.dhdl.xvg  ethanol.2.dhdl.xvg  ethanol.4.dhdl.xvg  
> ethanol.6.dhdl.xvg  ethanol.8.dhdl.xvg
> ethanol.1.dhdl.xvg  ethanol.3.dhdl.xvg  ethanol.5.dhdl.xvg  ethanol.7.dhdl.xvg
>
> gpc-logindm02-$ python 
> /project/p/pomes/cneale/GPC/exe/pymbar/trunk/examples/alchemical-free-energy/alchemical-gromacs.py
>  -d ANALYSIS -p ethanol. -t 300 -s 1000 -v
> Warning on use of the timeseries module: If the inherent timescales of the 
> system are long compared to those being analyzed, this statistical 
> inefficiency may be an underestimate.  The estimate presumes the use of many 
> statistically independent samples.  Tests should be performed to assess 
> whether this condition is satisfied.   Be cautious in the interpretation of 
> the data.
> Started on Thu Sep 26 22:32:05 2013
> Illegal instruction
>
> (omitting the -t -s and -v flags yields the same error).
>
> I realize that I'm hijacking the gromacs list a little bit here but there is 
> no mailing list for http://www.alchemistry.org/
> and it seems like a good idea to have these Q&A's archived.
>
> Thank you,
> Chris.
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] segfault when running the alchemistry.org ethanol solvation free energy tutorial

[Michael Shirts]

Sounds like we've resolved all the confusion.  Thanks for prompt help
in making this clearer and better.

On Thu, Sep 26, 2013 at 7:01 PM, Christopher Neale
 wrote:
> Agreed, the following parameters do not segfault in single or double 
> precision:
> sc-alpha = 0.5
> sc-power = 1
> sc-r-power   = 6
> Same goes for http://bugzilla.gromacs.org/issues/1306
>
> The following parameters give a segfault in single precision but are ok in 
> double precision
> sc-alpha = 0.001
> sc-power = 1
> sc-r-power   = 48
> Same goes for http://bugzilla.gromacs.org/issues/1306
>
> Sorry for the confusion earlier. I was using a compilation that I thought was 
> double precision
> but it was actually single. Recompiling in double precision gave me the 
> stability outlined above.
>
> Thank you for your assitance Mark and Michael.
>
> Chris.
>
> -- original message --
>
> Michael Shirts mrshirts at gmail.com
> Fri Sep 27 00:41:17 CEST 2013
> Previous message: [gmx-users] segfault when running the alchemistry.org 
> ethanol solvation free energy tutorial
> Next message: [gmx-users] segfault when running the alchemistry.org ethanol 
> solvation free energy tutorial
> Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
> And with this change, single is running fine as well.
>
> This was a known issue, but was only documented in the expanded.mdp
> files, which was an oversight.  After this, I switched so the default
> is less likely to cause problems.  Because of some theory improvements
> developed in the group in free energy calculation pathways, the
> sc-r-power=48 pathway will now be phased out anyway by 5.1.
>
> On Thu, Sep 26, 2013 at 6:37 PM, Michael Shirts  wrote:
>> I thought I had just managed to solve the issue :)
>>
>> If you look at the soft core parameters, there are two types listed --
>> one with sc-r-power = 48, and one with sc-r-power = 6.  The sc-r-power
>> are more stable with single precision calculations.
>>
>> I have changed the files on the website to make the single precision
>> ones the default.  Expanded.mdp warned about the issues with
>> precision, but left the sc-r-power in place; the Ethanol.mdp did not
>> warn about the potential issue.  Now both include the more efficient
>> path commented out.
>>
>> NOTE that this means the exact numbers of the tutorials are not quite
>> right anymore; the process is unchanged, as is the final answer, but
>> the intermediate dG's will be slightly different.
>>
>> However, I need to redo them anyway to make them easier in the next
>> 1-2 weeks, so I will update them then.
>>
>> If it still fails with double, that's a different issue -- because I'm
>> running them fine with double.
>>
>>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] segfault when running the alchemistry.org ethanol solvation free energy tutorial

[Michael Shirts]

And with this change, single is running fine as well.

This was a known issue, but was only documented in the expanded.mdp
files, which was an oversight.  After this, I switched so the default
is less likely to cause problems.  Because of some theory improvements
developed in the group in free energy calculation pathways, the
sc-r-power=48 pathway will now be phased out anyway by 5.1.

On Thu, Sep 26, 2013 at 6:37 PM, Michael Shirts  wrote:
> I thought I had just managed to solve the issue :)
>
> If you look at the soft core parameters, there are two types listed --
> one with sc-r-power = 48, and one with sc-r-power = 6.  The sc-r-power
> are more stable with single precision calculations.
>
> I have changed the files on the website to make the single precision
> ones the default.  Expanded.mdp warned about the issues with
> precision, but left the sc-r-power in place; the Ethanol.mdp did not
> warn about the potential issue.  Now both include the more efficient
> path commented out.
>
> NOTE that this means the exact numbers of the tutorials are not quite
> right anymore; the process is unchanged, as is the final answer, but
> the intermediate dG's will be slightly different.
>
> However, I need to redo them anyway to make them easier in the next
> 1-2 weeks, so I will update them then.
>
> If it still fails with double, that's a different issue -- because I'm
> running them fine with double.
>
>
>
>
>
> On Thu, Sep 26, 2013 at 6:32 PM, Christopher Neale
>  wrote:
>> My mistake  I still get a segfault even when using double precision. (EM 
>> doesn't help, nor does switching to Berendsen pressure coupling).
>>
>> Note that I can stop the segfault when running at init-lambda-state = 0 if I 
>> set:
>>
>> couple-lambda0   = none
>> couple-lambda1   = none
>>
>> instead of
>>
>> couple-lambda0   = vdw-q
>> couple-lambda1   = none
>>
>> This looks like http://bugzilla.gromacs.org/issues/1306
>>
>> Thank you,
>> Chris.
>>
>> -- original message --
>>
>> Thank you Mark.
>>
>> I actually found that it crashed wihout the -multi part (no hamiltonian 
>> exchange).
>> The command that I used was: mdrun -nt 1 -deffnm ethanol.1 -dhdl 
>> ethanol.1.dhdl.xvg
>>
>> If I use the double precision version, there is no segfault. That's a 
>> working solution, but it is worrysome.
>>
>> Thank you,
>> Chris.
>>
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] segfault when running the alchemistry.org ethanol solvation free energy tutorial

[Michael Shirts]

I thought I had just managed to solve the issue :)

If you look at the soft core parameters, there are two types listed --
one with sc-r-power = 48, and one with sc-r-power = 6.  The sc-r-power
are more stable with single precision calculations.

I have changed the files on the website to make the single precision
ones the default.  Expanded.mdp warned about the issues with
precision, but left the sc-r-power in place; the Ethanol.mdp did not
warn about the potential issue.  Now both include the more efficient
path commented out.

NOTE that this means the exact numbers of the tutorials are not quite
right anymore; the process is unchanged, as is the final answer, but
the intermediate dG's will be slightly different.

However, I need to redo them anyway to make them easier in the next
1-2 weeks, so I will update them then.

If it still fails with double, that's a different issue -- because I'm
running them fine with double.





On Thu, Sep 26, 2013 at 6:32 PM, Christopher Neale
 wrote:
> My mistake  I still get a segfault even when using double precision. (EM 
> doesn't help, nor does switching to Berendsen pressure coupling).
>
> Note that I can stop the segfault when running at init-lambda-state = 0 if I 
> set:
>
> couple-lambda0   = none
> couple-lambda1   = none
>
> instead of
>
> couple-lambda0   = vdw-q
> couple-lambda1   = none
>
> This looks like http://bugzilla.gromacs.org/issues/1306
>
> Thank you,
> Chris.
>
> -- original message --
>
> Thank you Mark.
>
> I actually found that it crashed wihout the -multi part (no hamiltonian 
> exchange).
> The command that I used was: mdrun -nt 1 -deffnm ethanol.1 -dhdl 
> ethanol.1.dhdl.xvg
>
> If I use the double precision version, there is no segfault. That's a working 
> solution, but it is worrysome.
>
> Thank you,
> Chris.
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] segfault when running the alchemistry.org ethanol solvation free energy tutorial

[Michael Shirts]

Just to be clear, is this the expanded ensemble version of the calculation?

On Thu, Sep 26, 2013 at 5:25 PM, Mark Abraham  wrote:
> I found the -multi version of that tutorial a bit temperamental...
> Michael Shirts suggested that double precision is more reliable for
> expanded ensemble. Hopefully he can chime in in a day or two.
>
> Mark
>
> On Thu, Sep 26, 2013 at 9:00 PM, Christopher Neale
>  wrote:
>> Dear Users:
>>
>> Has anyone successfully run the free energy tutorial at 
>> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy
>>  ?
>>
>> I just tried it and I get a segmentation fault immediately (see output at 
>> the end of this post).
>>
>> I get a segfault with both 4.6.3 and 4.6.1.
>>
>> Note that if I modify the .mdp file to set free-energy = no , then the 
>> simulation runs just fine. (I have, of course, set init-lambda-state in the 
>> .mdp file that I downloaded from the aforementioned site and I get a 
>> segfault with any value of init-lambda-state from 0 to 8).
>>
>> gpc-f103n084-$ mdrun -nt 1 -deffnm ethanol.1 -dhdl ethanol.1.dhdl.xvg
>>  :-)  G  R  O  M  A  C  S  (-:
>>
>>   GROup of MAchos and Cynical Suckers
>>
>> :-)  VERSION 4.6.3  (-:
>>
>> Contributions from Mark Abraham, Emile Apol, Rossen Apostolov,
>>Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,
>>  Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans,
>> Peter Kasson, Carsten Kutzner, Per Larsson, Pieter Meulenhoff,
>>Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz,
>> Michael Shirts, Alfons Sijbers, Peter Tieleman,
>>
>>Berk Hess, David van der Spoel, and Erik Lindahl.
>>
>>Copyright (c) 1991-2000, University of Groningen, The Netherlands.
>>  Copyright (c) 2001-2012,2013, The GROMACS development team at
>> Uppsala University & The Royal Institute of Technology, Sweden.
>> check out http://www.gromacs.org for more information.
>>
>>  This program is free software; you can redistribute it and/or
>>modify it under the terms of the GNU Lesser General Public License
>> as published by the Free Software Foundation; either version 2.1
>>  of the License, or (at your option) any later version.
>>
>> :-)  mdrun  (-:
>>
>> Option Filename  Type Description
>> 
>>   -s  ethanol.1.tpr  InputRun input file: tpr tpb tpa
>>   -o  ethanol.1.trr  Output   Full precision trajectory: trr trj cpt
>>   -x  ethanol.1.xtc  Output, Opt. Compressed trajectory (portable xdr format)
>> -cpi  ethanol.1.cpt  Input, Opt.  Checkpoint file
>> -cpo  ethanol.1.cpt  Output, Opt. Checkpoint file
>>   -c  ethanol.1.gro  Output   Structure file: gro g96 pdb etc.
>>   -e  ethanol.1.edr  Output   Energy file
>>   -g  ethanol.1.log  Output   Log file
>> -dhdl ethanol.1.dhdl.xvg  Output, Opt! xvgr/xmgr file
>> -field  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>> -table  ethanol.1.xvg  Input, Opt.  xvgr/xmgr file
>> -tabletf  ethanol.1.xvg  Input, Opt.  xvgr/xmgr file
>> -tablep  ethanol.1.xvg  Input, Opt.  xvgr/xmgr file
>> -tableb  ethanol.1.xvg  Input, Opt.  xvgr/xmgr file
>> -rerun  ethanol.1.xtc  Input, Opt.  Trajectory: xtc trr trj gro g96 pdb cpt
>> -tpi  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>> -tpid ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>>  -ei  ethanol.1.edi  Input, Opt.  ED sampling input
>>  -eo  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>>   -j  ethanol.1.gct  Input, Opt.  General coupling stuff
>>  -jo  ethanol.1.gct  Output, Opt. General coupling stuff
>> -ffout  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>> -devout  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>> -runav  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>>  -px  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>>  -pf  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>>  -ro  ethanol.1.xvg  Output, Opt. xvgr/xmgr file
>>  -ra  ethanol.1.log  Output, Opt. Log file
>>  -rs  ethanol.1.log  Output, Opt. Log file
>>  -rt  ethanol.1.log  Output, Opt. Log file
>> -mtx  ethanol.1.mtx  Output, Opt. Hessian matrix
>>  -dn  ethanol.1.ndx  Output, Opt. Index file
>> -multidir ethanol.1  Input, Opt., Mult. Run directory
>> -membed  ethanol.1.dat  Input, Opt.  Generic data file
>>  -mp

Re: [gmx-users] Re: Need protein-ligand free energy calculation tutorial

[Michael Shirts]

There are some starter files here:

http://www.gromacs.org/Documentation/Tutorials/GROMACS_USA_Workshop_and_Conference_2013/An_introduction_to_free_energy_calculations%3a_Michael_Shirts%2c_Session_2A

Which can be used in conjunction with the Alchemistry.org instructions.

But it needs to be updated a bit more . . .

On Thu, Sep 19, 2013 at 3:38 AM, hsp85  wrote:
> Actually,
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/
> there are two different manuals - one about free energy calculation and
> other about protein ligand complexes.
>
> This page is also can be helpfull
> http://www.alchemistry.org/wiki/Category:Free_Energy_How-to%27s
>
> Sergey Filkin
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/Need-protein-ligand-free-energy-calculation-tutorial-tp5011299p5011302.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] A question about velocity rescaling thermostat and velocity verlet integrator

[Michael Shirts]

The only integrators with stochastic force components are sd and bd.

vrescale has a small stochastic component, but that is for the target
kinetic energy, and is not a random force acting on each particle.

On Thu, Aug 29, 2013 at 3:15 PM, Ali Sinan Saglam
 wrote:
> Hi,
>
> I was planning to use the velocity verlet ingetrator (md-vv) with velocity
> rescaling thermostat (tcoupl = v-rescale). I know that velocity rescaling
> thermostat is a stochastic thermostat and uses a random seed and generally
> with the stochastic thermostats one would use the sd integrator that uses
> ld-seed = -1 option for a random seed.
>
> I was just wondering if one uses the velocity verlet integrator with a
> stochastic thermostat does the ld-seed option work as it does with the sd
>  (or bd) integrator?  (the manual currently doesn't say anything about this
> http://manual.gromacs.org/online/mdp_opt.html#run)
>
> Best,
> Ali Sinan Saglam
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Long range Lennard Jones

[Michael Shirts]

IPS in CHARMM involves additional calculations beyond a simple
homogeneous approximation -- roughly equivalent to PME for dispersion,
though its a bit messier.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2723858/



On Thu, Aug 29, 2013 at 7:56 AM, Justin Lemkul  wrote:
>
>
> On 8/29/13 1:18 AM, Gianluca Interlandi wrote:
>>
>> Justin,
>>
>> I respect your opinion on this. However, in the paper indicated below by
>> BR
>> Brooks they used a cutoff of 10 A on LJ when testing IPS in CHARMM:
>>
>> Title: Pressure-based long-range correction for Lennard-Jones interactions
>> in
>> molecular dynamics simulations: Application to alkanes and interfaces
>> Author(s): Lague, P; Pastor, RW; Brooks, BR
>> Source: JOURNAL OF PHYSICAL CHEMISTRY B Volume: 108 Issue: 1 Pages:
>> 363-368
>> Published: JAN 8 2004
>>
>
> I cannot say for sure whether the DispCorr implementation in Gromacs is the
> same or comparable to IPS in CHARMM.  You would have to test that, and also
> test a more complex system than simple liquids.
>
>
>> There is also a paper by Piana and Shaw where different cutoffs for
>> non-bonded
>> are tested with CHARMM22 on Anton:
>>
>> http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0039918
>>
>> They found some subtle differences, in particular for cutoffs shorter than
>> 9 A.
>> However, Anton uses abrupt truncation (no switching) and I believe that
>> the
>> differences they found at cutoffs > 9 A would be much smaller if they had
>> used a
>> finer mesh (as they show at the 8 A cutoff). I always use
>> fourierspacing=1.0
>>
>
> Note that Shaw's group is using their own modified force field, CHARMM22*,
> so that factors in here, as well.  I do recall that paper, though I haven't
> read it in a while, so details may be fuzzy.  Wasn't the principal point to
> test methods for long-range electrostatics and the influence of cutoffs in
> that context?  I seem to recall LJ taking a backseat there.  It is certainly
> true that short-range Coulomb cutoffs are more flexible when using Ewald
> methods.
>
>
>> I agree though that it strongly depends on the system and I have always
>> run
>> control simulations but never found significant differences in the case of
>> just
>> proteins.
>>
>
> That's good to know that you've tested different setups.  Are the water
> interaction energies comparable?  Have you tested model compounds or just
> full proteins?  Over how long of a time period?
>
>
>> Finally, I have not tested it in gromacs, but in NAMD there is a
>> performance
>> gain of 25% when using the shorter cutoff. This is a factor to consider.
>> When I
>> asked for Teragrid supercomputing allocations back in 2006 and 2007 and I
>> suggested 10/12/14 cutoff, the reviewers always complained and cut my
>> requested
>> time of 20% with the justification that I must use a shorter cutoff.
>>
>
> That's unfortunate.  In my opinion, it is shameful that desired performance
> exceeds proven integrity of the data.  If all we're after is performance, I
> can write up a dozen papers in the next few months with completely useless
> data, but they will be produced quickly! ;)
>
>
> -Justin
>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
>
> ==
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: NPT-REMD

[Michael Shirts]

Can you clarify - Do you mean that different replicas have different
average pressures?

WITHIN each replica, the +/- 2000 bar changing from step to step is
very common for using an atomic virial like gromacs does.

The AVERAGES of EACH replica should each be the average pressure they
are set as (+/- much less than 2000 -- maybe 10-20 bar).

You could also try md-vv + MTTK + Nose-Hoover.  I have had decent luck
with this combination and tested it for a while.

You should also try running without any exchange.  It could be that
the problem is not replica exchange, but something else.



On Sun, Aug 25, 2013 at 8:36 PM, rahul.seth.gromacs
 wrote:
> Thanks everyone for your suggestions. Pressures for some of my replicas
> fluctuate between -2000 bars to +2000 bars. I do not think that can be
> termed as about right isn't it? The averages I am talking are over 1ns.
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/NPT-REMD-tp5010721p5010728.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] NPT-REMD

[Michael Shirts]

Pressure should fluctuate significantly.  The estimator for the
pressure that is generally used is very noisy. The question is, do the
pressure averages over, say, 500 ps or 1 ns look about right?

On Sun, Aug 25, 2013 at 5:07 PM, Mark Abraham  wrote:
> On Sun, Aug 25, 2013 at 6:22 PM, rahul seth
>  wrote:
>> Hi All,
>>
>> I have been trying to perform NPT-REMD with a Protein Substrate and water.
>> I am trying to use md+nose-hoover thermostat and parrinello-rahman
>> barostat. However, I am not quite sure about the tau-t and the tau-p that I
>> should be using here. I paste a part of my mdp file below:
>> define   =  -DPOSRES_subs
>> constraints=  all-bonds
>> pbc   =  xyz
>> integrator   =  md
>> ;ld_seed= %8i
>> dt =  0.002; ps !
>> nsteps   =  25   ; 50 ns
>> ;
>> nstcomm =  1
>> nstcalcenergy   =  10
>> nstxout   =  0
>> nstxtcout=  5000 ;every 10 ps
>> nstvout   =  0
>> nstfout   =  0
>> nstlog=  1000
>> nstenergy   =  1000
>> nstlist=  10
>> ns_type =  grid
>> rlist   =  1.0
>> vdwtype =  cut-off
>> rvdw  =  1.0
>> coulombtype =  pme
>> rcoulomb=  1.0
>> fourierspacing  = 0.12
>> pme_order   = 4
>> ewald_rtol   = 1e-5
>> optimize_fft =  yes
>> tcoupl= nose-hoover
>> nsttcouple  = 5
>> tc-grps   =  subs Protein SOL
>> tau_t  = 0.2 0.2 0.2
>> ref_t   =  %8.3f %8.3f %8.3f
>> ; Energy monitoring
>> energygrps  =  subs Protein SOL
>> ; Isotropic pressure coupling is now on
>> Pcoupl   = parrinello-rahman
>> nstpcouple  = 1
>> refcoord-scaling   = all
>> Pcoupltype  = isotropic
>> tau_p  = 2.0
>> compressibility = 4.5e-5
>> ref_p= 1.01325
>>
>>
>> Although the temperatures in this case stays close to the desired values,
>> the pressures of the individual replicas fluctuate significantly and hence
>> I do not think I have weird replica exchange probabilities.
>
> Seemingly large
> http://www.gromacs.org/Documentation/Terminology/Pressure fluctuations
> are normal - see link. You should satisfy yourself that your .mdp file
> produces ensembles that are respectable before you get involved with
> REMD, and gathering statistics for that will take longer for pressure
> than temperature.
>
>> I have varied tau-t and tau-p from 0.2 to 1 and 1 to 5 respectively. It
>> doesn't really seem to improve anything. Can anyone suggest a set of
>> optimized values for these parameters? I understand these maybe context
>> dependent but to give you a rough idea I have about ~18,000 atoms with 170
>> for atoms 1156 for substrate and the rest for water.
>
> Your literature survey of similar simulations should give you an idea
> of usage in the field, but (e.g. per GROMACS manual thermostat
> section) I would think tau-t for N-H should not go below 1. Otherwise,
> first observe whether you have a real problem before making haphazard
> changes :-)
>
> Mark
>
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Early registration period ending for 2013 GROMACS USA Workshop and Conference

[Michael Shirts]

Hi,  Rajat-

We are working on ways to make the workshop materials available
online, either streaming or posted afterwards.  The tutorials will
have online components as well. We were unable to obtain travel
funding for the conference this year, unfortunately!

> I would be willing to pay the whole workshop fee (and more) if you could 
> kindly video the whole workshop and make it available to us. Thank you.

We are not planning on charging for the material, but we should make
it easier to donate if you wanted to after we get this material out!


On Tue, Aug 20, 2013 at 12:14 AM, rajat desikan  wrote:
> Hi Prof. Shirts,
> I am sure that I speak for a lot of GROMACS users when I say that I want to
> attend the workshop, but I cannot. I would be willing to pay the whole
> workshop fee (and more) if you could kindly video the whole workshop and
> make it available to us. Thank you.
>
>
> On Tue, Aug 20, 2013 at 10:50 AM, Michael Shirts <
> michael.shi...@virginia.edu> wrote:
>
>> Dear GROMACS users-
>>
>> I'd like to remind you all about the 2013 GROMACS USA Workshop and
>> Conference at the University of Virginia Sept 13th-15th.  There are
>> still registration slots available, but the early registration
>> deadline of Aug 22nd is coming up in just a few days; after the 22nd,
>> the registration price will rise from the very low $60 for the two day
>> conference to the moderately-low-but-why-pay-more $95.
>>
>> As mentioned in previous emails to the list, the workshop and
>> conference will consist of two days of tutorials, discussions of
>> future plans for GROMACS, face time with many of the main GROMACS
>> developers, plenary software and application sessions, and an optional
>> last day of development planning to which attendees are also invited
>> to help influence the future directions of GROMACS.
>>
>> Please visit http://faculty.virginia.edu/gromacsworkshop for
>> registration and for much more information about the workshop,
>> including travel logistics. The website was also recently has been
>> updated with more specifics about the program.
>>
>> For specific questions about registering or logistics after visiting
>> the website, please write to
>> gromacsworkshop-registrat...@virginia.edu.
>>
>> Sincerely,
>> The 2013 GROMACS USA Workshop and Conference Steering Committee
>> Michael Shirts (chair)
>> Angel Garcia
>> Berk Hess
>> Yu-Shan Lin
>> Erik Lindahl
>> Peter Kasson
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>
>
>
> --
> Rajat Desikan (Ph.D Scholar)
> Prof. K. Ganapathy Ayappa's Lab (no 13),
> Dept. of Chemical Engineering,
> Indian Institute of Science, Bangalore
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Early registration period ending for 2013 GROMACS USA Workshop and Conference

[Michael Shirts]

Dear GROMACS users-

I'd like to remind you all about the 2013 GROMACS USA Workshop and
Conference at the University of Virginia Sept 13th-15th.  There are
still registration slots available, but the early registration
deadline of Aug 22nd is coming up in just a few days; after the 22nd,
the registration price will rise from the very low $60 for the two day
conference to the moderately-low-but-why-pay-more $95.

As mentioned in previous emails to the list, the workshop and
conference will consist of two days of tutorials, discussions of
future plans for GROMACS, face time with many of the main GROMACS
developers, plenary software and application sessions, and an optional
last day of development planning to which attendees are also invited
to help influence the future directions of GROMACS.

Please visit http://faculty.virginia.edu/gromacsworkshop for
registration and for much more information about the workshop,
including travel logistics. The website was also recently has been
updated with more specifics about the program.

For specific questions about registering or logistics after visiting
the website, please write to
gromacsworkshop-registrat...@virginia.edu.

Sincerely,
The 2013 GROMACS USA Workshop and Conference Steering Committee
Michael Shirts (chair)
Angel Garcia
Berk Hess
Yu-Shan Lin
Erik Lindahl
Peter Kasson
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] LINCS Constraints - all-bonds or h-bonds?

[Michael Shirts]

I don't go beyond 2 fs with either all- bonds or h-bonds. Things like kinetic 
energy start being subtly off.

H-bonds has less chance of failing with large numbers of constraints- less 
iteration required, especially if bond system cross parallelization boundaries.

If your molecules are < 10 atoms, it probably doesn't matter either way.

Sent from my iPhone

On Aug 15, 2013, at 9:11, "Barnett, James W."  wrote:

> Searching through this mailing list it seems like some have stated that 
> with a 2 fs timestep (dt=0.002), constraints=h-bonds is fine in general.
> 
> The questions I have are:
> 
> 1) What are some personal opinions on when it is ok to switch to h-bonds 
>   from all-bonds for LINCS constraints? Is 2 fs and h-bonds a general 
> practice?
> 
> 2) Also, if typically 2 fs and h-bonds are ok, what time-step do users 
>   (or you personally) generally go to with all-bonds? 
> 
> I am speaking generally here of course. Thanks for your responses.
> 
> -- 
> Wes Barnett | jbarn...@tulane.edu
> 
> -- 
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Unkwown Keyword HILLS

[Michael Shirts]

This is a plumed error, not a gromacs error. Gromacs never handles those 
keywords. 

Sent from my iPhone

On Aug 14, 2013, at 1:40, Albert  wrote:

> Does anybody have any idea what's the problem?
> 
> I use the tutorial example and I don't know why it doesn't work.
> 
> THX
> 
> 
> On 08/13/2013 07:19 PM, Albert wrote:
>> Dear:
>> 
>> I am trying to run plumed with gromacs plugin. Here is my plumed.dat file 
>> which I defined two dihedral angels as cvs:
>> 
>> *HILLS HEIGHT 0.3 W_STRIDE 450
>> WELLTEMPERED SIMTEMP 310 BIASFACTOR 1.96
>> TORSION LIST 1 4 65 344 SIGMA 0.12
>> TORSION LIST 2 46 80 656 SIGMA 0.12
>> 
>> ENDMETA*
>> 
>> I am using plumed-1.3+gromacs-4.6.2 with command:
>> 
>> 
>> *mpirun -np 24 mdrun_mpi -s md.tpr -plumed plumed.dat -g md.log -v -x md.xtc 
>> -o md.trr -e md.edr*
>> 
>> 
>> but it failed with messages:
>> 
>> 
>> *starting mdrun 'protein'
>> 1 steps, 20.0 ps.
>> ! PLUMED ERROR: Line 1 Unkwown Keyword HILLS
>> 
>> ! ABORTING RUN
>> ! PLUMED ERROR: Line 1 Unkwown Keyword HILLS
>> 
>> ! ABORTING RUN
>> --*
>> 
>> thank you very much
>> 
>> Albert
> 
> -- 
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www interface 
> or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Unkwown Keyword HILLS

[Michael Shirts]

This is

Sent from my iPhone

On Aug 14, 2013, at 1:40, Albert  wrote:

> Does anybody have any idea what's the problem?
> 
> I use the tutorial example and I don't know why it doesn't work.
> 
> THX
> 
> 
> On 08/13/2013 07:19 PM, Albert wrote:
>> Dear:
>> 
>> I am trying to run plumed with gromacs plugin. Here is my plumed.dat file 
>> which I defined two dihedral angels as cvs:
>> 
>> *HILLS HEIGHT 0.3 W_STRIDE 450
>> WELLTEMPERED SIMTEMP 310 BIASFACTOR 1.96
>> TORSION LIST 1 4 65 344 SIGMA 0.12
>> TORSION LIST 2 46 80 656 SIGMA 0.12
>> 
>> ENDMETA*
>> 
>> I am using plumed-1.3+gromacs-4.6.2 with command:
>> 
>> 
>> *mpirun -np 24 mdrun_mpi -s md.tpr -plumed plumed.dat -g md.log -v -x md.xtc 
>> -o md.trr -e md.edr*
>> 
>> 
>> but it failed with messages:
>> 
>> 
>> *starting mdrun 'protein'
>> 1 steps, 20.0 ps.
>> ! PLUMED ERROR: Line 1 Unkwown Keyword HILLS
>> 
>> ! ABORTING RUN
>> ! PLUMED ERROR: Line 1 Unkwown Keyword HILLS
>> 
>> ! ABORTING RUN
>> --*
>> 
>> thank you very much
>> 
>> Albert
> 
> -- 
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www interface 
> or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Hamiltonian replica exchange not working in 4.6

[Michael Shirts]

Hi, Sanku-

The way to invoke Hamiltonian replica exchange has changed to be a bit
more flexible.  We should go back and make sure that this legacy way
is supported (I thought this invocation was supported, but apparently
it isn't), but what you should be able to do to get it working quickly
is include in all the .mdp files the line:

fep-lambdas  0 0.08 . . . . . . 1.0

I.e. list all of the lambdas that all of the files use.  This is
described in the documentation.

It should run correctly then.  If not, then bug me on the list again.


On Tue, Aug 6, 2013 at 6:52 PM, Sanku M  wrote:
> Dear Gromacs user,
>   I am not sure whether this is any potential bug. But I found that in 
> gromacs4.6.3, hamiltonian replica exchange simulation is not
> working. The SAME simulation works perfectly fine in gromacs 4.5.4
> But, when I try to run the simulation using gromacs4.6.3, I get following 
> error:
> "The properties of the 10 systems are all the same, there is nothing to 
> exchange
> For more information and tips for troubleshooting, please check the GROMACS"
>
>
> Here is the details on the error:
>I have first defined A and B state by scaling the relevant parameter.. 
> Then I have  generated different tpr file by using mdp file which only differ 
> by 'init-lambda' parameter:
> For example, I have 10 different mdp files: Free energy part of two of them 
> are shown below:
> for hremd0.mdp, the relevant part for free energy part is
>   ; Free energy control stuff
> free_energy  = yes
> init-lambda  = 0
> delta_lambda = 0
>
> for hremd1.mdp
> ; Free energy control stuff
> free_energy  = yes
> init-lambda  = .080
> delta_lambda = 0
> ..
> and so on..
> The rest of mdp portions are same for all the mdp files...
>
>
>
> However, with gromacs4.5.4, the same simulation works without problem. I do 
> not get the error that all the systems are identical.
>
> In both cases, I am using same system with OPLS forcefield..
> The way I run the simulations is:
> mpirun -np 10 mdrun_463 -multi 10 -v -s hremd -replex 250 -c &>log_repl
>
>
> Any help will be appreciated.
>
> Sanku
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Unphysical conformations in decoupled free energy simulation

[Michael Shirts]

Hi, all-

That particular redmine was in the case of couple-intramol NOT being
specified.  So it may be different.

Can you upload samples files for this test case to redmine.gromacs.org
with specific instructions?  Please upload both the 'behaving
correctly' and 'not behaving correctly' versions.

Note that in may cases, you can get the correct behavior by changing
the topologies themselves, though this is more complicated.

On Mon, Aug 5, 2013 at 10:25 AM, Justin Lemkul  wrote:
>
>
> On 8/5/13 9:52 AM, Joerg Sauter wrote:
>>
>> On 08/05/2013 02:15 PM, Justin Lemkul wrote:
>>>
>>>
>>>
>>> On 8/5/13 8:12 AM, Joerg Sauter wrote:

 On 08/05/2013 01:54 PM, Justin Lemkul wrote:
>
>
>
> On 8/5/13 7:44 AM, Joerg Sauter wrote:
>>
>> Dear all,
>>
>> I am trying to compute the free energy of hydration for cellobiose (a
>> beta
>> (1-4)
>> glucose dimer) using BAR in Gromacs 4.6.1. However, I encounter a
>> problem.
>>
>> I find that the vacuum conformations of the molecule in a regular
>> vacuum
>> simulation differ from the conformations in the decoupled simulation
>> in the
>> free
>> energy case i.e., with the additional mdp entries:
>>
>> free-energy  = yes
>> init-lambda = 0
>> couple-lambda0   = none
>> couple-moltype   = solute
>> couple-intramol  = no
>>
>> Here is a histogram of the dihedral angles of the glycosidic linkage
>> in
>> the vacuum case
>> https://dl.dropboxusercontent.com/u/70358077/reg.pdf (it stays in the
>> global
>> free energy minimum)
>> and this is the decoupled free energy case (same starting conformation
>> in the
>> global free energy minimum)
>> https://dl.dropboxusercontent.com/u/70358077/fe.pdf
>>
>> I understand that Gromacs replaces the intramolecular interactions
>> with
>> explicit
>> pair
>> interactions. Therefore, I had to increase table-extension but that
>> did not
>> change much. I was thinking that maybe this could be a problem
>> specific to
>> this
>> topology (a GLYCAM06h conversion from Amber), however, I do not see
>> how
>> this can
>> occur.
>> I hope someone has an idea what is going wrong.
>>
>> mdp file: https://dl.dropboxusercontent.com/u/70358077/md.mdp (same
>> behaviour
>> for longer runs)
>> top file: https://dl.dropboxusercontent.com/u/70358077/cellob.top
>> gro file: https://dl.dropboxusercontent.com/u/70358077/cellob.gro
>>
>> This is sort of a minimal example. The same problem occurs
>> in a regular simulation when decoupling from water using multiple
>> lambdas.
>>
>
> Have you tried running with rlist=rcoulomb=rvdw = 0 rather than 99? If
> nothing
> else, I would try again with version 4.6.3.  There have been tons of
> fixes to
> the free energy code since 4.6.1, though I don't know offhand which, if
> any,
> would be affecting your results. Regardless, the current version should
> always
> be the "best" version ;)
>
> -Justin
>
 Hi Justin,

 I just tried it with version 4.6.3. and rlist=rcoulomb=rvdw = 0.
 Unfortunately,
 I get the same results.

>>>
>>> OK, good to know.  We can rule out 4.6.1-era bugs, then.  Does
>>> visualization
>>> or an analysis of energy terms stored in the .edr file provide any clues
>>> as to
>>> the source of the problem?
>>>
>>> -Justin
>>>
>> Visualization does not show any abnormalities beyond odd confomations in
>> the
>> free energy coupling case.
>>
>> Regarding the energies: I am not sure whether this is normal behaviour,
>> but, in
>> the free energy coupling case I do not see energy values for LJ 14 or
>> Coulomb 14
>> interactions. I thought that by selecting couple-intramol=0 all
>> intramolecular
>> interactions would be replaced by pairtype (14) interactions.
>>
>> In addition, for Coulomb (SR) I get in the free energy case
>> LJ (SR)-22.8084  0.0765.90879 -0.379238
>> (kJ/mol)
>> in comparisson to
>> LJ (SR) 16.8295   0.8113.2303 4.51371
>> (kJ/mol)
>> in the regular vacuum simulation. So this looks like the 14 interactions
>> get
>> shifted to the LJ/Coulomb terms in the coupling simulation?! (However, the
>> numbers do not add up)
>>
>
> This is smelling buggy to me, and there have been issues with 1-4
> interactions in the free energy code before
> (http://redmine.gromacs.org/issues/1225), so there could be a lingering
> problem.  Hopefully Michael will chime in.  It might be worthwhile to open
> up a bug report on redmine.gromacs.org with all of the information you have
> provided (input files, output data, and the very nice description of the
> issue you have posted).
>
>
> -Justin
>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> 

Re: [gmx-users] Re: restraint-lambdas for position restraints in hamiltonian exchange

[Michael Shirts]

That is correct. Such a functionality wouldn't be that hard to
implement - but there are a long list of easy functionalities to be
implemented. You can submit a request to redmine.gromacs.org so that
the request is archived.

On Sat, Aug 3, 2013 at 6:12 PM, Dejun Lin  wrote:
> So I take it that in the position restraint case (not COM-pulling), where
> the reference positions are determined by the starting structure instead of
> a B-state topology, the reference positions won't be swapped ?
>
>
> 2013/8/3 Michael Shirts-2 [via GROMACS] <
> ml-node+s5086n5010324...@n6.nabble.com>
>
>> Short answer is anything that has a B state parameter can be included
>> in in Hamiltonian exchange.
>>
>> If it's pull code or explicit restraints, it's controlled by restraint
>> lambda.
>>
>> > I went through the manual and couldn't find any definite answers to the
>> > following questions.
>> >
>> > First, I wonder if the reference positions of position restraints, not
>> just
>> > the force constants, of different replicas are exchanged in hamiltonian
>> > exchange based on restraint-lambdas? For example, if I have 1 molecule
>> that
>> > has 2 structures, say, A and B and the following 2-component
>> > restraint-lambdas:  1.0 1.0 with replica 0 starting in A and replica 1
>> in B,
>> > would the exchange between replica 0 and 1 yield anything meaningful ?
>>
>> Right now, the pull code doesn't have B state positions, just force
>> constants. It CAN be done using distance restraints, but then you have
>> to make the atoms involved part of the same molecule.  Adding B state
>> distances to pull code was recently requested, and will probably be
>> straightforward to put in soon.
>>
>> > In
>> > other word, would the relaxation towards the other structure occur in
>> both
>> > states once an exchange was accepted?
>> >
>> > Another question is what types of restraints can restraint-lambdas act
>> on?
>> > For example, bond/angle/distance/position restraints, COM-pulling
>> potential
>> > ?
>>
>> Anything with a B state in gromacs topologies.
>> --
>> gmx-users mailing list[hidden 
>> email]<http://user/SendEmail.jtp?type=node&node=5010324&i=0>
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to [hidden 
>> email]<http://user/SendEmail.jtp?type=node&node=5010324&i=1>.
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>>
>> --
>>  If you reply to this email, your message will be added to the discussion
>> below:
>>
>> http://gromacs.5086.x6.nabble.com/restraint-lambdas-for-position-restraints-in-hamiltonian-exchange-tp5010322p5010324.html
>>  To unsubscribe from restraint-lambdas for position restraints in
>> hamiltonian exchange, click 
>> here<http://gromacs.5086.x6.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=5010322&code=ZGVqdW4ubGluQGdtYWlsLmNvbXw1MDEwMzIyfDE0MzAxNDIxNA==>
>> .
>> NAML<http://gromacs.5086.x6.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>>
>
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/restraint-lambdas-for-position-restraints-in-hamiltonian-exchange-tp5010322p5010325.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] restraint-lambdas for position restraints in hamiltonian exchange

[Michael Shirts]

Short answer is anything that has a B state parameter can be included
in in Hamiltonian exchange.

If it's pull code or explicit restraints, it's controlled by restraint lambda.

> I went through the manual and couldn't find any definite answers to the
> following questions.
>
> First, I wonder if the reference positions of position restraints, not just
> the force constants, of different replicas are exchanged in hamiltonian
> exchange based on restraint-lambdas? For example, if I have 1 molecule that
> has 2 structures, say, A and B and the following 2-component
> restraint-lambdas:  1.0 1.0 with replica 0 starting in A and replica 1 in B,
> would the exchange between replica 0 and 1 yield anything meaningful ?

Right now, the pull code doesn't have B state positions, just force
constants. It CAN be done using distance restraints, but then you have
to make the atoms involved part of the same molecule.  Adding B state
distances to pull code was recently requested, and will probably be
straightforward to put in soon.

> In
> other word, would the relaxation towards the other structure occur in both
> states once an exchange was accepted?
>
> Another question is what types of restraints can restraint-lambdas act on?
> For example, bond/angle/distance/position restraints, COM-pulling potential
> ?

Anything with a B state in gromacs topologies.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Expanded ensemble simulation died with fatal error: Something wrong in choosing new lambda state with a Gibbs move

[Michael Shirts]

Hi Dejun-

The basic problem is that for this particular configuration, the
current state is the only state with nonzero weight. Note that the
state with the second highest weight has weight 10^-7.  When it tries
to compare weights in single precision, it has a numerical overflow
and fails.

A few things:

1. This really should be more robust, so that it will realize it's
supposed to stay in the most likely state, since that's the only state
with nonzero weight.  I have a fix that I've been working on for
exactly this problem, but it's not quite ready yet.  Hopefully in the
next couple of days.

2. This problem is very unlikely to occur in double precision, if you
can afford the performance hit in the meantime.

3. If this is a typical average difference, exchanges will be very
unlikely.  You should probably choose your lambda intervals to be a
bit closer together at the end range.

Hopefully this will give you enough information to move forward for
the time being until a better fix is implemented.


On Wed, Jul 31, 2013 at 2:15 PM, Dejun Lin  wrote:
> Hi all,
>
> I'm running an expanded ensemble simulation using gromacs 4.6.3 and it
> crashed with the error:
>
> Fatal error:
> Something wrong in choosing new lambda state with a Gibbs move -- probably
> underflow in weight determination.
> Denominator is:   0 1.002384e+00
>   idEnumerator  weights
>   0 -9.1451739502e+02 0.00e+00 0.00e+00
>   1 -9.128174e+02 0.00e+00 0.00e+00
>   2 -8.8548516846e+02 0.00e+00 0.00e+00
>   3 -8.7096899414e+02 0.00e+00 0.00e+00
>   4 -8.5645288086e+02 0.00e+00 0.00e+00
>   5 -8.4193676758e+02 0.00e+00 0.00e+00
>   6 -8.2742059326e+02 0.00e+00 0.00e+00
>   7 -8.1290447998e+02 0.00e+00 0.00e+00
>   8 -7.9838836670e+02 0.00e+00 0.00e+00
>   9 -7.8387219238e+02 0.00e+00 0.00e+00
>  10 -7.6935607910e+02 0.00e+00 0.00e+00
>  11 -7.5483990479e+02 0.00e+00 0.00e+00
>  12 -7.4032379150e+02 0.00e+00 0.00e+00
>  13 -7.2580767822e+02 0.00e+00 0.00e+00
>  14 -7.1129150391e+02 0.00e+00 0.00e+00
>  15 -6.9677539062e+02 0.00e+00 0.00e+00
>  16 -6.8225927734e+02 0.00e+00 0.00e+00
>  17 -6.6774316406e+02 0.00e+00 0.00e+00
>  18 -6.5322698975e+02 0.00e+00 0.00e+00
>  19 -6.3871087646e+02 0.00e+00 0.00e+00
>  20 -6.2419470215e+02 0.00e+00 0.00e+00
>  21 -6.0967858887e+02 0.00e+00 0.00e+00
>  22 -5.9516247559e+02 0.00e+00 0.00e+00
>  23 -5.8064630127e+02 0.00e+00 0.00e+00
>  24 -5.6613018799e+02 0.00e+00 0.00e+00
>  25 -5.5161407471e+02 0.00e+00 0.00e+00
>  26 -5.3709790039e+02 0.00e+00 0.00e+00
>  27 -5.2258178711e+02 0.00e+00 0.00e+00
>  28 -5.0806564331e+02 0.00e+00 0.00e+00
>  29 -4.9354953003e+02 0.00e+00 0.00e+00
>  30 -4.7903335571e+02 0.00e+00 0.00e+00
>  31 -4.6451724243e+02 0.00e+00 0.00e+00
>  32 -4.518311e+02 0.00e+00 0.00e+00
>  33 -4.3548400879e+02 0.00e+00 0.00e+00
>  34 -4.2096792603e+02 0.00e+00 0.00e+00
>  35 -4.0645178223e+02 0.00e+00 0.00e+00
>  36 -3.9193563843e+02 0.00e+00 0.00e+00
>  37 -3.8107025146e+02 0.00e+00 0.00e+00
>  38 -3.6290338135e+02 0.00e+00 0.00e+00
>  39 -3.4838723755e+02 0.00e+00 0.00e+00
>  40 -3.3387109375e+02 0.00e+00 0.00e+00
>  41 -3.1935494995e+02 0.00e+00 0.00e+00
>  42 -3.0483883667e+02 0.00e+00 0.00e+00
>  43 -2.9032269287e+02 0.00e+00 0.00e+00
>  44 -2.7580654907e+02 0.00e+00 0.00e+00
>  45 -2.6129040527e+02 0.00e+00 0.00e+00
>  46 -2.4677430725e+02 0.00e+00 0.00e+00
>  47 -2.3225816345e+02 0.00e+00 0.00e+00
>  48 -2.1774200439e+02 0.00e+00 0.00e+00
>  49 -2.0322586060e+02 0.00e+00 0.00e+00
>  50 -1.8970976257e+02 0.00e+00-1.00e+00
>  51 -1.7419361877e+02 0.00e+00 0.00e+00
>  52 -1.5967747498e+02 0.00e+00 0.00e+00
>  53 -1.4516131592e+02 0.00e+00 0.00e+00
>  54 -1.3064523315e+02 0.00e+00 0.00e+00
>  55 -1.1612908173e+02 0.00e+00 0.00e+00
>  56 -1.0161293030e+02 7.0064923216e-45 0.00e+00
>  57 -8.7096786499e+01 1.4939846888e-38 0.00e+00
>  58 -7.2580688477e+01 3.0102835162e-32 0.00e+00
>  59 -5.8064540863e+01 6.0658294505e-26 0.00e+00
>  60 -4.3548393250e+01 1.865341e-19 0.00e+00
>  61 -2.9032243729e+01 2.4629560544e-13 0.00e+

[gmx-users] Reminder about US GROMACS workshop + soliciting presenters for talks and tutorials

[Michael Shirts]

Dear gmx-users list members-

I'd like to remind you all about the 2013 GROMACS Workshop and
Conference at the University of Virginia Sept 13th-15th; registration
costs will rise from $60 to $95 on Aug 22nd. The original invitation,
with link to the website, is at the end of this email.

We are issuing a last call for contributed 20 min talks, with
decisions on the schedule made next week, and would like to invite you
to apply to present your GROMACS-related research. Email
gromacsworkshop-registrat...@virginia.edu if you wish to submit a talk
for consideration; see the preliminary schedule on the conference
webpage for more information on the design of the talks.

We still also want to give a opportunity to advanced users who may be
interested in either presenting specific tutorial subjects at the
workshop or alternatively having an online tutorial hosted in the
GROMACS website. See the list of tutorial topics on the workshop
schedule, or propose your own. We want to try to leverage the energy
and experience of the GROMACS community to enrich the tutorial
experience, rather than have the same people up there the whole
weekend!

Sincerely,

The 2013 GROMACS USA Workshop and Conference Steering Committee

Michael Shirts (chair)
Angel Garcia
Berk Hess
Yu-Shan Lin
Erik Lindahl
Peter Kasson

Original announcement:

~~~

We are pleased to announce the 2013 GROMACS USA Workshop and
Conference at the University of Virginia in Charlottesville, Virginia,
on September 13-15th. This will be the first full GROMACS workshop
held here in the U.S. The workshop and conference will consist of two
days of tutorials, discussions of future plans for GROMACS, face time
with a large percentage of the developers, plenary software and
application sessions, and a last day of development planning to which
attendees are also invited to help influence the future directions of
GROMACS.


Please visit http://faculty.virginia.edu/gromacsworkshop for
registration and for much more information about the program.

For specific questions about registering or logistics after visiting
the website, please write to
gromacsworkshop-registrat...@virginia.edu.


Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shi...@virginia.edu
(434)-243-1821
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: Question on Hessian and forces not being preserved on restart with Hessian calculation?

[Michael Shirts]

Problem partly addressed.  If I run normal modes in -nt 1, then I get
the same force as after the minimization.

I'll file a redmine.

On Sat, Jul 20, 2013 at 9:43 AM, Michael Shirts  wrote:
> To follow up -- if I try to minimize again using -t, I get the same
> low forces as in the minimization in the previous step.  So it appears
> to be something with what do_nm is doing, not with errors in the
> output structure.
>
> On Sat, Jul 20, 2013 at 9:37 AM, Michael Shirts  wrote:
>> When I minimize a structure, I can get down to the force max being <0.01
>>
>> Low-Memory BFGS Minimizer converged to Fmax < 0.01 in 6839 steps
>> Potential Energy  = -5.12340607768673e+03
>> Maximum force =  6.68907856457542e-03 on atom 3029
>> Norm of force =  2.19978176343026e-03
>>
>> kJ/nm.  However, when I try to perform a normal mode analysis, it
>> complains that the maximum energy is >30 kJ/mn.
>>
>>
>> Non-cutoff electrostatics used, forcing full Hessian format.
>> Allocating Hessian memory...
>>
>> Maximum force: 3.91984e+01
>> Maximum force probably not small enough to ensure that you are in an
>> energy well. Be aware that negative eigenvalues may occur when the
>> resulting matrix is diagonalized.
>>
>> I've restarted from the binary .trr output (-t option) , and I'm using
>> double precision. Any suggestions as to why the force max is not the
>> asme on restarting?  Only difference is going from l-bfgs to nm.
>>
>> I've seen some older comments on this, but no resolution and nothing recent.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: Question on Hessian and forces not being preserved on restart with Hessian calculation?

[Michael Shirts]

To follow up -- if I try to minimize again using -t, I get the same
low forces as in the minimization in the previous step.  So it appears
to be something with what do_nm is doing, not with errors in the
output structure.

On Sat, Jul 20, 2013 at 9:37 AM, Michael Shirts  wrote:
> When I minimize a structure, I can get down to the force max being <0.01
>
> Low-Memory BFGS Minimizer converged to Fmax < 0.01 in 6839 steps
> Potential Energy  = -5.12340607768673e+03
> Maximum force =  6.68907856457542e-03 on atom 3029
> Norm of force =  2.19978176343026e-03
>
> kJ/nm.  However, when I try to perform a normal mode analysis, it
> complains that the maximum energy is >30 kJ/mn.
>
>
> Non-cutoff electrostatics used, forcing full Hessian format.
> Allocating Hessian memory...
>
> Maximum force: 3.91984e+01
> Maximum force probably not small enough to ensure that you are in an
> energy well. Be aware that negative eigenvalues may occur when the
> resulting matrix is diagonalized.
>
> I've restarted from the binary .trr output (-t option) , and I'm using
> double precision. Any suggestions as to why the force max is not the
> asme on restarting?  Only difference is going from l-bfgs to nm.
>
> I've seen some older comments on this, but no resolution and nothing recent.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Question on Hessian and forces not being preserved on restart with Hessian calculation?

[Michael Shirts]

When I minimize a structure, I can get down to the force max being <0.01

Low-Memory BFGS Minimizer converged to Fmax < 0.01 in 6839 steps
Potential Energy  = -5.12340607768673e+03
Maximum force =  6.68907856457542e-03 on atom 3029
Norm of force =  2.19978176343026e-03

kJ/nm.  However, when I try to perform a normal mode analysis, it
complains that the maximum energy is >30 kJ/mn.


Non-cutoff electrostatics used, forcing full Hessian format.
Allocating Hessian memory...

Maximum force: 3.91984e+01
Maximum force probably not small enough to ensure that you are in an
energy well. Be aware that negative eigenvalues may occur when the
resulting matrix is diagonalized.

I've restarted from the binary .trr output (-t option) , and I'm using
double precision. Any suggestions as to why the force max is not the
asme on restarting?  Only difference is going from l-bfgs to nm.

I've seen some older comments on this, but no resolution and nothing recent.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] segfault with an otherwise stable system when I turn on FEP (complete decoupling)

[Michael Shirts]

Chris, can you post a redmine on this so I can look at the files?

Also, does it crash immediately, or after a while?

On Thu, Jul 18, 2013 at 2:45 PM, Christopher Neale
 wrote:
> Dear Users:
>
> I have a system with water and a drug (54 total atoms; 27 heavy atoms). The 
> system is stable when I simulate it for 1 ns. However, Once I add the 
> following options to the .mdp file, the run dies after a few ps with a 
> segfault.
>
> free-energy = yes
> init-lambda = 1
> couple-lambda0 = vdw-q
> couple-lambda1 = none
> couple-intramol = no
> couple-moltype = drug
>
> I do not get any step files or any lincs warnings. If I look at the .xtc and 
> .edr files, there is no indication of something blowing up before the 
> segfault. I have also verified that the drug runs without any problem in 
> vacuum. I get the same behaviour if I remove constraints and use a timestep 
> of 0.5 fs. The segfault is reproducible with v4.6.1 and v4.6.3. I am using 
> the charmm FF, but I converted all UB angles in my drug to type-1 angles and 
> still got the segfault. I also get the segfault with particle decomposition 
> and/or while running a single thread. I am currently using the SD integrator, 
> but I get the same segfault with md and md-vv. Couple-intramol=yes doesn't 
> resolve it, neither does using separate T-coupling groups for the water and 
> drug. Neither does turning off pressure coupling.
>
> Here is the .mdp file that works fine, but gives me a segfault when I add the 
> free energy stuff (above):
>
> constraints = all-bonds
> lincs-iter =  1
> lincs-order =  6
> constraint_algorithm =  lincs
> integrator = sd
> dt = 0.002
> tinit = 0
> nsteps = 10
> nstcomm = 1
> nstxout = 0
> nstvout = 0
> nstfout = 0
> nstxtcout = 500
> nstenergy = 500
> nstlist = 10
> nstlog=0 ; reduce log file size
> ns_type = grid
> vdwtype = cut-off
> rlist = 0.8
> rvdw = 0.8
> rcoulomb = 0.8
> coulombtype = cut-off
> tc_grps =  System
> tau_t   =  1.0
> ld_seed =  -1
> ref_t = 310
> gen_temp = 310
> gen_vel = yes
> unconstrained_start = no
> gen_seed = -1
> Pcoupl = berendsen
> pcoupltype = isotropic
> tau_p = 4
> compressibility = 4.5e-5
> ref_p = 1.0
>
> I do realize that some of these settings are not ideal for a production run. 
> I started with the real Charmm cutoffs + PME, etc, (which also gives the 
> segfault) but this is what I am using right now for quick testing.
>
> The only thing keeping me from filing a redmine issue is that if I remove my 
> drug and do the FEP on one of the water molecules (using the FEP code listed 
> above), I have no segfault. Therefore it is clearly related to the drug, 
> whose parameters I built so I may have caused the problem somehow. 
> Nevertheless, the drug runs fine in water and in vacuum without the FEP code, 
> so I can't imagine what could be causing this segfault (also, the fact that 
> it's a segfault means that I don;t get any useful info from mdrun as to what 
> might be going wrong).
>
> Thank you,
> Chris.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: gmx 4.6.1, Expanded ensemble: weird balancing factors

[Michael Shirts]

Unfortunately, the semiiso fix for MTTK is not trivial, and will have
to wait until 5.0.

Berendsen volume + expanded ensemble would likely give very bad
results, especially since surface area fluctuations are important.

Replica exchange works with sd and parrinello-Rahman, if that helps!

I am working on checking various combinations of expanded ensemble and
integrators that can more easily be supported, but that may take a few
weeks.







On Wed, Jul 17, 2013 at 6:22 PM, Dejun Lin  wrote:
> The test I did was using gmx-4.6.3. I'm currently working on a membrane
> system, which requires semiisotropic pressure coupling. I know that MTTK in
> gmx-4.6.3 doesn't support semiiso yet so the only combination available for
> expanded ensemble is md-vv + v-rescale tcoupl + berendsen pcoupl, which I'm
> testing right now. A preliminary run gave me some sane dG values although
> it's not optimal for data collection as you pointed out. I wonder if it's a
> trivial fix that might have been done to add semiiso to MTTK?
>
> Thanks again for your help!
>
>
> 2013/7/17 Michael Shirts 
>
>> > It seems to have something to do with the integrator/pressure-coupling.
>>
>> that is what I expected based on some preliminary testing earlier.
>>
>> When
>> > I ran the tutorial on
>> >
>> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Ethanol_solvation_with_expanded_ensemble
>> ,
>> > everything seems fine
>>
>> OK good!
>>
>> > unless I switched to sd integrator instead of md-vv
>>
>> There may be a problem with sd and expanded ensemble. This MAY have
>> been fixed in 4.6.3, but I'll need a little bit of time to check.
>>
>> > (and Pcoupl from MTTK to berendsen),
>>
>> Berendsen is provably wrong. The ensemble is incorrect.  It is evil
>> whenever a true distribution is needed. There are warnings that should
>> be printed when you set up the system.  Perhaps they should be made
>> stronger with expanded ensemble!
>>
>> in which case SHAKE didn't seem to be
>> > able to settle the constraints.
>>
>> Expanded ensemble is rough on constraints.  The integration is
>> formally correct, but it is harder to converge pressure, especially
>> with constraints.  The larger the system, however, the more stable it
>> is.
>>
>> > file from the tutorial. Any idea?
>>
>> Well, one solution is to run your molecule with the parameters in the
>> .mdp!  Would that solve your problem for now?
>>
>> Longer term. I will check on the SD to see if it is fixed in 4.6.3 --
>> I seem to re call this, but I'm not sure.   However, I am not sure
>> that it can be made formally correct.  I will either fix this or make
>> the warnings more clear -- or not permit it at all.
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: gmx 4.6.1, Expanded ensemble: weird balancing factors

[Michael Shirts]

> It seems to have something to do with the integrator/pressure-coupling.

that is what I expected based on some preliminary testing earlier.

When
> I ran the tutorial on
> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Ethanol_solvation_with_expanded_ensemble,
> everything seems fine

OK good!

> unless I switched to sd integrator instead of md-vv

There may be a problem with sd and expanded ensemble. This MAY have
been fixed in 4.6.3, but I'll need a little bit of time to check.

> (and Pcoupl from MTTK to berendsen),

Berendsen is provably wrong. The ensemble is incorrect.  It is evil
whenever a true distribution is needed. There are warnings that should
be printed when you set up the system.  Perhaps they should be made
stronger with expanded ensemble!

in which case SHAKE didn't seem to be
> able to settle the constraints.

Expanded ensemble is rough on constraints.  The integration is
formally correct, but it is harder to converge pressure, especially
with constraints.  The larger the system, however, the more stable it
is.

> file from the tutorial. Any idea?

Well, one solution is to run your molecule with the parameters in the
.mdp!  Would that solve your problem for now?

Longer term. I will check on the SD to see if it is fixed in 4.6.3 --
I seem to re call this, but I'm not sure.   However, I am not sure
that it can be made formally correct.  I will either fix this or make
the warnings more clear -- or not permit it at all.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: window exchange umbrella sampling

[Michael Shirts]

4.5.7 does not support Hamiltonian exchange.  It says all properties
are the same, because all the temperatures and pressures are the same
-- it won't switch the umbrellas.

On Wed, Jul 17, 2013 at 3:30 PM, Parisa  wrote:
> Hi Michael,
>
> I think that this is an issue with the gromacs version I am using (4.6.1). I
> switched to 4.5.7 and now, I get the error:
>
> Fatal error:
> The properties of the ... systems are all the same, there is nothing to
> exchange
>
> I don't understand what I am missing here. I suppose if I have 3 replicas to
> fix a molecule at positions 1.0 and 2.0 and 3.0 while the force constant is
> being exchanged between them in a range of 300 to 200, I can define 3 mdp
> files as below:
>
> For replica 1 at position 1.0:
>
> pull_group0  = POPC
> pull_weights0=
> pull_pbcatom0= 0
> pull_group1  = DNA
> pull_weights1=
> pull_pbcatom1= 0
> pull_vec1= 0.0 0.0 0.0
> pull_init1   = 0.0 0.0 1.0
> pull_rate1   = 0
> pull_k1  = 300
> pull_kB1 = 200
> free_energy  = yes
> restraint-lambdas  =  1.0 0.5 0.0
> init-lambda-state   =  0
> nstdhdl = 10
>
> For replica 2 at position 2.0:
>
> pull_group0  = POPC
> pull_weights0=
> pull_pbcatom0= 0
> pull_group1  = DNA
> pull_weights1=
> pull_pbcatom1= 0
> pull_vec1= 0.0 0.0 0.0
> pull_init1   = 0.0 0.0 2.0
> pull_rate1   = 0
> pull_k1  = 300
> pull_kB1 = 200
> free_energy  = yes
> restraint-lambdas  =  1.0 0.5 0.0
> init-lambda-state   =  1
> nstdhdl = 10
>
> For replica 3 at position 3.0:
>
> pull_group0  = POPC
> pull_weights0=
> pull_pbcatom0= 0
> pull_group1  = DNA
> pull_weights1=
> pull_pbcatom1= 0
> pull_vec1= 0.0 0.0 0.0
> pull_init1   = 0.0 0.0 3.0
> pull_rate1   = 0
> pull_k1  = 300
> pull_kB1 = 200
> free_energy  = yes
> restraint-lambdas  =  1.0 0.5 0.0
> init-lambda-state   =  2
> nstdhdl = 10
>
> Could you please comment if I am missing something here? Unfortunately,
> there is not much about this posted before.
>
> Many thanks,
>
> Parisa
>
>
>
>
>
>
>
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/window-exchange-umbrella-sampling-tp5009894p5009947.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: window exchange umbrella sampling

[Michael Shirts]

Ah, this is a force field issue -- urey-bradley terms are not
supported free energy calculations. However, since only restraints are
changing, this warning doesn't really need to be there.

It would be relatively simple to put in a check to allow this to work,
but it might take a week or two to get around to it.

I'll file a redmine in the meantime, in case someone else can get to it first!

http://redmine.gromacs.org/issues/1302



On Tue, Jul 16, 2013 at 5:32 PM, Parisa  wrote:
> Thanks for quick reply. I fixed the issue with kB1. Below is what I have
> changed in the .mdp file. I use this to run it "mpirun -np 3  mdrun_mpi -s
> prefix_.tpr -multi 3 -replex 6 -dhdl a.dhdl". But I get this error (I am
> using version 4.6.1):
>
> Fatal error:
> Function type U-B not implemented in ip_pert
> For more information and tips for troubleshooting, please check the GROMACS
>
> Here is the input file:
> .
> .
> .
> pull_group0  = POPC
> pull_weights0=
> pull_pbcatom0= 0
> pull_group1  = DNA
> pull_weights1=
> pull_pbcatom1= 0
> pull_vec1= 0.0 0.0 0.0
> pull_init1   = 0.0 0.0 0.0
> pull_rate1   = 0
> pull_k1  = 300
> pull_kB1 = 200
>
> free_energy  = yes
> restraint-lambdas  =  1.0 0.5 0.0
> init-lambda-state   =  0
> nstdhdl = 10
> .
> .
> .
>
> Any idea about it?
>
> Thanks,
>
> Parisa
>
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/window-exchange-umbrella-sampling-tp5009894p5009899.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: gmx 4.6.1, Expanded ensemble: weird balancing factors

[Michael Shirts]

That's a good question.  My understanding and experience is the
difference in required overlap in replica exchange and expanded
ensemble methods is not significant.  There are probably some more
detailed studies out there. Expanded ensemble can be somewhat more
efficient (see a paper by Park and Pande) But you really do want to
have at least 30% of the states moving each time in any case.

However, you DO need less overlap to get fairly good free energy
estimates.  I.e. the overlap for BAR to work fairly well is lower than
the what you need for good mixing.  That's more anecdotal -- I'd have
to dig a bit to get you good results on that . . .

On Tue, Jul 16, 2013 at 3:19 PM, Dejun Lin  wrote:
> Hi Michael,
>
> Thanks for the reply. Just a quick follow-up. Do you think the overlap of
> energy histogram between different lambdas matter for lambda-dynamics in
> general? (Maybe not in this particular case or in the case of the tutorial
> you just posted.) I wonder if we need as many intermediate lambdas as in
> the case of replica exchange since the weights compensate the difference in
> energy.
>
> Thanks,
> Dejun
>
>
> 2013/7/16 Michael Shirts 
>
>> Hi, all-
>>
>> This not a problem with W-L, but is instead something that is wrong
>> with a particular combination of mdp options that are not working for
>> expanded ensemble simulations.  W-L can equilibrate to incorrect
>> distributions because it decreases the weights too fast (more on that
>> later), but that is not the problem here. The option wl-oneovert
>> allows switching to a slower scaling that is more likely to converge.
>>
>> The reason that it is not a W-L problem is that after the WL weights
>> are equilibrated, it is equilibrium sampling in the expanded ensemble.
>>  If the W-L equilibration was a problem, then the distribution of
>> states would not be flat.  They are flat, so therefore the expanded
>> ensemble dynamics are wrong.
>>
>> I have example expanded ensemble setups.
>>
>>
>> http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Ethanol_solvation_with_expanded_ensemble
>>
>> These mdp files should work. Note that you should be able to swap out
>> any top, and it will still work.  If you get CORRECT results (and it
>> should take just a few ns to see) with these tops, then I will go
>> through and try to figure out which differences are causing the
>> problems.
>>
>> Thanks for the patience while we test this new functionality.
>> Frequently when one puts new code in the hands of others it breaks in
>> ways that were not foreseen!
>>
>> On Tue, Jul 16, 2013 at 12:53 PM, Dejun Lin  wrote:
>> > I have a similar question here. Can anyone tell if the Wang-Landau
>> algorithm
>> > in lambda space is robust in that it doesn't depend on the choice of the
>> > convergence factor (weight-equil-wl-delta), the flatness of the histogram
>> > (wl-ratio) and/or the frequency of trial move (nstexpanded), especially
>> in
>> > the case of barely overlapping energy distribution corresponding to
>> > different lambdas? I can imagine in the extreme case, with only 2
>> lambdas (
>> > l = 1 or 0), the gap between the 2 energy distributions might be so large
>> > such that for most of the time, trial moves from the "center" of
>> > distribution 1 would frequently end up in the "tail" of distribution 0.
>> In
>> > this case, the Wang-Landau weights would be biased if there's not enough
>> > time for the system to relax back to equilibrium.
>> >
>> > Thanks,
>> > Dejun
>> >
>> >
>> >
>> > --
>> > View this message in context:
>> http://gromacs.5086.x6.nabble.com/gmx-4-6-1-Expanded-ensemble-weird-balancing-factors-tp5007681p5009892.html
>> > Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
>> > --
>> > gmx-users mailing listgmx-users@gromacs.org
>> > http://lists.gromacs.org/mailman/listinfo/gmx-users
>> > * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> > * Please don't post (un)subscribe requests to the list. Use the
>> > www interface or send it to gmx-users-requ...@gromacs.org.
>> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't po

Re: [gmx-users] window exchange umbrella sampling

[Michael Shirts]

You need to have different pull parameters at the end states. Right
now, pull-kB1 is not defined in your code, so there is nothing to
interpolate to: it assume pull-kB1 = pull-k1.

Longer scale -- one would want to define reference distances that
change with lambda within the same simulation, but looking over the
code, I'm not seeing anything.  One would want to add a pull_vect1B to
change the location of the restraint as a function of lambda, but this
isn't yet defined.  That should be something we look at in the future
. . .

On Tue, Jul 16, 2013 at 2:41 PM, Parisa Akhshi  wrote:
> I am trying to use Hamiltonian replica exchange as implemented in Gromacs.
> The idea is to use different pulling umbrella forces and then obtain a PMF
> profile. I am specifically interested in exchanging force constants between
> windows.
>
> To do a quick test, I successfully ran replica exchange in which the
> temperature is exchanged. In this case, I prepared for example 8 tpr files
> at 8 different temperatures and run them using
>
> mpirun -np 8  mdrun_mpi -s prefix_.tpr -multi 8 -replex 100
>
> And it ran just fine.
>
> For window exchange umbrella sampling, though, I am confused with the input
> preparation. I am not sure how to specify in the .mdp file how the values
> of pull force should be used to define different force constants for
> different replicas. Also, when submitting the job, how should I ask the
> program to exchange pull forces between windows? I tried using  this:
>
> http://lists.gromacs.org/pipermail/gmx-users/2011-April/060448.html
>
> But, it gives me the error:
>
> Fatal error:
> The properties of the ... systems are all the same, there is nothing to
> exchange
>
> I noticed a previous discussion on below link:
>
> http://gromacs.5086.x6.nabble.com/Hamiltonian-replica-exchange-umbrella-sampling-with-gmx-4-6-td5006187.html
>
> However, I am not sure how to use the restraint-lambda parameter exactly.
> Is there any example, ... how to use it?
>
> I have copied the related portion of one of the the .mdp files below for
> moving MOL2 with respect to MOL1
>
> Thanks for your help,
>
> Parisa
>
>
>
>
> .
> .
> .
> pull_group0  = MOL1
> pull_weights0=
> pull_pbcatom0= 0
> pull_group1  = MOL2
> pull_weights1=
> pull_pbcatom1= 0
> pull_vec1= 0.0 0.0 0.0
> pull_init1   = 0.0 0.0 0.0
> pull_rate1   = 0
> pull_k1  = 300
>
> free_energy  = yes
> restraint-lambdas  =  0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0
> .
> .
> .
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: gmx 4.6.1, Expanded ensemble: weird balancing factors

[Michael Shirts]

Hi, all-

This not a problem with W-L, but is instead something that is wrong
with a particular combination of mdp options that are not working for
expanded ensemble simulations.  W-L can equilibrate to incorrect
distributions because it decreases the weights too fast (more on that
later), but that is not the problem here. The option wl-oneovert
allows switching to a slower scaling that is more likely to converge.

The reason that it is not a W-L problem is that after the WL weights
are equilibrated, it is equilibrium sampling in the expanded ensemble.
 If the W-L equilibration was a problem, then the distribution of
states would not be flat.  They are flat, so therefore the expanded
ensemble dynamics are wrong.

I have example expanded ensemble setups.

http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Ethanol_solvation_with_expanded_ensemble

These mdp files should work. Note that you should be able to swap out
any top, and it will still work.  If you get CORRECT results (and it
should take just a few ns to see) with these tops, then I will go
through and try to figure out which differences are causing the
problems.

Thanks for the patience while we test this new functionality.
Frequently when one puts new code in the hands of others it breaks in
ways that were not foreseen!

On Tue, Jul 16, 2013 at 12:53 PM, Dejun Lin  wrote:
> I have a similar question here. Can anyone tell if the Wang-Landau algorithm
> in lambda space is robust in that it doesn't depend on the choice of the
> convergence factor (weight-equil-wl-delta), the flatness of the histogram
> (wl-ratio) and/or the frequency of trial move (nstexpanded), especially in
> the case of barely overlapping energy distribution corresponding to
> different lambdas? I can imagine in the extreme case, with only 2 lambdas (
> l = 1 or 0), the gap between the 2 energy distributions might be so large
> such that for most of the time, trial moves from the "center" of
> distribution 1 would frequently end up in the "tail" of distribution 0. In
> this case, the Wang-Landau weights would be biased if there's not enough
> time for the system to relax back to equilibrium.
>
> Thanks,
> Dejun
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/gmx-4-6-1-Expanded-ensemble-weird-balancing-factors-tp5007681p5009892.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Larger number of decimal places for coordinates with velocities

[Michael Shirts]

Have you checked out the -ndec option for trjconv?  If you have a high
precision format (.trr, or .xtc if they are stored with sufficient
precision) you can print out a .gro file (that gromacs can read) with
higher precision.

Gromacs can read .gro files with increased precisions in the
coordinates -- as long as the precision does not CHANGE anywhere in
the file.

On Fri, Jul 5, 2013 at 11:16 AM, C.M.Sampson  wrote:
> Dear all,
>
> I use Gromacs 4.5.5 with the AMBER ff99SB force field.
>
> I'm working on a method that requires me to run short NVE simulations
> one after the other, but randomly generating velocities at the start of
> each simulation. i.e. I would run 10 ps, use the final structure with
> new random velocities and run another 10 ps.
>
> I had issues with the final potential energy of the previous simulation
> not matching the first potential energy of the current simulation, but
> managed to fix that by converting the coordinates from the .xtc file to
> a .gro file.
>
> My problem is that when I then put the velocities into my new .gro file
> the simulation breaks after the first step and the kinetic energy is way
> too high:
>
> "
>Step   Time Lambda
>   00.00.0
>
>Energies (kJ/mol)
>  Bond  Angle Ryckaert-Bell.  LJ-14 Coulomb-14
> 5.92468e+031.01216e+037.68471e-016.38163e-013.95610e+00
>  LJ (SR)   Coulomb (SR)   Coul. recip.  PotentialKinetic En.
>  8.66770e+03   -2.27543e+035.12418e+046.45763e+042.11278e+07
>Total EnergyTemperature Pressure (bar)
> 2.11924e+077.81749e+051.08075e+07
> "
>
> The temperature should be 300K and if I use a .gro file with less
> decimal places these same velocities work.
>
> I thought it had to be the way I had formatted my .gro file, but found
> the following:
>
> "
> This format is fixed, ie. all columns are in a fixed position.
> Optionally (for now only yet with trjconv) you can write gro files with
> any number of decimal places, the format will then be n+5 positions with
> n decimal places (n+1 for velocities) in stead of 8 with 3 (with 4 for
> velocities). Upon reading, the precision will be inferred from the
> distance between the decimal points (which will be n+5). Columns contain
> the following information (from left to right):
>   * residue number (5 positions, integer)
>   * residue name (5 characters)
>   * atom name (5 characters)
>   * atom number (5 positions, integer)
>   * position (in nm, x y z in 3 columns, each 8 positions with 3
> decimal places)
>   * velocity (in nm/ps (or km/s), x y z in 3 columns, each 8
> positions with 4 decimal places)
> "
>
> within http://manual.gromacs.org/online/gro.html .
>
> An extract of my .gro file can be seen below:
>
> "
> Generated by trjconv : 2168 system t=  15.0
>  2168
> 1ETH C11   2.735383   2.672010   1.450194  0.2345 -0.1622
> 0.2097
> 1ETHH112   0.015804   2.716597   1.460588  0.8528 -0.7984
> 0.6605
> 1ETHH123   2.744822   2.565544   1.409227 -2.3812  2.8618
> 1.8101
> "
>
> I was wondering if there's a way to use a larger number of decimal
> places for the coordinates with velocities?
>
> Best Wishes
>
> Chris Sampson
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] a question concerning on entropy

[Michael Shirts]

No.

This is a statistical mechanical issue, not a GROMACS issue.  For
interacting systems, entropy is a quantity describing the system as a
whole, and cannot be defined for different parts of the system, at
least not in any way such that the individual components can be added
together.

I'm also not certain the md.edr file will give the entropy of the system . . .

On Thu, Jul 4, 2013 at 2:28 PM, Albert  wrote:
> Hello :
>
>  I've got a question about the the entropy. As we all know that in the
> md.edr file it will give us the entropy value of the system along the
> simulations.
>
> However, my system is a protein/membrane system, and I am only would like to
> make statics for the protein/water related entropy. I am just wondering, is
> it possible to do this?
>
> thank you very much
>
> Albert
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: 1-4 interactions free energy calculations

[Michael Shirts]

Hi, Sonia-

Gromacs 4.6.2 (some bug fixes vs 4.6.1) with the example files you
point out from Alchemistry.org should work well for expanded ensemble.
 David Mobley and I have been validating expanded ensemble and replica
exchange, and the files posted there now are stable for all sizes of
systems, and give reasonable results.  expanded ensemble with NPT is
somewhat unstable for smaller systems (1000 total molecules), though
is much more stable for larger systems -- the posted parameter should
work.

However, note that you have to be careful with expanded ensemble, as
the weights can sometime converge to a poor value.

I'm working on some more guidance on expanded ensemble besides the
examples on Alchemistry.org, but it may take a few weeks.  Check back
in a few weeks for some more information on Alchemistry.org.

On Tue, Jun 25, 2013 at 4:36 PM, Sonia Aguilera
 wrote:
> Hi Justin,
>
> Thank you for your answer. I’m performing several tests to see what is the
> best for my system. The reason I decided not to couple the intramolecular
> interactions is because I think that the annihilation of the molecule will
> lead to very extreme configurations and that will affect my phase-space
> overlap. I know I can just increase the number of intermediate states, but I
> wanted to first test my theory. Also, I have limited resources to run.
>
> There is a gap between the indirect and direct calculation methods  for my
> system. If I set up couple-intramol=no, then I can´t use domain
> decomposition but particle decomposition in all simulations (charge groups
> are too big, and there is no domain decomposition… ) and I got the 1-4
> interactions warning. If I restrain the h-bonds, then the minimization
> converges but not to Fmax less than 100 in only 499 steps and 12 minutes (8
> cores machines). If I continue and run the NVT stabilization, then it takes
> 17h20min to complete only 100 ps, and I still got the 1-4 interactions
> warning. I still have to run the NPT and 3 ns MD, which will take too much
> time. I think it worth it  because I think it will help with the phase-space
> overlap, and I will only have to perform 20 series of simulations (vdw and
> coulomb separately) to get the free energy change.
>
> On the other hand, if I set up couple-intramol=yes, then I have to run 20
> extra simulations in vacuo to counter the annihilation effect. I think that
> it will also affect the phase-space overlap. The good thing is that I can
> use domain decomposition and the minimization converges to machine precision
> (takes 7000 steps and around 40 minutes in a machine with the same
> characteristics that the previous one). Also, the same 100-ps NVT takes only
> 3h20min. If I think of the overall picture, it seems that the total time I
> will spend doing dd (coupleintramol=yes) rather than pd (coupleintramol=no)
> will be less. However I’m worried about the phase-space overlap and the
> available resources I have for running since it makes difficult to just
> increase the intermediate states.
>
> Also, I think I can try doing the same but with the 4.6.1 version. It looks
> like I can now couple both vdw and coulomb in the same simulation
> (http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy).
> I have also wanted to try expanded ensemble to improve my sampling, but I
> have not found that much documentation about the implementation on gromacs.
> It would be great if you can provide an example.
>
> So, do you know another way to improve my sampling and the phase-space
> overlap that can help me to solve my problem?  Also, I’m working with a
> charged molecule. Would it help in any aspect if I neutralize the system?
> Thank you very much in advance,
>
> Best regards,
>
> Sonia Aguilera
> Graduate Assistant
> Department of Chemical Engineering
> Universidad de los Andes
>
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/1-4-interactions-free-energy-calculations-tp5009364p5009378.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Enthalpy Confusion

[Michael Shirts]

> or should i be doing < U+*ref_p > = ?

More specifically,  + *ref_p = H

 isn't really meaningful thing.  I mean, you can define something
such that  = H, but that's not really thermodynamics.

> example system gives  = -1168 kJ/mol and i find  = -725 kJ/mol either

Interesting.  What material at what phase conditions?  For liquids,
the PV contribution should be very small.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] 2013 GROMACS USA Workshop and Conference

[Michael Shirts]

Dear GROMACS users list members:

We are pleased to announce the 2013 GROMACS USA Workshop and
Conference at the University of Virginia in Charlottesville, Virginia,
on September 13-15th.  This will be the first full GROMACS workshop
held here in the U.S.

The workshop and conference will consist of two days of tutorials,
discussions of future plans for GROMACS, face time with a large
percentage of the developers, plenary software and application
sessions, and a last day of development planning to which attendees
are also invited to help influence the future directions of GROMACS.

Please visit http://faculty.virginia.edu/gromacsworkshop for
registration and for much more information about the program.

For specific questions about registering or logistics after visiting
the website, please write to
gromacsworkshop-registrat...@virginia.edu.

Sincerely,
The 2013 GROMACS USA Workshop and Conference Steering Committee
Michael Shirts (chair)
Angel Garcia
Berk Hess
Yu-Shan Lin
Erik Lindahl
Peter Kasson
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Enthalpy Confusion

[Michael Shirts]

If you are computing enthaply in the NPT ensemble, P is constant, and
is the applied pressure.

The "pressure" quantity calculated from the KE and the virial is not
the pressure.  It is a quantity that when averaged over time is equal
the pressure.  Only the average is meaningful macroscopically.

If you are computing enthalpy in another ensemble (which is possible,
though it may be harder to interpret) then you would use the average
pressure from this quantity.



On Tue, Jun 11, 2013 at 3:08 PM, Jeffery Perkins  wrote:
> that's what i thought, and what i tried to do, my pressure is a bit higher
> then that, we want a Lennard-Jones liquid so it's running at 1000+ bar, and
> while I agree that gromacs is giving H as Etot + pV it appears that when i
> calculate pV i get a different value from what g_energy returns for it I did
> p in Pa, V in m^3, as the whole box, to get J of energy, and then multiply
> by 6.02E23 to get J/mol of my box, and then convert down to kJ/mol to be the
> same units as g_energy.
>
> When you say the applied pressure you mean that p = ref_p? or the calculated
> pressure from the virial and Ke?
>
> Thanks for the help,
>
> Jeffery
>
>
>
> --
> View this message in context: 
> http://gromacs.5086.x6.nabble.com/Enthalpy-Confusion-tp5009053p5009058.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Free Energy Calculations in Gromacs

[Michael Shirts]

An important final point is that you can always see EXACTLY what
grompp is putting into the B state by running gmxdump on the resulting
tpr.  It's a LOT of information, but all in text all the interactions
are listed explicitly there.

On Mon, Jun 10, 2013 at 6:20 PM, Justin Lemkul  wrote:
>
>
> On 6/10/13 2:50 PM, JW Gibbs wrote:
>>
>> Hi,
>>
>> I have been trying to perform the simulations using the amber forcefield,
>> in
>> which the [ pairtypes ] directive is not defined explicitly in the
>> ffnonbonded.itp file, but rather the 1-4 interactions are generated as per
>> the defaults section in the forcefield.itp file. I was wondering what
>> happens to these 1-4 interactions at lambda=1 state?
>>
>> Suppose 1 corresponds to CA and 4 corresponds to NA at state A (lambda =
>> 0)
>> and the CA_per and NA_per are the corresponding atomtypes at state B
>> (lambda
>> = 1).
>>
>> It is defined in the topology file that
>>
>> ...
>> ...
>> [ pairs ]
>>
>> 1  4   1
>> 
>> 
>>
>> So does it mean that at state B (lambda = 1), the 1-4 nonbonded
>> interactions
>> will be calculated between CA_per and  NA_per?
>>
>> Although the [ pairs_nb ] section is described in the manual, I would
>> appreciate if someone elaborates a little more.
>>
>
> The actual answer depends on exactly how you're treating the system.  Are
> you using [pairs_nb] in your topology or just [pairs]?  The outcome will be
> different depending on which you're using.  Also note that, as the manual
> says, you don't really need to mess with [pairs_nb] at all; you can achieve
> unscaled internal interactions using "couple-intramol = no" in the .mdp
> file.  Without seeing the .mdp file, it's even more difficult to know what
> you're doing and what you should expect.
>
> -Justin
>
> --
> 
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Difference between the electrostatic treatments PME/Cut-offs and Reaction Field

[Michael Shirts]

> It should also be noted (and obvious now that I actually look into it) that 
> using dispersion correction results in both the latent heat of vapourisation 
> and density of the alkanes being over estimated (for both Cut-off and 
> Reaction Field, and by the same amount).

That may not be quite the right way to think about it.  The dispersion
correction gives a result that (for homogeneous fluids) is essentially
cutoff independent.  So it's a little "truer" calculation in that it
is relatively independent of your potential cutoff, which is a
nonphysical quantity that's really an artifact of the simulation.

So saying that the properties are overestimated with a dispersion
correction, you're essentially saying that the latent heat of
vaporization and the densities would also be overestimated if your ran
with longer cutoffs, and that they are only right if run with quite
short cutoffs.

On Tue, Jun 4, 2013 at 10:51 PM, Dallas Warren  wrote:
> "Problem" solved.
>
> Just as well I held off reporting the "issue" in full until I had explored 
> everything, I would ended have look a bit stupid ;-)  But asking the initial 
> question helped direct my thinking to the reason.
>
> The issue was the difference in van der Waals cut off between the PME/Cut-off 
> and Reaction Field methods that I was using, 0.9/0.9/0.9 vs 0.8/1.4/1.4  The 
> difference was hidden in the results until you turned off Dispersion 
> Correction, which was what was confusing me and final led to realizing what 
> is going on.  My forehead hurt from slapping it after coming to the 
> realization ...
>
> It should also be noted (and obvious now that I actually look into it) that 
> using dispersion correction results in both the latent heat of vapourisation 
> and density of the alkanes being over estimated (for both Cut-off and 
> Reaction Field, and by the same amount).
>
> Catch ya,
>
> Dr. Dallas Warren
> Drug Discovery Biology
> Monash Institute of Pharmaceutical Sciences, Monash University
> 381 Royal Parade, Parkville VIC 3052
> dallas.war...@monash.edu
> +61 3 9903 9304
> -
> When the only tool you own is a hammer, every problem begins to resemble a 
> nail.
>
>
>> -Original Message-
>> From: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> boun...@gromacs.org] On Behalf Of Mark Abraham
>> Sent: Thursday, 30 May 2013 9:32 PM
>> To: Vitaly Chaban; Discussion list for GROMACS users
>> Subject: Re: [gmx-users] Difference between the electrostatic
>> treatments PME/Cut-offs and Reaction Field
>>
>> Things should be identical - any quantity computed from a zero charge
>> has
>> to be zero :-).
>>
>> Mark
>>
>> On Thu, May 30, 2013 at 1:26 PM, Dr. Vitaly Chaban
>> wrote:
>>
>> > Hmmm...
>> >
>> > And what does the observed difference look like, numerically?
>> >
>> >
>> >
>> >
>> >
>> > On Thu, May 30, 2013 at 10:14 AM, Mark Abraham
>> > > >wrote:
>> >
>> > > No charges = no problem. You can trivially test this yourself with
>> mdrun
>> > > -rerun ;-)  Manual 4.1.4 talks about what RF is doing.
>> > >
>> > > Mark
>> > >
>> > >
>> > > On Thu, May 30, 2013 at 6:38 AM, Dallas Warren
>> > > > >wrote:
>> > >
>> > > > In a system that has no charges, should we observe a difference
>> between
>> > > > simulations using PME/Cut-offs or Reaction Field?
>> > > >
>> > > > >From my understanding there should not be, since there are no
>> charges
>> > > > which treatment you use shouldn't' make a difference.
>> > > >
>> > > > However, it does and I am trying to work out why.
>> > > >
>> > > > Any suggestions on the reason?
>> > > >
>> > > > What is it that Reaction Field is doing, does it influence
>> anything
>> > other
>> > > > than long range charge interactions?
>> > > >
>> > > > Catch ya,
>> > > >
>> > > > Dr. Dallas Warren
>> > > > Drug Discovery Biology
>> > > > Monash Institute of Pharmaceutical Sciences, Monash University
>> > > > 381 Royal Parade, Parkville VIC 3052
>> > > > dallas.war...@monash.edu
>> > > > +61 3 9903 9304
>> > > > -
>> > > > When the only tool you own is a hammer, every problem begins to
>> > resemble
>> > > a
>> > > > nail.
>> > > > --
>> > > > gmx-users mailing listgmx-users@gromacs.org
>> > > > http://lists.gromacs.org/mailman/listinfo/gmx-users
>> > > > * Please search the archive at
>> > > > http://www.gromacs.org/Support/Mailing_Lists/Search before
>> posting!
>> > > > * Please don't post (un)subscribe requests to the list. Use the
>> > > > www interface or send it to gmx-users-requ...@gromacs.org.
>> > > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> > > >
>> > > --
>> > > gmx-users mailing listgmx-users@gromacs.org
>> > > http://lists.gromacs.org/mailman/listinfo/gmx-users
>> > > * Please search the archive at
>> > > http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> > > * Please don't post (un)subscribe requests to the list. Use the
>> > > www interface or send it to gmx-users-requ...@gromacs.o

Re: [gmx-users] Re: Nose-Hover chains for membrane protein simulation

[Michael Shirts]

And I my point was I didn't think that there was going to be a
measurable ergodicity difference between NH chains and NH.  Given the
size of the system, the chaoticness of atomic motions will likely give
you configurational sampling indistinguishable from full ergodicity.
Most of the errors that are solved by NH chains are for small toy
systems.

On Sun, Jun 2, 2013 at 2:17 AM, James Starlight  wrote:
> Michael,
>
>
> thanks for suggestions.
>
> the main reason of ussage N-H with chains is the assumption that simple N-H
> does not provide ergodicity of the system assuming that I'd like to sample
> all temperature acceptable conformations on the 100 ns trajectory.
>
> But as I understood the chain regime does not compatible with the membrane
> protein simulation due to the artifacts arising with MTTK batostat.
>
> James
>
> 2013/6/1 Michael Shirts 
>
>> I can't think of any instance where nose-hoover chains provides an
>> advantage over nose-hoover in a large system -- all the demonstrations
>> of superiority are in model systems that are not particularly chaotic.
>>  As the system gets more chaotic, it matters less.
>>
>> I would go with md, nose-hoover (w/o chains), and Parrinello-Rahman
>> with semiisotropic scaling.
>>
>> On Sat, Jun 1, 2013 at 1:07 AM, James Starlight 
>> wrote:
>> > Performing 5ns simulation after 2 ns equilibration in NPT ensemble (MTTK
>> > barostat 5ps coupling ) I've observed non-physical behaviour of my system
>> > with the constant drift of the protein molecule as the rigid body in the
>> > y-z plane
>> >
>> > Energy  Average   Err.Est.   RMSD  Tot-Drift
>> >
>> ---
>> > Pressure   -760.137 --193.776218.059
>>  (bar)
>> >
>> >
>> > >From manual I've noticed that MMTK doest not support *semiisotropic
>> > scalling.  *Doest it means that with the Nose-hover chains and md-vv I
>> > should use only weaker coupling during productions runs (I cant use
>> > Parinello;s barostat with such options too)
>> >
>> > Thanks for help
>> >
>> > James
>> >
>> >
>> >
>> > 2013/5/31 James Starlight 
>> >
>> >> Dear Gromacs users!
>> >>
>> >> I'd like to perform simulation of the membrane protein in lipid-water
>> >> system using Nose-Hover with chains.
>> >>
>> >> From manual I've found that with such thermostat I should use (1) md-vv
>> >> integrator (2) MTTK  instead of Parinello's batostat  and (3) shake
>> instead
>> >> of LINCS.
>> >>
>> >>
>> >> How doest such options compatible with the simulation of membrane
>> proteins
>> >> in general ? On what other options should I pay attention during
>> simulation
>> >> of membrane protein with NH chains ?
>> >>
>> >>
>> >>
>> >> Thanks for help,
>> >> James
>> >>
>> > --
>> > gmx-users mailing listgmx-users@gromacs.org
>> > http://lists.gromacs.org/mailman/listinfo/gmx-users
>> > * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> > * Please don't post (un)subscribe requests to the list. Use the
>> > www interface or send it to gmx-users-requ...@gromacs.org.
>> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Nose-Hover chains for membrane protein simulation

[Michael Shirts]

I can't think of any instance where nose-hoover chains provides an
advantage over nose-hoover in a large system -- all the demonstrations
of superiority are in model systems that are not particularly chaotic.
 As the system gets more chaotic, it matters less.

I would go with md, nose-hoover (w/o chains), and Parrinello-Rahman
with semiisotropic scaling.

On Sat, Jun 1, 2013 at 1:07 AM, James Starlight  wrote:
> Performing 5ns simulation after 2 ns equilibration in NPT ensemble (MTTK
> barostat 5ps coupling ) I've observed non-physical behaviour of my system
> with the constant drift of the protein molecule as the rigid body in the
> y-z plane
>
> Energy  Average   Err.Est.   RMSD  Tot-Drift
> ---
> Pressure   -760.137 --193.776218.059  (bar)
>
>
> >From manual I've noticed that MMTK doest not support *semiisotropic
> scalling.  *Doest it means that with the Nose-hover chains and md-vv I
> should use only weaker coupling during productions runs (I cant use
> Parinello;s barostat with such options too)
>
> Thanks for help
>
> James
>
>
>
> 2013/5/31 James Starlight 
>
>> Dear Gromacs users!
>>
>> I'd like to perform simulation of the membrane protein in lipid-water
>> system using Nose-Hover with chains.
>>
>> From manual I've found that with such thermostat I should use (1) md-vv
>> integrator (2) MTTK  instead of Parinello's batostat  and (3) shake instead
>> of LINCS.
>>
>>
>> How doest such options compatible with the simulation of membrane proteins
>> in general ? On what other options should I pay attention during simulation
>> of membrane protein with NH chains ?
>>
>>
>>
>> Thanks for help,
>> James
>>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] issue in replica exchange

[Michael Shirts]

Summarizing!

On Fri, May 3, 2013 at 12:31 AM, XAvier Periole  wrote:
>
> Are confirming that you reproduce the problem with gmx-4.6.1 or simply 
> summarizing in case we lose track :))
>
> On May 2, 2013, at 23:31, Michael Shirts  wrote:
>
>> So to summarize -- the problem appears to be with particle decomposition.
>>
>> On Thu, May 2, 2013 at 4:15 PM, XAvier Periole  wrote:
>>>
>>> I'll look at the 4.6.1 version next week, I could install it but I got a 
>>> conflict between the environmental variable defining openMP variable but I 
>>> turned it off during compilation …
>>>
>>> You could try to run on particle decomposition to see if you get a problem 
>>> … it should one quite quick.
>>>
>>> On May 2, 2013, at 2:36 PM, Michael Shirts  wrote:
>>>
>>>> Both.  So if 4.6.1 doesn't work, I want to know so we can patch it
>>>> before 4.6.2 comes out.  If it does work, then there is probably stuff
>>>> that can be backported.
>>>>
>>>> On Thu, May 2, 2013 at 8:32 AM, XAvier Periole  wrote:
>>>>>
>>>>> You mean working with or working on the code?
>>>>>
>>>>> I'll try gmx-4.6.1
>>>>>
>>>>> On May 2, 2013, at 2:26 PM, Michael Shirts  wrote:
>>>>>
>>>>>> Quick check here -- is 4.6 behaving correctly?  I actually spent some
>>>>>> time working on REMD in 4.6, and it seems to be behaving  correctly in
>>>>>> my hands with temperature and pressure control.
>>>>>>
>>>>>> Thanks for any additional info on this!
>>>>>>
>>>>>> On Thu, May 2, 2013 at 8:18 AM, Mark Abraham  
>>>>>> wrote:
>>>>>>> On Thu, May 2, 2013 at 12:58 PM, XAvier Periole  
>>>>>>> wrote:
>>>>>>>
>>>>>>>>
>>>>>>>> I saw that redmine report, which could be related but it seems to 
>>>>>>>> happen
>>>>>>>> only for runs done outside the domain and particle decompositions.
>>>>>>>>
>>>>>>>> I'll fill up a red mine.
>>>>>>>>
>>>>>>>> Anything I could do to help speeding the fix?
>>>>>>>
>>>>>>> What'd be really nice is some thought on how one can demonstrate that 
>>>>>>> the
>>>>>>> implementation of the exchange matches what would be expected from the
>>>>>>> theory. For T-exchange under NVT, it is sufficient to rescale velocities
>>>>>>> and quantities derived from them by the correct factor. That includes
>>>>>>> various things like T-coupling history and integrator half-step 
>>>>>>> quantities
>>>>>>> (and does REMD with leap-frog make sense anyway?). For NPT, there's
>>>>>>> probably also some P-coupling quantities to scale, and the box to 
>>>>>>> exchange.
>>>>>>> Anything I've missed? Hopefully virial contributions don't matter either
>>>>>>> way?
>>>>>>>
>>>>>>> Perhaps a decent first step is to hack the code to do a "self 
>>>>>>> exchange," by
>>>>>>> clearing the entire state and rebuilding with what would/should be 
>>>>>>> received
>>>>>>> from an exchange with a hypothethetical replica in an identical
>>>>>>> pre-exchange state. Only if the code can do that (i.e. mdrun -reprod
>>>>>>> produces a trajectory indistinguishable from a run that does not attempt
>>>>>>> this self exchange) is it worth considering proper state exchanges, and 
>>>>>>> the
>>>>>>> process of making the code do the former should illustrate what is 
>>>>>>> required
>>>>>>> for the latter.
>>>>>>>
>>>>>>> Mark
>>>>>>> --
>>>>>>> gmx-users mailing listgmx-users@gromacs.org
>>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>>>> * Please search the archive at 
>>>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>>>>> * Please don't post (un)subscribe requests to the list. Use the

Re: [gmx-users] issue in replica exchange

[Michael Shirts]

So to summarize -- the problem appears to be with particle decomposition.

On Thu, May 2, 2013 at 4:15 PM, XAvier Periole  wrote:
>
> I'll look at the 4.6.1 version next week, I could install it but I got a 
> conflict between the environmental variable defining openMP variable but I 
> turned it off during compilation …
>
> You could try to run on particle decomposition to see if you get a problem … 
> it should one quite quick.
>
> On May 2, 2013, at 2:36 PM, Michael Shirts  wrote:
>
>> Both.  So if 4.6.1 doesn't work, I want to know so we can patch it
>> before 4.6.2 comes out.  If it does work, then there is probably stuff
>> that can be backported.
>>
>> On Thu, May 2, 2013 at 8:32 AM, XAvier Periole  wrote:
>>>
>>> You mean working with or working on the code?
>>>
>>> I'll try gmx-4.6.1
>>>
>>> On May 2, 2013, at 2:26 PM, Michael Shirts  wrote:
>>>
>>>> Quick check here -- is 4.6 behaving correctly?  I actually spent some
>>>> time working on REMD in 4.6, and it seems to be behaving  correctly in
>>>> my hands with temperature and pressure control.
>>>>
>>>> Thanks for any additional info on this!
>>>>
>>>> On Thu, May 2, 2013 at 8:18 AM, Mark Abraham  
>>>> wrote:
>>>>> On Thu, May 2, 2013 at 12:58 PM, XAvier Periole  wrote:
>>>>>
>>>>>>
>>>>>> I saw that redmine report, which could be related but it seems to happen
>>>>>> only for runs done outside the domain and particle decompositions.
>>>>>>
>>>>>> I'll fill up a red mine.
>>>>>>
>>>>>> Anything I could do to help speeding the fix?
>>>>>>
>>>>>
>>>>> What'd be really nice is some thought on how one can demonstrate that the
>>>>> implementation of the exchange matches what would be expected from the
>>>>> theory. For T-exchange under NVT, it is sufficient to rescale velocities
>>>>> and quantities derived from them by the correct factor. That includes
>>>>> various things like T-coupling history and integrator half-step quantities
>>>>> (and does REMD with leap-frog make sense anyway?). For NPT, there's
>>>>> probably also some P-coupling quantities to scale, and the box to 
>>>>> exchange.
>>>>> Anything I've missed? Hopefully virial contributions don't matter either
>>>>> way?
>>>>>
>>>>> Perhaps a decent first step is to hack the code to do a "self exchange," 
>>>>> by
>>>>> clearing the entire state and rebuilding with what would/should be 
>>>>> received
>>>>> from an exchange with a hypothethetical replica in an identical
>>>>> pre-exchange state. Only if the code can do that (i.e. mdrun -reprod
>>>>> produces a trajectory indistinguishable from a run that does not attempt
>>>>> this self exchange) is it worth considering proper state exchanges, and 
>>>>> the
>>>>> process of making the code do the former should illustrate what is 
>>>>> required
>>>>> for the latter.
>>>>>
>>>>> Mark
>>>>> --
>>>>> gmx-users mailing listgmx-users@gromacs.org
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>> * Please search the archive at 
>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>>> * Please don't post (un)subscribe requests to the list. Use the
>>>>> www interface or send it to gmx-users-requ...@gromacs.org.
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>> --
>>>> gmx-users mailing listgmx-users@gromacs.org
>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>> * Please search the archive at 
>>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>> * Please don't post (un)subscribe requests to the list. Use the
>>>> www interface or send it to gmx-users-requ...@gromacs.org.
>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> --
>>> gmx-users mailing listgmx-users@gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> * Please search the archive at 
>>> http://www.gromacs.org/Support/Mailing_Li

Re: [gmx-users] issue in replica exchange

[Michael Shirts]

Both.  So if 4.6.1 doesn't work, I want to know so we can patch it
before 4.6.2 comes out.  If it does work, then there is probably stuff
that can be backported.

On Thu, May 2, 2013 at 8:32 AM, XAvier Periole  wrote:
>
> You mean working with or working on the code?
>
> I'll try gmx-4.6.1
>
> On May 2, 2013, at 2:26 PM, Michael Shirts  wrote:
>
>> Quick check here -- is 4.6 behaving correctly?  I actually spent some
>> time working on REMD in 4.6, and it seems to be behaving  correctly in
>> my hands with temperature and pressure control.
>>
>> Thanks for any additional info on this!
>>
>> On Thu, May 2, 2013 at 8:18 AM, Mark Abraham  
>> wrote:
>>> On Thu, May 2, 2013 at 12:58 PM, XAvier Periole  wrote:
>>>
>>>>
>>>> I saw that redmine report, which could be related but it seems to happen
>>>> only for runs done outside the domain and particle decompositions.
>>>>
>>>> I'll fill up a red mine.
>>>>
>>>> Anything I could do to help speeding the fix?
>>>>
>>>
>>> What'd be really nice is some thought on how one can demonstrate that the
>>> implementation of the exchange matches what would be expected from the
>>> theory. For T-exchange under NVT, it is sufficient to rescale velocities
>>> and quantities derived from them by the correct factor. That includes
>>> various things like T-coupling history and integrator half-step quantities
>>> (and does REMD with leap-frog make sense anyway?). For NPT, there's
>>> probably also some P-coupling quantities to scale, and the box to exchange.
>>> Anything I've missed? Hopefully virial contributions don't matter either
>>> way?
>>>
>>> Perhaps a decent first step is to hack the code to do a "self exchange," by
>>> clearing the entire state and rebuilding with what would/should be received
>>> from an exchange with a hypothethetical replica in an identical
>>> pre-exchange state. Only if the code can do that (i.e. mdrun -reprod
>>> produces a trajectory indistinguishable from a run that does not attempt
>>> this self exchange) is it worth considering proper state exchanges, and the
>>> process of making the code do the former should illustrate what is required
>>> for the latter.
>>>
>>> Mark
>>> --
>>> gmx-users mailing listgmx-users@gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> * Please search the archive at 
>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> * Please don't post (un)subscribe requests to the list. Use the
>>> www interface or send it to gmx-users-requ...@gromacs.org.
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] issue in replica exchange

[Michael Shirts]

Quick check here -- is 4.6 behaving correctly?  I actually spent some
time working on REMD in 4.6, and it seems to be behaving  correctly in
my hands with temperature and pressure control.

Thanks for any additional info on this!

On Thu, May 2, 2013 at 8:18 AM, Mark Abraham  wrote:
> On Thu, May 2, 2013 at 12:58 PM, XAvier Periole  wrote:
>
>>
>> I saw that redmine report, which could be related but it seems to happen
>> only for runs done outside the domain and particle decompositions.
>>
>> I'll fill up a red mine.
>>
>> Anything I could do to help speeding the fix?
>>
>
> What'd be really nice is some thought on how one can demonstrate that the
> implementation of the exchange matches what would be expected from the
> theory. For T-exchange under NVT, it is sufficient to rescale velocities
> and quantities derived from them by the correct factor. That includes
> various things like T-coupling history and integrator half-step quantities
> (and does REMD with leap-frog make sense anyway?). For NPT, there's
> probably also some P-coupling quantities to scale, and the box to exchange.
> Anything I've missed? Hopefully virial contributions don't matter either
> way?
>
> Perhaps a decent first step is to hack the code to do a "self exchange," by
> clearing the entire state and rebuilding with what would/should be received
> from an exchange with a hypothethetical replica in an identical
> pre-exchange state. Only if the code can do that (i.e. mdrun -reprod
> produces a trajectory indistinguishable from a run that does not attempt
> this self exchange) is it worth considering proper state exchanges, and the
> process of making the code do the former should illustrate what is required
> for the latter.
>
> Mark
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Free Energy Calculations in Gromacs

[Michael Shirts]

You have to change atom types.  For example:

[ atomtypes ]
;name  bond_typemasscharge   ptype  sigma  epsilon
h1h1 0.  0.  A  2.47135e-01  6.56888e-02
h1_pert h1 0.  0.  A  2.47135e-01  3.56888e-02
; perturbed

The mass and charge can be zero, because they will be defined in the [
atoms ] part

[ atoms ]
;   nrtype  resnr residue  atom  cgnr  chargemasstypeB
 chargeB  MassB
1 h1  1 BLAHH1a 1 -0.0891.008h1_pert
 -0.030  1.008;  perturbed



On Sat, Apr 20, 2013 at 4:04 PM, HANNIBAL LECTER  wrote:

> Hi,
>
> I am trying to perform expanded ensemble simulations between 2 states A
> (lambda=0) and B (lambda =1) with 6 intermediate lambda values.
>
> In state B the Hamiltonian is rescaled, such that the epsilons of the vdW
> interactions, the charges, the bond, angle and dihedral constants are a
> multiple of the similar terms of State A.
>
> It's not quite clear to me (going through the archives of the gmx-users
> mailing list how to perform these calculations. One way to do this, is to
> use a single topology file in which the charges and the other terms are
> specified for both states A and B. However, it is not clear as to how
> should I scale the epsilons in the topology file. (My atoms are not mutated
> in state B. They are the same atoms with scaled epsilons.)
>
> Secondly, since I have the topology files for states A and B, is there a
> way I could perform the simulations (any particular way in grompp) where
> both the topology files can be preprocessed designating the end states and
> using the mdp options, the simulations corresponding to the intermediate
> lambda values performed??
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Hamiltonian replica exchange umbrella sampling with gmx 4.6

[Michael Shirts]

Hi, Dejun-

Right now, the vector of lambda parameters is simply vdw, coul, bonded,
restraint, temperature.  You can't have, say, 2 different coul vectors or
two different restraint vectors for different restraints.  But you can
change any of these components.

You define the vector manually by writing out the list.  So to do 2D in vdw
and restraints, you would define states:

vdw-lambdas=  0.0 0.0 0.0 0.5 0.5 0.5 1.0 1.0 1.0
restraint-lambdas  =  0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0

Multidimensional movement in parameter space with replica exchange can be
done with the md setting '-nex Nswaps'.  This does multiple pair swaps,
instead of just neighbor pair swaps.  I'd suggest something like Nswaps =
N^3, where N is the number of lambda states.  This is a close approximation
to Gibbs sampling (http://jcp.aip.org/resource/1/jcpsa6/v135/i19/p194110_s1),
where _all_ permutations are selected based on their Boltzmann weight.
 It's an approximation because with a finite number of swaps, you don't
quite get uncorrelated moves in state space, but it is rigorously correct
thermodynamically.  It is more powerful than randomly selecting a dimension
to move in, since it allows moves in all dimensions simultaneously with
proper weight.

Note that the number of lambdas states is hardcoded as 1024, but that can
(I think) be edited as desired without causing direct problems (other than
perhaps needing to make the .mdp lines longer).

We are working on more general multistate lambda vector formalism for 5.0.
 If you have suggestions, let me know.

I'd be happy to look over input files or give additional advice on a
specific setup.


On Tue, Apr 2, 2013 at 12:52 PM, Dejun Lin  wrote:

> Hi Michael,
>
> Do the codes now support walking in multidimensional parameter space? i.e.,
> a state is defined by a set of lambda parameters {l1,l2,l3,...,ln} and a MC
> move is attempted along one of the parameter, which is randomly picked.
>
> Thanks,
> Dejun
>
>
>
> --
> View this message in context:
> http://gromacs.5086.n6.nabble.com/Hamiltonian-replica-exchange-umbrella-sampling-with-gmx-4-6-tp5006187p5006842.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Hamiltonian replica exchange umbrella sampling with gmx 4.6

[Michael Shirts]

Great -- if it doesn't seem to be working the way it should after some
playing around, then submit it as a redmine issue, and I'll take a
look.

On Fri, Mar 8, 2013 at 2:51 AM, Joakim Jämbeck  wrote:
> Dear Michael,
>
> Thank you for your reply.
>
> Yes, it is relatively clear now. Will try to play around with this later 
> today.
>
> Best,
> Joakim
>
>
>> Date: Thu, 7 Mar 2013 08:55:31 -0500
>> From: Michael Shirts 
>> Subject: Re: [gmx-users] Hamiltonian replica exchange umbrella
>>   sampling with   gmx 4.6
>> To: Discussion list for GROMACS users 
>>
>> Hi, Joakim-
>>
>> Hamiltonian exchange only should work if there is a lambda coupling
>> parameter that defines the potential at each state.  You need to
>> define your pulling potential so that the coupling-lambda parameter
>> can be used to define the different pulling location centers along
>> your trajectory.  Does this make it clearer?
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Hamiltonian replica exchange umbrella sampling with gmx 4.6

[Michael Shirts]

Hi, Joakim-

Hamiltonian exchange only should work if there is a lambda coupling
parameter that defines the potential at each state.  You need to
define your pulling potential so that the coupling-lambda parameter
can be used to define the different pulling location centers along
your trajectory.  Does this make it clearer?


On Thu, Mar 7, 2013 at 7:26 AM, Joakim Jämbeck  wrote:
> Dear gmx-users,
>
> I am currently trying to runt Hamiltonian replica exchange umbrella sampling 
> in hope to do some better sampling.
>
> I have generated a number of tpr-files along my reaction coordinate and they 
> all run fine in independent simulations. The issue comes when I would like 
> the replica exchange to start.
>
> The following line is used to initiate the exchange:
>
> mdrun_mpi -replex 1000 -s md -pf pullf -px pullx -multi 40 -maxh 0.5
>
> All replicas have the same temperature and the following error is what I face 
> seconds after submitting the job:
>
> The properties of the 40 systems are all the same, there is nothing to 
> exchange
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
>
> I could simply change the temperature between the replicas by 0.001 K and it 
> would run I think. But that is not very elegant.
>
> Does anyone have any suggestions?
>
> Thanks in advance!
>
> Best regards,
> Joakim--
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] CMAP and Free Energy

[Michael Shirts]

There is no theoretical reason to exclude it. The CMAP code is routed
differently in the logic, and was put in at the same time the free
energy code was, so it's just software engineering issues.

This is a good candidate for inclusion in 5.0 if enough people request.

On Wed, Feb 27, 2013 at 1:04 PM, francesco oteri
 wrote:
> Thank you so much!
> I am gonna file a feature request
>
>
> 2013/2/27 Mark Abraham 
>
>> Unfortunately, there's lots of areas of complex software projects where the
>> interaction of some new feature with some old feature is not of enough
>> interest to the people who maintained/developed those features for them to
>> put work into them. You can file a feature request at
>> redmine.gromacs.orgso that people know there is some interest in it -
>> but don't hold your
>> breath waiting!
>>
>> Offhand, I can't think of a theoretical reason why that'd be a problem, but
>> I'm an expert in neither FE calcs nor the details of CMAP :-)
>>
>> Mark
>>
>> On Wed, Feb 27, 2013 at 12:08 PM, francesco oteri <
>> francesco.ot...@gmail.com
>> > wrote:
>>
>> > Dear gromacs users,
>> > can someone of you explain me why CMAP is not managed in  Free Energy?
>> > I guess that the force derived by CMAP can be traited as any other force.
>> > Is there any theoretical problem in doing:
>> > CMAP(lambda) = (lambda)*CMAP +(1-lambda)*CMAP
>> > ?
>> >
>> > Francesco
>> > --
>> > gmx-users mailing listgmx-users@gromacs.org
>> > http://lists.gromacs.org/mailman/listinfo/gmx-users
>> > * Please search the archive at
>> > http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> > * Please don't post (un)subscribe requests to the list. Use the
>> > www interface or send it to gmx-users-requ...@gromacs.org.
>> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>
>
>
> --
> Cordiali saluti, Dr.Oteri Francesco
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] g_bar for larger systems (protein-protein interaction)

[Michael Shirts]

My personal opinion is that for large protein-protein calculations,
the free energy should be computed through potential of mean force
calculations, NOT alchemical methods, using the endpoints (properly
corrected) to determine the free energy of association.

There are a number of tutorials and example papers on how to compute
these potentials of mean force, so I won't go into those details here.

On Mon, Feb 25, 2013 at 3:49 PM, Ricardo O. S. Soares
 wrote:
> Hello dear users,
>
> I just took a look at the new gromacs paper at Bioinformatics, and I noticed
> a new tool g_bar, which is used for free energy calculations. I also found a
> nice tutorial by Justin Lemkul. While I'm adapting the workflow to my
> system, I'd like to know from anyone with more expertise in this field, if
> this method can be successfully applied to larger molecules as ligands,
> other than small organic ligands. What is the effectiveness in determining a
> protein-protein (300+ residues) energy of binding, if any?
>
> Thanks,
>
> Ricardo.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] The time for the temperature and pressure coupling

[Michael Shirts]

Ah, now perhaps I see that I misread the question - it could have been
phrased more clearly.  If Erik understood it correctly, then the
answer to the question is: It depends on the integrator.  The
simulation is not constrained to a particular temperature or pressure
- rather, the dynamics are modified such that the ensemble is
consistent with the external temperature and pressure are specified.
Exactly where in the integration routine the dynamics are modified
depends on the method used.  See the manual for more of the details!

On Thu, Feb 7, 2013 at 3:59 AM, Erik Marklund  wrote:
> Hi,
>
> Perhaps a side point: Temperature and pressure can not be seen as
> constraints to the system at any given instant in the sense that e.g. the
> instantaneous kinetic energy perfectly match the temperature at every time
> step just because you have a thermostat. Time and ensemble averages will,
> however, reflect the temperature and pressure coupling.
>
> Erik
>
>
> On Feb 6, 2013, at 11:28 PM, Bao Kai wrote:
>
>> Dear Gromacs Team,
>>
>> I have a small question related to the scheme of the MD in Gromacs.
>>
>> When are the temperature and pressure constrains are enforced, before
>> the update of the velocity and position or after?
>>
>> Thank you very much.
>>
>> Best Regards,
>> Kai
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use thewww interface
> or send it to gmx-users-requ...@gromacs.org.
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] The time for the temperature and pressure coupling

[Michael Shirts]

The coordinates and velocities that are printed (and that are used to
calculate the properties like energy, virial, etc) are always
consistent with the constraints.  The exact order of how things are
done often depends on the integrator.  For example, velocity scaling
can be done before or after constraints, because there are no
velocities parallel to the bond vector, so after velocity scaling, the
constraints are still obeyed.

On Wed, Feb 6, 2013 at 5:28 PM, Bao Kai  wrote:
> Dear Gromacs Team,
>
> I have a small question related to the scheme of the MD in Gromacs.
>
> When are the temperature and pressure constrains are enforced, before
> the update of the velocity and position or after?
>
> Thank you very much.
>
> Best Regards,
> Kai
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Expanded Ensemble and Gromacs 4.6

[Michael Shirts]

Hi, Joakim-

Expanded ensemble is still a bit experimental.  I don't immediately
see any problem that jump right out, but if you go to
http://redmine.gromacs.org/ and file a bug report, including giving
example files that cause the problem, I can take a look at it.

On Tue, Feb 5, 2013 at 6:00 AM, Joakim Jämbeck  wrote:
> Hi,
>
> I am trying to use the Expanded Ensemble (EE) method to compute the free 
> energy of solvation of a small organic molecule.
>
> Basically I am playing around but I cannot get the simulations to run.
>
> Here are my EE and free energy settings:
>
> %---
> free-energy = expanded  ; Expanded ensemble
> couple-moltype  = C1X   ; Molecule to introduce
> couple-lambda0  = none  ; Go from no interactions with 
> solvent...
> couple-lambda1  = vdw-q ; ...to full interactions
> couple-intramol = no; Do not decouple internal 
> interactions
> init_lambda_state   = 0 ; Start from the first column of the 
> lambda vectors
> delta-lambda= 0 ; No increments in lambda
> nstdhdl = 100   ; Frequency for writing dH/dlambda
> coul-lambdas= 0.0 0.25 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
> vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0
> dhdl-print-energy   = yes   ; Include total energy in the dhdl 
> file
> sc-alpha= 0.5   ; Soft core alpha parameter
> sc-power= 1 ; Power of lambda in soft core 
> function
> sc-r-power  = 6 ; Power of radial term in the soft 
> core function
> sc-sigma= 0.3   ; Soft core sigma
> sc-coul = no; No soft core for Coulomb
> separate-dhdl-file  = yes   ; Seperate dhdl files
> dhdl-derivatives= yes   ; Print derivatives of the Hamiltonian
>
> ; expanded ensemble variables
> nstexpanded = 100   ; Number of steps between attempts to 
> change the Hamiltonian
> lmc-stats   = wang-landau   ; WL algorithm to explore state space
> lmc-move= gibbs ; Decides which state to move to
> lmc-weights-equil   = number-steps  ; EE weight updating stops after...
> weight-equil-number-steps = 500 ; ...10ns of simulation
>
> ; Seed for Monte Carlo in lambda space
> lmc-seed= 1993
> lmc-repeats = 1
> lmc-gibbsdelta  = -1
> lmc-forced-nstart   = 0
> symmetrized-transition-matrix = no
> nst-transition-matrix   = -1
> mininum-var-min = 100
> wl-scale= 0.6
> wl-ratio= 0.8
> init-wl-delta   = 1
> wl-oneovert = yes
> %---
>
> Besides this I use LINCS to constrain all bonds, sd-integrator and a "normal" 
> cut-off scheme with a 1 fs time step.
>
> However, once I try to run the files on the cluster I always end up with 
> LINCS warnings and after a few seconds the program crashes due to too many 
> LINCS warnings. If I increase the time step to 2 fs I run into problems 
> SETTLE for the water. I always start from a nice structure that is taken as 
> the final snap shot from a simple 1 ns MD simulation of my system.
>
> What could be the problem?
>
> If I use free_energy = yes instead of EE things work fine.
>
> Did I perhaps mess up the EE settings or something?
>
> Thanks in advance.
>
> Best,
> Joakim
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] meta-dynamics in gromacs-4.6

[Michael Shirts]

I assume PLUMED will be implemented for Gromacs 4.6, as many PLUMED
developers use Gromacs.  Perhaps any PLUMED lurkers on the list can
speak up. . . .

On Wed, Jan 16, 2013 at 9:20 AM, Mark Abraham  wrote:
> The GROMACS team has no plans for that. The usual problem here is that
> everybody would like every algorithm included, but that developers with
> time and experience are scarce :-) It's an open source project though, so
> anyone can do whatever they like. We're prepared to consider inclusions to
> the main project.
>
> Note that we're doing a major rework of the code base to C++ in the next
> year (or more), so people implementing new features may wish to consider
> that in their choice to write code :-)
>
> Mark
>
> On Tue, Jan 15, 2013 at 2:26 PM, James Starlight 
> wrote:
>
>> Dear Gromacs developers!
>>
>>
>> There is well-known plugin plumed which can be used for implementation
>> of meta-dynamics simulation if Gromacs-4.5. I wounder to know if it
>> possible to include some meta-dynamics options in the new gromacs
>> release ( similar to inclusion of essential dynamics sampling in
>> previous gromacs versions) ?
>>
>>
>> Thanks for attention,
>> James
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Re: Is vacuum simulation NVT?

[Michael Shirts]

> Could you tell me is there any difference of different Tau_t ussage (
> inverse friction in case of Stochastic dynamics) for simulation of
> water-soluble as well as membrane-proteins ? In the first case I'm
> using tau_t 2ps that is lower than internal water friction. In the
> second case one part of protein could be in the membrane ( e.g
> helixes) but other ( e.g loops) in water media. Both of that solvents
> have different characteristic viscousity. So what Tau_t should be used
> in stochastic dynamics for such biphastic systems?

I don't think there is a "right" sd tau_t.  You probably want at
tau_t^-1 that is lower (tau_t longer) than the natural viscosity of
either system, such that the viscosity of the actual intermolecular
interactions dominates, not the unphysical viscosity of the
thermostat.  I don't know how much lower is good enough, though.  You
may have to experiment.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Is vacuum simulation NVT?

[Michael Shirts]

> In the absence of PBC, you simply have an infinite system.  In a loose
> sense, that may be NVT, but V is infinite, so whether or not you can
> consider that to be constant or not is theoretical math above what I know :)

A real molecule in vacuum is usually NVE -- it is not coupled to the
environment, and thus must have conserved energy.  You can certainly
add a thermostat, and then it will be NVT, though it won't be very
much like a real isolated gas molecule.  If you try to run NPT, the
simulation will likely crash because of numerical instabilities, and
there's not much of a point, since you are essentially either 1) in
the ideal gas limit if running with no periodic boundary conditions 2)
in some sort of weird superdilute crystal that really doesn't resemble
anything real if run with a periodic boundary conditions

If you are subtracting out the center of mass motion, then V is not
infinite -- you remove the center of mass degree of freedom, and thus
you have a very different ensemble than if you include the center of
mass motion.  You would need to multiply by V to get the partition
function for an actual gas.

Note that I would strongly suggest sd as the integrator/thermostat,
since there are real issues with ergodicity in systems with only a few
degrees of freedom
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Question about conserved energy in MTTK

[Michael Shirts]

Hi, all-

I would recommend using Parrinellio-Rahman + Nose-Hoover md + at least
until 4.6.

A random-walk drift in the conserved energy is actually what MTTK
gives -- it's not as conserved as, say, energy conservation, it just
has an expectation value of zero drift over time, which means that the
RMSD will increase with time according to sqrt(dt).

But if you are seeing constant drift, something is wrong.

On Wed, Nov 28, 2012 at 7:20 AM, tarak karmakar  wrote:
> Hi,
>
> I was also facing the same problem. If you check your pressure during
> this NPT run, u can see that it got increased to a higher value. I had
> posted the same problem few days back, u can follow the thread. It
> seems MTTK is not stable enough and is not performing well in this
> context. So I have moved to Leap-frog, Nose-Hoover, Parinello-Rahman
> combination for the NPT simulation. There is one paper as well by
> Prof. Shirts in JCTC.
>
> Cheers,
> Tarak
>
> On Thu, Nov 22, 2012 at 2:08 PM, Shun Sakuraba  
> wrote:
>> Dear list,
>>
>> I am trying to use MTTK barostat in GROMACS 4.5.5.
>> After analyzing the result for a while, I found that the conserved energy 
>> (not total energy) of MTTK is drifting during the simulation.
>> The .xvg, .edr files are uploaded at [1] and [2]. It is drifting with a 
>> constant ratio of ca. -185 kJ/mol/ps.
>>
>> I cannot believe this is an expected behavior, so could anyone point out 
>> where I am wrong in my simulation setup? I found similar report at [3] but 
>> seems it was when 4.5 was in pre-release stage.
>>
>> Thanks in advance for your help!
>>
>> * Simulation detail
>> The system consists of 1000 SPC-E water molecules, and the time step is set 
>> to 0.5 fs, just in case the long timestep harms the conserevation (c.f. 
>> [3]). The interaction energy is set to switching version, just in case, too. 
>> Changing these parameter does not seem to improve the conservation.
>> The double precision version of GROMACS is used (single precision version 
>> also has the same problem).
>> The system has been pre-equilibrated with Berendsen pressure coupling 
>> simulation with the same pressure and temperature.
>>
>> [1] https://www.dropbox.com/s/16bgfhvavqcw1f8/mttk05.xvg
>> [2] https://www.dropbox.com/s/tg8butw39rmz8qk/mttk05.edr
>> [3] 
>> http://lists.gromacs.org/pipermail/gmx-developers/2010-January/003979.html
>>
>> == .mdp file contents follow
>>
>> integrator = md-vv
>> define =
>>
>> dt = 0.0005
>> nsteps = 100 ; 500 ps
>>
>> coulombtype = PME-Switch
>> vdwtype = Switch
>> pbc = xyz
>>
>> rlist = 1.2
>> rcoulomb = 1.0
>> rcoulomb_switch = 0.9
>> rvdw = 1.0
>> rvdw_switch = 0.9
>> nstlist = 1
>>
>> tinit = 0
>> tcoupl = nose-hoover
>> tc_grps = System
>> tau_t = 0.5
>> ref_t = 300.0
>> nsttcouple = 1
>>
>> pcoupl = MTTK
>> pcoupltype = isotropic
>> compressibility = 4.5e-5
>> ref_p = 1.01325
>> tau_p = 0.5
>> refcoord_scaling = no
>> nstpcouple = 1
>>
>> constraints = hbonds
>> constraint_algorithm = LINCS
>>
>> nstxtcout = 100
>> nstlog = 100
>> nstenergy = 100
>> nstvout = 0
>> nstxout = 1000
>>
>> --
>> Shun SAKURABA, Ph.D.
>> Postdoc @ Molecular Modeling & Simulation Group, Japan Atomic Energy Agency
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at 
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: pressure_coupling

[Michael Shirts]

It's in review with JCTC right now.

On Thu, Nov 22, 2012 at 2:19 PM, ABEL Stephane 175950
 wrote:
> Hello,
>
> This is a very nice and interesting work, Michael. Thank you for the efforts 
> you made in writing this paper. I hope you will publish it.
>
> Best
>
> Stephane
>
>
> 
>
> Hi, all-
>
> There are some issues with MTTK + constraints that are being worked
> out for 4.6.  The good thing is, I have developed some sensitive tests
> of the correct volume distribution (see
> http://arxiv.org/abs/1208.0910) and the errors in PR are very, very
> small. I would recommend using md + PR for projects with code before
> 4.6.
>
> On Thu, Nov 22, 2012 at 4:27 AM, Florian Dommert
>  wrote:
>>> -Urspr?ngliche Nachricht-
>>> Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>>> boun...@gromacs.org] Im Auftrag von tarak karmakar
>>> Gesendet: Donnerstag, 22. November 2012 10:15
>>> An: Discussion list for GROMACS users
>>> Betreff: Re: [gmx-users] pressure_coupling
>>>
>>> U r right FLorian
>>> I have also tried playing around the tau_p but in vain.
>>> Even in absence of any constraints, it is giving almost same result.
>>> Em thinking to move again to Leap-Frog, NH , PR.  I see people generally
>> use this
>>> combination a lot.
>>> Thanks
>>>
>>> Tarak
>>
>> Yes, that is right. The reason might be, that it is stable, working with
>> Leap-Frog and implemented in Gromacs. However actually PR does not produce
>> an NpT but an isoenthalpic ensemble. It also does not conform to both
>> pressure virial theorems (see appendix of the Nose paper cited in the
>> Gromacs manual). For this reason it would be very very good, if MTTK would
>> work in Gromacs, because it fulfills all requirements for an NpT ensemble.
>> On the other hand side the deviations vanish in the thermodynamic limit, so
>> if your system is large enough, there should be no significant difference.
>>
>> /Flo
>>
>>>
>>> On Wed, Nov 21, 2012 at 8:08 PM, Florian Dommert >> stuttgart.de> wrote:
>>> >
>>> >
>>> > ---
>>> > Florian Dommert
>>> > Dipl. Phys.
>>> >
>>> > Institut f?r Computerphysik
>>> > Universit?t Stuttgart
>>> > Allmandring 3
>>> > D-70569 Stuttgart
>>> >
>>> > Tel.: 0711-68563613
>>> > Fax: 0711-68563658
>>> >
>>> >
>>> >> -Urspr?ngliche Nachricht-
>>> >> Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>>> >> boun...@gromacs.org] Im Auftrag von tarak karmakar
>>> >> Gesendet: Mittwoch, 21. November 2012 15:03
>>> >> An: Discussion list for GROMACS users
>>> >> Betreff: Re: [gmx-users] pressure_coupling
>>> >>
>>> >> Thanks for the information Flo.
>>> >> Before doing NPT I have already equilibrated my system by heating it
>>> >> from
>>> > 0K to
>>> >> 300K in 300 ps, then the pressure has reached to 1 bar. Now while
>>> >> doing
>>> > NPT I'm
>>> >> getting the excess pressure.
>>> >> Is there any problem with the coupling constant ? I am checking it by
>>> > taking
>>> >> different tau_p values. Let's see.
>>> >>
>>> >
>>> > I don't think that playing around with the coupling constant will help
>> you.
>>> > You can set it to extreme values, but you won't see any difference.
>>> > The coupling constant determines, how fast the system pressure should
>>> > relax to the reference pressure. I would see a better possibility to
>>> > play around by simulating for a longer time. Then observing the
>>> > variation of the pressure in time, the size of the fluctuation and the
>>> > excess pressure. Perhaps something will change, but I don't think so.
>>> > I play around with the coupling constants but observed no change.
>>> >
>>> > Maybe, but this is really speculation, there is a problem with the
>>> > combination of constraints and MTTK. Please check the archives of the
>>> > user and developer list to obtain more information.
>>> >
>>> > /Flo
>>> >
>>> >>
>>> >> On Wed, Nov 21, 2012 at 1:16 AM, Florian Dommert >> >> stuttgart.de> wrote:
>>> >> >> -Urspr?ngliche Nachricht-
>>> >> >> Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>>> >> >> boun...@gromacs.org] Im Auftrag von Justin Lemkul
>>> >> >> Gesendet: Dienstag, 20. November 2012 18:33
>>> >> >> An: Discussion list for GROMACS users
>>> >> >> Betreff: Re: [gmx-users] pressure_coupling
>>> >> >>
>>> >> >>
>>> >> >>
>>> >> >> On 11/20/12 12:29 PM, tarak karmakar wrote:
>>> >> >> > Thanks Justin for the quick reply.
>>> >> >> > Is there any problem with the algorithms ??
>>> >> >> >
>>> >> >> > I have used Velocity Verlet , Nose-Hoover and MTTK combination.
>>> >> >> > SHAKE has been used to constrains the H-covalent bonds.
>>> >> >> > tau_t = 1 ps
>>> >> >> > tau_P = 1 ps
>>> >> >> > I got the mean pressure at ~130 bar.
>>> >> >> >
>>> >> >> > Previously with the same initial coordinates I have used
>>> >> >> > Leap-Frog, NH, Parinello-Rehman with LINCS to constrain H-covalent
>>> bonds.
>>> >> >> > tau_t was 0.1 ps
>>> >> >> > and tau P was 2 ps.
>>> >> >> > The I have seen the pressure fluctuating around 1 

Re: [gmx-users] pressure_coupling

[Michael Shirts]

Hi, all-

There are some issues with MTTK + constraints that are being worked
out for 4.6.  The good thing is, I have developed some sensitive tests
of the correct volume distribution (see
http://arxiv.org/abs/1208.0910) and the errors in PR are very, very
small. I would recommend using md + PR for projects with code before
4.6.

On Thu, Nov 22, 2012 at 4:27 AM, Florian Dommert
 wrote:
>> -Ursprüngliche Nachricht-
>> Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> boun...@gromacs.org] Im Auftrag von tarak karmakar
>> Gesendet: Donnerstag, 22. November 2012 10:15
>> An: Discussion list for GROMACS users
>> Betreff: Re: [gmx-users] pressure_coupling
>>
>> U r right FLorian
>> I have also tried playing around the tau_p but in vain.
>> Even in absence of any constraints, it is giving almost same result.
>> Em thinking to move again to Leap-Frog, NH , PR.  I see people generally
> use this
>> combination a lot.
>> Thanks
>>
>> Tarak
>
> Yes, that is right. The reason might be, that it is stable, working with
> Leap-Frog and implemented in Gromacs. However actually PR does not produce
> an NpT but an isoenthalpic ensemble. It also does not conform to both
> pressure virial theorems (see appendix of the Nose paper cited in the
> Gromacs manual). For this reason it would be very very good, if MTTK would
> work in Gromacs, because it fulfills all requirements for an NpT ensemble.
> On the other hand side the deviations vanish in the thermodynamic limit, so
> if your system is large enough, there should be no significant difference.
>
> /Flo
>
>>
>> On Wed, Nov 21, 2012 at 8:08 PM, Florian Dommert > stuttgart.de> wrote:
>> >
>> >
>> > ---
>> > Florian Dommert
>> > Dipl. Phys.
>> >
>> > Institut für Computerphysik
>> > Universität Stuttgart
>> > Allmandring 3
>> > D-70569 Stuttgart
>> >
>> > Tel.: 0711-68563613
>> > Fax: 0711-68563658
>> >
>> >
>> >> -Ursprüngliche Nachricht-
>> >> Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> >> boun...@gromacs.org] Im Auftrag von tarak karmakar
>> >> Gesendet: Mittwoch, 21. November 2012 15:03
>> >> An: Discussion list for GROMACS users
>> >> Betreff: Re: [gmx-users] pressure_coupling
>> >>
>> >> Thanks for the information Flo.
>> >> Before doing NPT I have already equilibrated my system by heating it
>> >> from
>> > 0K to
>> >> 300K in 300 ps, then the pressure has reached to 1 bar. Now while
>> >> doing
>> > NPT I'm
>> >> getting the excess pressure.
>> >> Is there any problem with the coupling constant ? I am checking it by
>> > taking
>> >> different tau_p values. Let's see.
>> >>
>> >
>> > I don't think that playing around with the coupling constant will help
> you.
>> > You can set it to extreme values, but you won't see any difference.
>> > The coupling constant determines, how fast the system pressure should
>> > relax to the reference pressure. I would see a better possibility to
>> > play around by simulating for a longer time. Then observing the
>> > variation of the pressure in time, the size of the fluctuation and the
>> > excess pressure. Perhaps something will change, but I don't think so.
>> > I play around with the coupling constants but observed no change.
>> >
>> > Maybe, but this is really speculation, there is a problem with the
>> > combination of constraints and MTTK. Please check the archives of the
>> > user and developer list to obtain more information.
>> >
>> > /Flo
>> >
>> >>
>> >> On Wed, Nov 21, 2012 at 1:16 AM, Florian Dommert > >> stuttgart.de> wrote:
>> >> >> -Ursprüngliche Nachricht-
>> >> >> Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-
>> >> >> boun...@gromacs.org] Im Auftrag von Justin Lemkul
>> >> >> Gesendet: Dienstag, 20. November 2012 18:33
>> >> >> An: Discussion list for GROMACS users
>> >> >> Betreff: Re: [gmx-users] pressure_coupling
>> >> >>
>> >> >>
>> >> >>
>> >> >> On 11/20/12 12:29 PM, tarak karmakar wrote:
>> >> >> > Thanks Justin for the quick reply.
>> >> >> > Is there any problem with the algorithms ??
>> >> >> >
>> >> >> > I have used Velocity Verlet , Nose-Hoover and MTTK combination.
>> >> >> > SHAKE has been used to constrains the H-covalent bonds.
>> >> >> > tau_t = 1 ps
>> >> >> > tau_P = 1 ps
>> >> >> > I got the mean pressure at ~130 bar.
>> >> >> >
>> >> >> > Previously with the same initial coordinates I have used
>> >> >> > Leap-Frog, NH, Parinello-Rehman with LINCS to constrain H-covalent
>> bonds.
>> >> >> > tau_t was 0.1 ps
>> >> >> > and tau P was 2 ps.
>> >> >> > The I have seen the pressure fluctuating around 1 bar( as
>> >> >> > expected) So can you please inform me from where this problem is
>> >> >> > coming - algorithms and/ tau_t and tau_P parameters ?
>> >> >> >
>> >> >>
>> >> >> I have no personal experience with the md-vv/MTTK combination.
>> >> >> The way to test if there is a bug or something is to take an
>> >> >> equilibrated system (as
>> >> > suggested
>> >> >> before) and continue it with the desired parameters.  If they
>> >> >> deviate or 

Re: [gmx-users] Re: Fast exchanges for REMD

[Michael Shirts]

> However, the time value (4 in this example) is limited to 6 digits.

Sounds like this should be increased?  There's a pending change to
replica exchange, so this could be added to 4.6 without disrupting the
release timing.

On Wed, Sep 26, 2012 at 11:22 AM, Andreas Zink  wrote:
> Dear all,
>
> I could finally demux my REMD trajectories with high EAF. They look fine,
> but I'm not 100% sure about it.
>
> Unfortunately, there seems to be a "bug" in mdrun. As you migh know, the log
> files contain the exchange attempts like:
>>
>> Replica exchange at step 2000 time 4
>> Repl ex  0123 x  45
>> Repl pr.00   1.0
>
> However, the time value (4 in this example) is limited to 6 digits. It's not
> really a bug because I think GROMACS recommends EAFs of max. 1/ps. But if
> you exchange every 0.5 ps (EAF= 2/ps) you can run the simulation for max.
> 9 ps only. Otherwise your log will simply chop the time after the
> decimal point, e.g.
>
>> Replica exchange at step ... time 9.5
>> ...
>>
>> Replica exchange at step ... time 10
>> ...
>>
>> Replica exchange at step ... time 10
>> ...
>>
>> Replica exchange at step ... time 11
>> ...
>
> I have written a Perl script which fixes the time values, based on the given
> number of steps and stepsize. Additionally the demux.pl script was changed
> because it also chops after 6 digits. I simply changed:
>>
>> printf(NDX "%-20g",$t);
>
> in the "pr_order" and "pr_revorder" subroutines to:
>>
>> printf(NDX "%-20.2f",$t);
>
> and it works just fine.
>
> Like I said the trajectories look fine, but I'm not really sure if it's
> actually correct that way. I would be happy if anyone would like to discuss
> this :)
>
>
> Cheers,
>
> Andi
>
>
>
>
> Am 24.09.2012 08:49, schrieb Andreas Zink:
>
>> Dear all,
>>
>> I've done some REMD simulations using a quite high exchange attempt
>> frequency (10 attempts per ps) as proposed by Sindhikara et al. ("Exchange
>> Often and Properly in Replica Exchange Molecular Dynamics",J. Chem. Theory
>> Comput. 2010, 6, 2804–2808 ).
>> Unfortunately, I have now recognized that the demux perl script cannot
>> account for an EAF which is higher than the saving frequency in the
>> trajectory.
>>
>> Comments from demux.pl:
>> # If your exchange was every N ps and you saved every M ps you can make
>> for
>> # the missing frames by setting extra to (N/M - 1). If N/M is not integer,
>> # you're out of luck and you will not be able to demux your trajectories
>> at all.
>>
>> In my case exchanges every 0.1 ps and saved every 5 ps
>>
>> Changing the demux.pl script, so that it writes the "replica_index.xvg"
>> with a higher precision (time in ps) should be no problem. However, my
>> question is: will this work together with trjcat? Does trjcat search for the
>> timeframe given in the first column of "replica_index.xvg", or does it links
>> each line to one saved timeframe? If so, could I just delete the additional
>> lines in "replica_index.xvg"?
>>
>> I would be really happy if someone could help me with this!
>> Thanks!
>>
>> Andi
>>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Experiences with Gromacs scaling on US supercomputer centers?

[Michael Shirts]

Hi all,

I'd be interested to know about people's experiences with Gromacs on
US national computing centers.  Which machines have it set up to scale
the
best?  We're putting in an XSEDE request soon, and I'm trying to
figure out which resource to request.  Our system is
semi-coarse-grained, using
reaction field electrostatics, so we don't have to worry about PME.
Of course, couldn't hurt to know about PME scaling as well.   I'm
interesting
scalings with 100K - 300K atoms.

Of course, best performance will probably change with 4.6 because of
all the setup tweaks, but let's start with 4.5 scaling info!

Best,

Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shi...@virginia.edu
(434)-243-1821
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: v-rescale

[Michael Shirts]

I've done some extensive testing (paper on testing method in the
works) and vrescale gives a very accurate ensemble very well for NVT.
Parrinello-Rahman and MTTK are the only algorithms that are correct
for NPT.  Berendsen barostat is not.  Note that there is a bug with
vrescale + md-vv + that is fixed in 4.6 and (hopefully) for the next
patch of 4.5.

For equilibrating a system, Berendsen for both temperature and
pressure is the best bet.  They artificially minimize fluctuations,
which is great for equilibration, bad for data collection.

On Thu, Sep 20, 2012 at 6:05 PM, Mark Abraham  wrote:
> On 20/09/2012 9:35 AM, Peter C. Lai wrote:
>>
>> I am not sure where the idea of using berendsen barostat with the
>> v-rescale
>> thermostat for equilibration came from, however. Doesn't the typical
>> equilibration begin with v-rescale for temperature equilibration then
>> adding parinello-rahman barostat then switching to nose-hoover for
>> production
>> runs (as nose-hoover chains result in the correct canonical distribution)?
>
>
> N-H does have known issues, see
> http://link.aip.org/link/doi/10.1063/1.2989802 and
> http://link.aip.org/link/doi/10.1063/1.2408420. I am not aware of any
> shortcomings of the Bussi v-rescale thermostat in GROMACS.
>
> Mark
>
>
>>
>> On 2012-09-19 04:24:27PM -0700, Ladasky wrote:
>>>
>>> Dear Sara,
>>>
>>> I just had a problem with my simulations that I traced to the use of the
>>> V-rescale temperature algorithm.  Here is my recent post:
>>>
>>>
>>> http://gromacs.5086.n6.nabble.com/Re-Water-molecules-cannot-be-settled-why-td4999302.html;cid=1348087067061-71#a5001121
>>>
>>> V-rescale may be appropriate in certain simulations, but it is apparently
>>> NOT appropriate when used with Berendsen pressure coupling during the
>>> initial equilibration.  I don't know if that is related to your problem,
>>> but
>>> it's something that I just discovered the hard way.
>>>
>>>
>>>
>>> --
>>> View this message in context:
>>> http://gromacs.5086.n6.nabble.com/v-rescale-tp5001066p5001122.html
>>> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
>>> --
>>> gmx-users mailing listgmx-users@gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> * Please don't post (un)subscribe requests to the list. Use the
>>> www interface or send it to gmx-users-requ...@gromacs.org.
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] BAR / g_bar problems

[Michael Shirts]

Hi, David-

Perhaps we can work with you to compare the internal g_bar
implementation with our external BAR and MBAR implementation run from
the .dhdl.xvg data output in 4.6.  It would probably be easier to run
these calculations after the optimization updates are added in the
next few days(?), but we can plan now.  Note that withe errors this
big, 100-200 ps should be enough to see what's going on, so it can be
done rapidly.  Drop me a line off the list to figure out details?

Best,
Michael

On Fri, Aug 24, 2012 at 9:22 AM, David van der Spoel
 wrote:
> Hi,
>
> we have terrible problems to get reasonable results from BAR free energy
> calculations. We basically follow Justin's tutorial for solvation of
> methanol, ethanol, diethyl-ether and neopentane in water using OPLS but get
> very strange results. The free energy curves are here:
>
> http://folding.bmc.uu.se/images/koko.jpg
>
> Molecule  Energy Exper
> Methanol -8-21
> Ethanol  -9-21
> Diethylether -11-7
> Neopentane   -10   +11
>
> any clue?
> --
> David van der Spoel, Ph.D., Professor of Biology
> Dept. of Cell & Molec. Biol., Uppsala University.
> Box 596, 75124 Uppsala, Sweden. Phone:  +46184714205.
> sp...@xray.bmc.uu.sehttp://folding.bmc.uu.se
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Only plain text messages are allowed!
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] When to use Dispersion Correction for Lipid Bilayers

[Michael Shirts]

Short answer -- use a dispersion correction if the force field was
parameterized to include one.  Make sure you use the cutoff that the
lipid was parameterized for as well.

Long answer -- neither on nor off is correct for a lipid bilayer.  A
dispersion corrections corrects for the fact that you are neglecting
the r^-6 term at long range.  However, it assumes an isotropic
distribution outside the cutoff.  Lipid bilayers have long range
order, so a standard dispersion correction is inappropriate.  The
lipid properties will be cutoff dependent, which is not a very good
thing.  See  PJ in't Veld, AE Ismail, GS Grest, J. Chem. Phys. 127,
144711 (2007) for a solution, using Ewald summation for the
Lennard-Jones terms.  This is in the works for Gromacs (and has been
for a while), but might still be a while because it's a little lower
on the to do list.  Once methods like this are in, it will be possible
to parameterize lipids that behave correctly.

Best,
Michael

On Sun, Aug 19, 2012 at 3:16 PM, David Ackerman  wrote:
> Hello,
>
> I was wondering when it is appropriate to use Dispersion Correction
> for lipid bilayers, or which setting (no, EnerPres, or Ener) is best.
> I have seen it used in some discussions, and not used in others. As
> for myself, I have simulated a DPPC bilayer. With DispCorr=EnerPres, I
> get an area per lipid of ~.615 nm^2 (whereas other reported values are
> closer to .64 nm^2) and slower diffusion than is reported. However,
> when I don't have a DispCorr term, my area per lipid becomes ~.656
> nm^2 and my diffusion more closely matches other reported values.
>
> So, should DispCorr be used at all for bilayer simulations, and if so,
> which type is most appropriate?
>
> Thank you for your help,
> David
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Only plain text messages are allowed!
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] lincs with mttk

[Michael Shirts]

Good question.  Short answer, no -- LINCS doesn't play well with a
velocity verlet based pressure control algorithm.

Long answer: MTTK has ended up not being robust because you need to
solve a self consistent set of equations every timestep using the
pressure estimator, which is extremely noisy, so it crashes randomly
for no good reason.

The good news is that I've done some fairly extensive testing, and
leapfrog MD + Parrinello-Rahman + Nose-Hoover or vrescale is very
close to the correct distribution -- there's really not that much
difference for any practical purpose (I have a preprint with more
details).  The volume fluctuations are almost exactly right. So I
think that for anything that doesn't require getting free energies to
kJ/mol accuracy you are OK with that combination (just don't use
Berendsen anything).

Longer term (5.0), MTTK + shake will be removed -- its just too
numerically unstable, and the self-consistent iteration part means
it's very hard to parallelize.  MTTK will be left in, but only without
constraints.  Hopefully, for 5.0 we will also have RESPA-like
functionality, so that the bondeds can be performed at a higher
frequency than the nonbondeds, which will be another useful way to get
longer times steps. We also have planned to implement a Monte Carlo
barostat, which will give exactly the correct NPT distribution for any
integrator.

Best,
Michael

On Fri, Jul 27, 2012 at 2:01 PM, Katie Maerzke  wrote:
> Hi all -
>
> Is there any plan to get LINCS working with the MTTK barostat in the near 
> future?
>
> Thanks
> Katie
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Only plain text messages are allowed!
> * Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] BAR gives different result than TI

[Michael Shirts]

The implementation of BAR in gromacs is pretty hard for me to follow
because of how everything is stored noncompactly in the histogram.  In
4.6, both can be computed from the same dhdl.xvg file, so it might be
easier to track down possible bugs.

On Fri, Jun 29, 2012 at 2:24 PM, David van der Spoel
 wrote:
> Hi,
>
> we've been trying to do free energy calculations for solvation of two small
> molecules in water, n-butylamine (NBA) and diethyl-ether (DEE). For one of
> them the result with BAR (using Justin's tutorial) is significantly
> different from TI:
>   BAR TI  Exper. (kJ/mol)
> NBA   -11.1   -10.8   -17.9
> DEE   -11.0   -1.8 -7.4
>
> Since for NBA the results are pretty close it seems that the protocol should
> be right, and the numbers seem to have converged pretty good as well.
>
> Any clues why these results could be so different?
>
> --
> David van der Spoel, Ph.D., Professor of Biology
> Dept. of Cell & Molec. Biol., Uppsala University.
> Box 596, 75124 Uppsala, Sweden. Phone:  +46184714205.
> sp...@xray.bmc.uu.sehttp://folding.bmc.uu.se
>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Only plain text messages are allowed!
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] What does --enable-fahcore mean?

[Michael Shirts]

It only matters for running on Folding@Home.  For other users of
gromacs, it doesn't do anything.


On Tue, Jun 26, 2012 at 3:50 PM, Bao Kai  wrote:
> Hi, all,
>
> I am wondering what the --enable-fahcore option of configure means.  I
> got the explanation from configure --help of "create a library with
> mdrun functionality", while it is not very clear to me.
>
> Best Regards,
> Kai
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Configurational bias Monte Carlo

[Michael Shirts]

There's a fair amount of interest for more general support for Monte
Carlo methods in GROMACS 5.0.  However, there is no any current
support for configuration bias Monte Carlo (or any other kind of MC)
currently in the code.

On Fri, Jun 1, 2012 at 12:10 PM, Benjamin Haley  wrote:
> Hello,
>
>   I have found the excellent documentation for creating polymer
> chains in GROMACS.  I want to create several chains and pack
> them into a volume, adjusting torsion angles to minimize
> interactions with other chains (and self-interaction within a chain).
> One method for doing this is the configurational bias Monte Carlo
> sampling.  Can this (or something similar) be done in GROMACS?
>
> Thank you!
> Ben Haley
> Purdue University
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] REMD question

[Michael Shirts]

Gromacs already supports replica exchange -- what particularly are you
implementing?

Equilibration of pressure is always a good idea -- even if you are
running NVT simulations, you want to get them to be at the equilibrium
volume for your system and temperature choice, which will require
equilibration at constant pressure.

On Mon, May 28, 2012 at 4:37 PM, Nathalia Garces  wrote:
> Dear Gromacs Users,
>
>
>
> We are implementing REMD method in Gromacs in protein folding, in your web
> page you give some steps that don´t mention any step about NPT
> stabilization.  This step is necessary to run REMD simulations?
>
>
>
> Thank you in advance,
>
>
>
> Nathalia
>
>
>
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] poor performance in Hemiltonian Replica Exchange

[Michael Shirts]

4.6 will include Hamiltonian replia exchange functionality built into
the MPI version.

Currently the description of the error is very vague -- if you can
write up what exactly the numbers are, and what they should be, with
files that exactly replicate the error, then I can take a look.  But
unless I can reproduce the error you are describing out of the box,
its unlikely I will be able to find it.

Additionally, it would be easiest if the files were deposited as a
redmine bug report, so that the information is centrally located.

Best,
Michael

On Thu, May 10, 2012 at 12:23 PM, francesco oteri
 wrote:
> Dear gromacs users,
>
> I performed a Hemiltonian Replica Exchange (i.e. replica exchange where each
> replica has a init_lambda=0, delta_lambda=0 and init_lambda ranging
> uniformely from 0 to 1).
>
> Since I have only ten fixed discrete lambda, I run a Temperature Replica
> Exchange where, for each replica I generated a .top file with the parameter
> rescaled through a
> python script ( in practice I did through python the same thing gromacs is
> supposed to do with the H-REM previously described). Now gromacs complained
> because
> every replica has the same setup, so I changed the temperatures using very
> close values (300.0001K,
>  300.0002K,300.0003K,300.0004K,300.0005K,300.0006K,300.0007K,300.0008K, 300.0009K)
> With this setup the simulation runs fine and I expect to have similar
> result.
>
> Then I compared the results observing two phenomena:
>
> 1) In the second case exchange rate is 100%, while in the first case I have
> an exchange rate close to 30%.
> Does it rise  because the temperatures are too close?
>
> 2) The second setup is 3x faster!
> In particular I observe an imbalance between PME and force calculation
> ranging from 10% to 60%.
> I tried to run each replia indipendently (a different mdrun instance for
> each .tpr file) but still I observe the same performance slowdown.
> I guess the free energy impairs the efficient force calculation, but I dont
> understand why.
>
> Can someone explain me the two observations?
>
>
>
> Francesco
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] BUG: Free energy calculation

[Michael Shirts]

Sabine, thanks for filing it in redmine!  Having a record helps a lot.
 Can you also attach all your input files to the redmine filing?  It
can only really be debugged if the input files you used are included.

Best,
Michael
On Fri, Mar 23, 2012 at 9:11 AM, Justin A. Lemkul  wrote:
>
>
> Sabine Reisser wrote:
>>
>> Hi Mark,
>>
>> with FE, without PR : same error
>> without FE, with PR: stable
>> without FE, without PR: stable
>>
>> I've never had this error before.
>>
>> Logfile says:
>> [...]
>> Initializing Domain Decomposition on 2 nodes
>> Dynamic load balancing: no
>> Will sort the charge groups at every domain (re)decomposition
>> Initial maximum inter charge-group distances:
>>    two-body bonded interactions: 3.341 nm, LJC Pairs NB, atoms 22476 22761
>>  multi-body bonded interactions: 0.649 nm, Angle, atoms 1668 1686
>> Minimum cell size due to bonded interactions: 3.675 nm
>> Using 0 separate PME nodes
>> Optimizing the DD grid for 2 cells with a minimum initial size of 3.675 nm
>> The maximum allowed number of cells is: X 1 Y 1 Z 2
>> Domain decomposition grid 1 x 1 x 2, separate PME nodes 0
>> PME domain decomposition: 2 x 1 x 1
>> Domain decomposition nodeid 0, coordinates 0 0 0
>> [...]
>> Linking all bonded interactions to atoms
>> There are 55376 inter charge-group exclusions,
>> will use an extra communication step for exclusion forces for PME
>>
>> The initial number of communication pulses is: Z 1
>> The initial domain decomposition cell size is: Z 5.09 nm
>>
>> The maximum allowed distance for charge groups involved in interactions
>> is:
>>                 non-bonded interactions           1.200 nm
>>            two-body bonded interactions  (-rdd)   4.213 nm
>>          multi-body bonded interactions  (-rdd)   4.213 nm
>>
>>
>
> These quantities look very weird to me.  They indicate interactions that are
> very far apart are influencing one another.  Can you provide a complete .mdp
> file?  It seems like some aspect of the free energy settings (perhaps
> couple-intramol?) and DD aren't getting along.  The other possibility is to
> try particle decomposition instead of DD (i.e. mdrun -pd).
>
> -Justin
>
> --
> 
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] BUG: Free energy calculation

[Michael Shirts]

Hi, Sabine-

If you can go to http://redmine.gromacs.org/projects/gromacs/issues
and file a bug report (including attaching files), I can look at it.
If that's ends up not working well, you can send me the files off of
the list, but it's usually better to have things in the redmine system
so problems are documented.  I'm also working on the updates to free
energy code in 4.6, so I want to make sure this will be solved there
as well.

Best,
Michael

On Fri, Mar 23, 2012 at 7:16 AM, Sabine Reisser  wrote:
> Hi Mark,
>
> with FE, without PR : same error
> without FE, with PR: stable
> without FE, without PR: stable
>
> I've never had this error before.
>
> Logfile says:
> [...]
> Initializing Domain Decomposition on 2 nodes
> Dynamic load balancing: no
> Will sort the charge groups at every domain (re)decomposition
> Initial maximum inter charge-group distances:
>    two-body bonded interactions: 3.341 nm, LJC Pairs NB, atoms 22476 22761
>  multi-body bonded interactions: 0.649 nm, Angle, atoms 1668 1686
> Minimum cell size due to bonded interactions: 3.675 nm
> Using 0 separate PME nodes
> Optimizing the DD grid for 2 cells with a minimum initial size of 3.675 nm
> The maximum allowed number of cells is: X 1 Y 1 Z 2
> Domain decomposition grid 1 x 1 x 2, separate PME nodes 0
> PME domain decomposition: 2 x 1 x 1
> Domain decomposition nodeid 0, coordinates 0 0 0
> [...]
> Linking all bonded interactions to atoms
> There are 55376 inter charge-group exclusions,
> will use an extra communication step for exclusion forces for PME
>
> The initial number of communication pulses is: Z 1
> The initial domain decomposition cell size is: Z 5.09 nm
>
> The maximum allowed distance for charge groups involved in interactions is:
>                 non-bonded interactions           1.200 nm
>            two-body bonded interactions  (-rdd)   4.213 nm
>          multi-body bonded interactions  (-rdd)   4.213 nm
>
>
>
>
>
> There it stops. In the 1-thread case, the second part is replaced by
> "Initiating Steepest Descents" and then writing out the energies for every
> step.
>
> Gromacs version is 4.5.5.
>
> I also attached the whole logfile.
>
> Cheers
> Sabine
>
>
>
>
>
> On 03/23/2012 11:52 AM, Mark Abraham wrote:
>>
>> On 23/03/2012 9:17 PM, Sabine Reisser wrote:
>>
>>>
>>> Dear gromacs users/developers,
>>>
>>> when trying to couple in a peptide into a membrane with:
>>>
>>> ; Define position restraints for peptide
>>> define          = -DPOSRES
>>>
>>> ; couple in peptide
>>> free_energy     = yes
>>> init_lambda     = 0.05
>>> sc_alpha        = 0.7
>>> sc_power        = 1
>>> couple-moltype  = Protein
>>> couple-lambda0  = none
>>> couple-lambda1  = vdw-q
>>>
>>>
>>> grompp works fine, but mdrun (2 threads) gives me
>>>
>>> Making 1D domain decomposition 1 x 1 x 2
>>> *** glibc detected *** mdrun: realloc(): invalid next size:
>>> 0x7f0f30305810 ***
>>>
>>> and breaks up.
>>>
>>>
>>> When running "mdrun -nt 1 " on only one thread, it works fine.
>>>
>>> Is this a known bug?
>>>
>>
>> First, is it likely not to be a problem with your setup... is your
>> system stable in parallel without FE code? Without position restraints?
>> What does your .log file say? What GROMACS version is it?
>>
>> Mark
>>
>
>
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Exchange interval in REMD

[Michael Shirts]

> It also depends whether or not you use constant pressure, in which case it
> makes sense to increase the time to long enough to let the volume relax.
> I still do not understand why people do NVT REMD, because it makes all but
> one replica have a pressure that is not the ambient pressure.

If all you care about is speedup at a single temperature, then this
fact doesn't really matter all that much.

Note that if you have an incorrect barostat (for example, Berendsen),
then NPT REMD is likely to produce some very bad artifacts.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] amber99sb in GROMACS vs amber99sb in AMBER

[Michael Shirts]

> Is anyone aware of any benchmark analysis about the implementation of the
> amber99sb force field (or any of its variants: 99sb-ildn, 99sb-nmr) in
> GROMACS and AMBER. I am interested to know in what extend the energies
> correlate and if the results agree with experimental data.

Whether the results compare agree with experimental data is irrelevant
to the correctness of the implementation -- that has to do with the
validity of AMBER99sb to begin with.  The only question is, do GROMACS
and AMBER give the same energies for the same configurations?

I actually do not know the answer, and would be interested to hear if
it's been tested.  http://ffamber.cnsm.csulb.edu/ does not appear to
have the very comprehensive tests that appear for earlier models at
the current time.  I SUSPECT there should not be a problem, since
AMBER99sb just changed a couple of torsions, so errors are unlikely
(though in theory possible).
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] How DispCorr influsnces pressure

[Michael Shirts]

> Dispcorr is only for homogeneous liquids. It should not be used for
> membranes.

More precisely -- the dispersion correction is an analytical
correction to both the pressure and energy that is rigorously correct
in the limit of the radial distribution function g(r)=1 outside the
van der Waals cutoff.  For a lipid bilayer, this not a good
approximation, since the system is not isotropic.

Since there is not "right" way to do membrane systems currently, the
best bet is to go with the same cutoff treatment the lipids were
parameterized with.

DispCorr = Ener is probably still a good idea, since g(r)=1 isn't
perfect, but it's better than completely ignoring the long range
energy.  DispCorr = Ener will not change the configurations sampled,
so any phase properties, pressures, etc. will be unaffected, but the
results will be somewhat less cutoff dependent.

See:
http://pubs.acs.org/doi/abs/10.1021/jp0735987
for some additional information.

Members of the Gromacs team are working behind the scenes for various
possible solutions to the problem for nonisotropic systems.  Don't
expect anything until 5.0, though.

Best,
Michael
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] FEP

[Michael Shirts]

> So BAR is only
> referring to the mathematical code used to calculate the overall free
> energy for the FEP, correct?

Yes.  The information required is the same, with the exception that
exponential averaging requires energy differences from only one state,
whereas BAR requires energy differences from both directions.  Since
you are likely running simulations at a range of states anyway,
there's no extra cost.

> My question is, for a mutation of a -CH3 group to a -H group, is it
> better to simply run:
> [+ from (Lambda=0 ,  R-CH3, full charges and interactions -STATE A)
> --> (Lambda=1, R-CH, full charges and interactions -STATE B)]
>
> OR
>
> [1) from (Lambda=0 ,  R-CH3, STATE A : Charges and LJ Interactions: ON)
> 2)  (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated)
> 3)  (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF)
> 4)  (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to 
> -H)
> 5)  (Lambda=1, R-CH3, -CH3  STATE B : Charges and LJ Interactions: ON)

Currently, you'd need to run three simulations to do this (will likely
be simplified somewhat in 4.6).

Set 1: turn R-CH3 charges off in a way that preserves the total charge.
Set 2: change CH3 LJ to H LJ
Set 3: Turn R-H charges on in a way that preserves the total charge.

Note that the first and last transformations will need only one two
states -- their energies will be VERY small.

> Reason I'm asking is because when I try the first choice to do it
> STATE A to STATE B in one step, when I reach Lambda=0.85 and above on
> the NVT equilibration right after EM, I receive errors saying that
> bonds are moving way to far off their constraints which leads me to
> believe that the system is moving too far from where it was energy
> minimized.  Errors such as:
>
> Step 188, time 0.376 (ps)
> LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.17, max 0.000636 (between atoms 9 and 68)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>      9     68   31.2    0.   0.1110      0.
>
> Step 188, time 0.376 (ps)  LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.15, max 0.000531 (between atoms 9 and 68)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>      9     68   31.0    0.   0.1110      0.

You probably have problems with the charge + LJ terms on the hydrogen.
 Note that you could possibly use soft core for the Coloumb
interactions for hydrogen -- it's sometimes causes sampling problems,
but for H's it probably doesn't matter.  I haven't done this approach
much so, so I can't tell for sure if it will work.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] FEP

[Michael Shirts]

FEP is a poorly defined term.  It can either mean 1) making small
changes 'alchemical' changes in the molecules and computing the free
energies by any method (BAR, TI, etc), or 2) it can mean specifically
computing the free energies by exponentially averaging the potential
energy differences.  Basically, using the exponential averaging
formula is a bad idea -- if you have code that only uses that method,
you can get decent results out if you are careful, but TI, BAR, and
MBAR are all better choices.

On Mon, Oct 10, 2011 at 10:58 AM, Justin A. Lemkul  wrote:
>
>
> Steven Neumann wrote:
>>
>> Thank you guys! So, is there any tutorial in Gromacs for calculating free
>> energy of ligand binding using FEP?
>>
>
> TI or BAR are better methods for calculating binding free energies, I would
> think.  FEP is best for mutating between different species.
>
> -Justin
>
>> Steven
>>
>> On Mon, Oct 10, 2011 at 2:02 PM, Justin A. Lemkul > > wrote:
>>
>>
>>
>>    mohsen ramezanpour wrote:
>>
>>        Hi
>>        Please have a look at Dr.Justin tutorial page at the following
>> link:
>>
>>
>>  http://www.bevanlab.biochem.__vt.edu/Pages/Personal/justin/__gmx-tutorials/index.html
>>
>>  
>>
>>
>>    This tutorial is not for FEP explicitly, but may be of some use.
>>     There is discussion on using the BAR algorithm for binding free
>>    energy calculations.
>>
>>    -Justin
>>
>>        Cheers
>>
>>
>>        On Mon, Oct 10, 2011 at 12:27 PM, Steven Neumann
>>        mailto:s.neuman...@gmail.com>
>>        >__>
>>        wrote:
>>
>>           Hi Gmx Users,
>>                Can you suggest some reading and some tutorial in
>>        calculations of
>>           binding free energy (ligand binding) in Gromacs? ?I want to
>>        use Free
>>           Energy Perturbation method.
>>                Steven
>>           --
>>           gmx-users mailing list    gmx-users@gromacs.org
>>        
>>           >
>>
>>           http://lists.gromacs.org/__mailman/listinfo/gmx-users
>>        
>>           Please search the archive at
>>           http://www.gromacs.org/__Support/Mailing_Lists/Search
>>         before
>>        posting!
>>           Please don't post (un)subscribe requests to the list. Use the
>>           www interface or send it to gmx-users-requ...@gromacs.org
>>        
>>           >        >.
>>
>>           Can't post? Read
>>        http://www.gromacs.org/__Support/Mailing_Lists
>>        
>>
>>
>>
>>    --     ==__==
>>
>>    Justin A. Lemkul
>>    Ph.D. Candidate
>>    ICTAS Doctoral Scholar
>>    MILES-IGERT Trainee
>>    Department of Biochemistry
>>    Virginia Tech
>>    Blacksburg, VA
>>    jalemkul[at]vt.edu  | (540) 231-9080
>>    
>>    http://www.bevanlab.biochem.__vt.edu/Pages/Personal/justin
>>    
>>
>>    ==__==
>>
>>    --     gmx-users mailing list    gmx-users@gromacs.org
>>    
>>    http://lists.gromacs.org/__mailman/listinfo/gmx-users
>>    
>>    Please search the archive at
>>    http://www.gromacs.org/__Support/Mailing_Lists/Search
>>     before posting!
>>    Please don't post (un)subscribe requests to the list. Use the www
>>    interface or send it to gmx-users-requ...@gromacs.org
>>    .
>>    Can't post? Read http://www.gromacs.org/__Support/Mailing_Lists
>>    
>>
>>
>
> --
> 
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Pleas

Re: [gmx-users] Question about Justin's Free Energy Tutorial

[Michael Shirts]

I don't think this is a question about new free energy code -- I think
this is asking about the fact if you can do a free energy calculation
by specifying the A and B variables in the topology, instead of using
the MDP coupl-moltype arguments.  This is actually the way free energy
calculations were managed before 4.0, and is still supported.

The answer is yes, though it requires some experience with what's
going on -- and I don't have time to document it all quite right now,
as it's described relatively thoroughly in the manual.

The one difference is that the coupl-moltype keywords makes it easy to
decouple molecules (turn off the intermolecular but not intramolecular
energies) from their environment, which is very hard (and sometimes
impossible) to do by changing A and B states. If you want to turn off
all interactions, then you can do it very straightforwardly by
changing the top files.

Best,
Michael

On Thu, Oct 6, 2011 at 2:42 PM, Justin A. Lemkul  wrote:
>
>
> Fabian Casteblanco wrote:
>>
>> Hello Justin,
>>
>> I have a question about your tutorial.  If I want to mutate one small
>> group of a molecule, I would have to not provide 'couple_lambda0' and
>> 'couple_lambda1', correct?  I would essentially have to follow sec
>> 5.7.4 in the Gromacs manual and I have to actually provide all state A
>> variable and all state B variables.  Gromacs would calculate the new B
>> state parameters for bond lengths, angles, etc, correct?  Are there
>> any other major differences to account for?
>>
>
> I have never attempted FEP with the new free energy code.  Try a simple test
> system first and make sure it works as expected.
>
> -Justin
>
> --
> 
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Recommended parameters for NVE simulation of SPCE water

[Michael Shirts]

Hmm.  Even now, I'm noticing problems with what I sent:

It should be;

> rcoulomb                    = 1.3.

For PME, rlist should equal rcoul.  PME-switch can improve energy
conservation, but the wrong PME-switch parameters can affect the
results too much.  PME w/o switch should be appropriate for most
'standard' usage.  In the next month or so (by 4.6 at least), the
manual will include some more guidance on accuracy in cutoffs.

Summarizing again:
> For pretty good energy conservation, I would suggest:
>
> rlist                           = 1.3
> coulombtype              = PME
> rcoulomb                    = 1.3
> vdw-type                    = Switch
> rvdw-switch                = 1.0
> rvdw                          = 1.1
>
> This should work quite well -- you might get some drift after 1-2 ns,
> but not much.  I'm working on developing suggested PME parameters
> right now for highly quantitative work, but it's not quite ready yet.
>
>
> 
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shi...@virginia.edu
> (434)-243-1821
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Hamiltonian replica exchange?

[Michael Shirts]

Hamiltonian replica exchange is planned for 4.6, and is being beta
tested by some users.


On Tue, Aug 16, 2011 at 2:39 PM, Sanku M  wrote:
> Hi,
>    I was wondering whether hamiltonian replica exchange simulation has been
> implemented in latest version of  gromacs . Or, is there any other way of
> performing the hamiltonian replica exchange using some variants of
> lambda-dynamics ?
> Sanku
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Recommended parameters for NVE simulation of SPCE water

[Michael Shirts]

> 1.  NOTE 1 above suggests that I use vdwtype = Shift.  When I do this, do
> you recommend that I apply long range dispersion corrections for both energy
> and pressure, using DispCorr = EnerPres, or for only energy, using DispCorr
> = Ener?  Typically, for various (non-NVE) calculations, I have been using
> DispCorr = no, but I am not sure if this is a good idea.  Pages 97-98 of the
> Gromacs 4.5.4 manual seem to suggest that the energy correction due to
> DispCorr is small and usually only significant for free energy calculations
> (which I will not be doing here).  As a rule of thumb, do you typically turn
> dispersion corrections off?

For constant pressure simulations, or for reaching the constant
pressure equilibrium simulation, you should definitely include a
dispersion correction -- the density will be too large, and will be
cutoff dependent.

For constant volume simulations, the dispersion correction will be
constant.  It will thus NOT affect energy conservation, but WILL
affect average potential energy and average total energy,
significantly.

> 2.  NOTE 2 above suggests that I use either coulombtype = PME-Switch or
> coulombtype = Reaction-Field-zero.  Do you have any advice or
> recommendation?

For pretty good energy conservation, I would suggest:

rlist   = 1.3
coulombtype  = PME
rcoulomb= 1.1
vdw-type= Switch
rvdw-switch= 1.0
rvdw  = 1.1

This should work quite well -- you might get some drift after 1-2 ns,
but not much.  I'm working on developing suggested PME parameters
right now for highly quantitative work, but it's not quite ready yet.



Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shi...@virginia.edu
(434)-243-1821
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] autocorrelation time of dVpot/dlambda?

[Michael Shirts]

> I am doing free energy calculation in Gromacs and want to get an error
> estimate of my results. Is it possible to compute the autocorrelation time
> of dVpot/dlambda in Gromacs using a certain length of trajectory such as 10
> ns? Thanks a lot

The amount of simulation time required to compute the autocorrelation
time of dVpot/dlambda will depend on the timescale of the system.  You
should be able to use g_analyze to analyze the energy.xvg output, and
get an autocorrelation time.  If your simulation time is 50x longer
than the autocorrelation time, then you are very likely converged.
However, if there are very long correlation times, you still might not
have captured all of the time dependent variations in  Vpot/dlambda.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Possible free energy bug?

[Michael Shirts]

Hi, all-

Have you tried running

constraints = hbonds?

That might eliminate some of the constraint issues.  Much less likely
for LINCS to break or have DD issues if only the hbonds are
constrained.  2 fs is not that big a deal for the heteroatom bonds.

Best,
Michael

On Thu, Mar 10, 2011 at 8:04 PM, Justin A. Lemkul  wrote:
>
> Hi Matt,
>
> Thanks for the extensive explanation and tips.  I'll work through things and
> report back.  It will take a while to get things going through (unless one
> of the early solutions works!) since I have no admin access to install new
> compilers, libraries, etc. and for some reason the only thing I can ever get
> to work in my home directory is Gromacs itself.  The joys of an aging
> cluster.
>
> We recently got access to gcc-4.4.5 on Linux, but we're stuck with 3.3 on OS
> X, so there's at least a bit of hope for one partition.
>
> Thanks again.
>
> -Justin
>
> Matthew Zwier wrote:
>>
>> Hi Justin,
>>
>> I should have specified that the segfault happened for us after we got
>> similar warnings and errors (DD and/or LINCS), so the segfault may
>> have been tangential.  Given that everything about your system worked
>> before GROMACS 4.5, it's possible that your older compilers are
>> generating code that's incompatible with the GROMACS assembly loops
>> (which you are likely running with, as they are the default option on
>> most mainstream processors).  The bug you mentioned in your original
>> post also has my antennae twitching about bad machine code.
>>
>> If that's indeed happening, it's almost certainly some bizarre
>> alignment issue, something like half of a float is getting overwritten
>> on the way into or out of the assembly code, which corruption would
>> trigger the results you describe.  It's also distantly possible that
>> GROMACS is working fine, but your copy of FFTW or BLAS/LAPACK (more
>> likely the latter) has alignment problems.  One final possibility
>> (which would explain the failure on YellowDog but unfortunately not
>> the failure on OS X) is that GCC is generating badly-aligned code for
>> auto-vectorized Altivec loops, which is still a problem for Intel's
>> SIMD instructions on 32-bit x86 architectures even with GCC 4.4.  I've
>> also observed MPI gather/reduce operations to foul up alignment (or
>> rigidly enforce it where badly compiled code is relying on broken
>> alignment) under exceedingly rare circumstances, usually involving
>> different libraries compiled with different compilers (which is
>> generally a bad idea for scientific code anyway).
>>
>> Okay...so all of that said, there are a few things to try:
>>
>> 1) Recompile GROMACS using -O2 instead of -O3; that'll turn off the
>> automatic vectorizer (on Yellow Dog) and various other relatively
>> risky optimizations (on both platforms).  CFLAGS="-O2 -march=powerpc"
>> in the environment AND on the configure command line would do that.
>> Check your build logs to make sure it took, though, because if you
>> don't do it exactly right, configure will ignore your directives and
>> merrily set up GROMACS to compile with -O3, which is the most likely
>> culprit for badly-aligned code.
>>
>> 2) Recompile GROMACS specifying a forced alignment flag.  I have no
>> experience with PowerPC, but -malign-natural and -malign-power look
>> like good initial guesses.  That's probably going to cause more
>> problems than it solves, but if you have a screwy BLAS/LAPACK or MPI,
>> it might help. I only suggest it because if you've already tried #1,
>> it will only take another half hour or hour of your time to recompile
>> GROMACS again.  Other than that, tinkering with alignment flags is a
>> really easy way to REALLY break code, so you might consider skipping
>> this and moving straight on to #3.
>>
>> 3) Snag GCC 4.4.4 or 4.4.5 and compile it, and use that to compile
>> GROMACS, again with -O2.  GCC takes forever to compile, but beyond
>> that, it's not as difficult as it could be.  Nothing preventing you
>> from installing it in your home directory, either, assuming you set
>> PATH and LD_LIBRARY_PATH (or DYLD_LIBRARY_PATH on OS X) properly.  You
>> might need to snag a new copy of binutils as well, if gcc refuses to
>> compile with the system ld.  This option would also probably get you
>> threading, since you certainly have hardware support for it.
>>
>> 4) Rebuild your entire GROMACS stack, including FFTW, BLAS/LAPACK,
>> MPI, and GROMACS itself with the same compiler (preferably GCC from
>> #3) and the same compiler options (which again should be -O2, and
>> definitely NOT any sort of alignment flag).  Put them in their own
>> tree (like "/opt/sci"), and definitely not in /usr (which is generally
>> managed by the system) or /usr/local (which tends to accumulate
>> cruft).  ATLAS is a good choice for BLAS, and there are directions on
>> the ATLAS website for building a complete and optimized LAPACK based
>> on BLAS.
>>
>> In practice, I've found I've had to do #4 for every piece of
>> scienti

Re: [gmx-users] free energy

[Michael Shirts]

One other thing I would point out is that the solvation free energy is
dependent on concentration.  you will get a different result with 4
polymer chains vs 3 vs 3, etc.  Make sure you understand the
dependence.  Also, the free energy will depend on the polymer chain
length.

Polymer and finite concentration calculations are harder to interpret
than monomer and infinite dilution calculations.  Make sure you
understand the differences.   I'm not sure understand all of them,
though, so you'll have to think about it yourself!  Basically, you
need to make sure the physical picture of the molecules in gromacs is
the physical picture of the realistic molecular system itself.


On Mon, Feb 14, 2011 at 6:28 PM, Moeed  wrote:
>
> Dear experts,
>
> I am going to do solvation FE of polymer (polyethylene) in a hydrocarbon
> solvent. I have prepared a system consisting of 4 polymer chains and 480
> hexane molecules with the actual density of polymer solution (~ 0.5 g/cm3).
>
> 1- For such a study I dont know how many polymers I need to have in my
> system. If FE can be done with only one chain, am I making system bigger in
> vain? Does this matter affect the accuracy of results?
>
> 2- I have switched off electrostatics so I am using
>
> free_energy  =   yes
> init_lambda  =   0
> delta_lambda =   0
> sc_alpha =   0.5
> sc-power =   1
> sc_sigma =   0.3
> couple-lambda0   =   vdw
> couple-lambda1   =   none
> couple-intramol  =   no
>
> In David Mobley's turorial the last three lines are not included. I wanted
> to know if I am to run say 10 simulations for different lambda, what purpose
> does the last three lines serve in 4.0.7  ? I got very close values in that
> tutorial without these settings. ( I know what these lines mean, just
> curious how these three lines affect the results in 4 X +).
>
> Please let me know your comments/point of view about the system and setting
> I am using.
>
> Thanks
> Moeed
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Secondary structure loss in implicit solvent simulations

[Michael Shirts]

OK, that is what I was trying to figure out -- is the problem
reproducible on both GPU and CPU.  Now, you havent answered the direct
question, if the energies are the same for at least the first 5 steps
are so -- without that knowledge, then here might be different errors
occurring in GPU vs CPU.

At this point, the question is then whether this works with another
code with the same input parameters.  Sander (part of the AmberTools
system)is downloadable (I believe to anyone, not just academics), so
you can try running the same system on Sander, using the best fit to
the implicit solvent model in gromacs.

If THAT works, and Gromacs GPU fails, then it would appear to be a
problem with Gromacs implicit solvent implementation, and should be
posted to redmine as a bug, as well as described on the list with full
details (more than you have provided so far!) so that it can be
reproduced by the developers.

Best,

On Tue, Jan 25, 2011 at 5:05 AM, K. Singhal  wrote:
> Hi
>
> It's not necessarily GPU-specific, it's implicit solvent-specific. I don't 
> get these problems in explicit solvent simulations on CPU, only in implicit 
> solvent simulations both on GPU as well as CPU. One of the problems that I 
> can think of is unbalanced charges that I would have balanced out using NaCl 
> ions, but not any more.
>
> Regards
> Kush
>
>
>
> --
> Kushagra Singhal
> Promovendus, Computational Chemistry
> van 't Hoff Institute of Molecular Sciences
> Science Park 904, room C2.119
> 1098 XH Amsterdam, The Netherlands
> +31 205256965
> Universiteit van Amsterdam
> k.sing...@uva.nl
>
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Re: Secondary structure loss in implicit solvent simulations

[Michael Shirts]

> Do you get the same effects if you run a "normal" simulation on CPU and not
> GPU?  That information would be critical for properly diagnosing what's
> going on. If it's not GPU-specific, in all likelihood whatever you're doing
> is incorrect somewhere along the way.

Folllowing up on this -- if it's not the same trajectory as the CPU to
within single precision, then something is wrong.

So the energies should be the same to 5 digits for the first time
steps, and then gradually drift away. If it's not the same for the
first 5-10 steps, then there is something going wrong with the GPU
code or the setup for the GPU code.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Secondary structure loss in implicit solvent simulations

[Michael Shirts]

A few questions:

1) What force field are you using?
2) do you get the same answers with and without GPU acceleration?
3) How long does it take for secondary structure to disappear?  100's
of ps?  10's of ns?


On Mon, Jan 17, 2011 at 8:45 AM, K. Singhal  wrote:
> Hi
>
> I have a question regarding the implicit solvent implementation in GROMACS 
> 4.5.1 (with GPU acceleration). I have been trying to simulate molecules with 
> varying size (PYP with 125 AAs, BBA5 with ~20 AAs, and and Trigger Factor 
> with 432 AAs), but can't seem to be able to maintain the secondary structure 
> in any of them. Can any one suggest a particular set of parameters than need 
> to be taken care of or used with certain fixed values for better results?
>
> While I have used a varied range of values for almost all parameters, here is 
> the list of parameters recently and unsuccessfully used:
>   integrator           = sd
>   nsteps               = 1000
>   init_step            = 0
>   ns_type              = Grid
>   nstlist              = 10
>   ndelta               = 2
>   nstcomm              = 10
>   comm_mode            = Linear
>   delta_t              = 0.002
>   pme_order            = 4
>   ewald_rtol           = 1e-05
>   ewald_geometry       = 0
>   epsilon_surface      = 0
>   optimize_fft         = FALSE
>   ePBC                 = xyz
>   bPeriodicMols        = FALSE
>   bContinuation        = FALSE
>   bShakeSOR            = FALSE
>   etc                  = No
>   nsttcouple           = -1
>   epc                  = No
>   epctype              = Isotropic
>   nstpcouple           = -1
>   tau_p                = 1
>   andersen_seed        = 815131
>   rlist                = 1.1
>   rlistlong            = 2.5
>   rtpi                 = 0.05
>   coulombtype          = Cut-off
>   rcoulomb_switch      = 0
>   rcoulomb             = 2.5
>   vdwtype              = Cut-off
>   rvdw_switch          = 0
>   rvdw                 = 2.5
>   epsilon_r            = 1
>   epsilon_rf           = 1
>   implicit_solvent     = GBSA
>   gb_algorithm         = OBC
>   gb_epsilon_solvent   = 80
>   nstgbradii           = 1
>   rgbradii             = 1.1
>   gb_saltconc          = 0.02
>   gb_obc_alpha         = 1
>   gb_obc_beta          = 0.8
>   gb_obc_gamma         = 4.85
>   gb_dielectric_offset = 0.009
>   sa_algorithm         = Still
>   sa_surface_tension   = 2.092
>   shake_tol            = 7.5e-06
>   bd_fric              = 0.5
>   ld_seed              = 2010
>
> Thanks and Regards
> Kush
>
>
>
> --
> Kushagra Singhal
> Promovendus, Computational Chemistry
> Universiteit van Amsterdam
>
> --
> gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] What does "Conserved-En." stands for

[Michael Shirts]

For integration schemes like Nose-Hoover that are not NVE but have
some other quantity that is conserved, then "Conserved-En. provides
the value of this conserved quantity.


Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shi...@virginia.edu
(434)-243-1821

On Mon, Nov 29, 2010 at 12:34 PM, Jim (Rui) Qiao  wrote:
> Dear all:
>
> When using g_energy to analyze the energies of the MD system, an entry
> termed "Conserved-En." is provided. What does this term stands for?
> The simulation is in NVT ensemble.
>
> Thanks,
>
> RQ
> --
> gmx-users mailing list    gmx-us...@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists