[gmx-users] version 3.3.1 .tpr file with version 3.2 g_order

[Michael]

   The g_order tool in gromacs 3.3.1 is broken.   I know that there 
is a patch available that would fix it, but a colleague advised me that 
recompiling g_order would not be quick and easy.  I therefore thought 
that I might use g_order from gromacs 3.2, but it complains of a version 
3.3.1 run input file with the version 3.2 g_order program.  The error 
message is:


Fatal error: reading tpx file (./topol.2.tpr) version 40 with version 31 
program


Is there any way that results from a simulation in gromacs 3.3.1 can be 
analyzed with an older version tool?  For example, would it work to run 
the version 3.3.1 .tpr file through the version 3.2 tpbconv program to 
get a new .tpr file, or would it yield the same error? 


Thanks for any suggestions
--Michael Skaug
University of California Davis
___
gmx-users mailing listgmx-users@gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to [EMAIL PROTECTED]

Can't post? Read http://www.gromacs.org/mailing_lists/users.php

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

[gmx-users] PMF from pull code, unexpected results

[Michael Brunsteiner]

Dear all,

I would like to calculate the potential of mean force between two molecules 
in aqueous solution using the pull code. For a start i performed a number 
of calulations with a couple of very simple model systems with settings
loosely based on the example in the gmx tutorial section (see mdp file included
at the bottom). My box is 8x8x8 nm**3, pmf is calculated from some r_min up to 
approx 3.8 nm. Performing simulations with two t-buntanol molecules in implicit
solvent, anlyzing the output (x and f*.xvg files) with g_wham, i get rather
unexpected results:

The PMFs look very differnt depending on whether I use "pulldim  Y Y Y" 
or "pulldim  "N N Y":
http://www.brunsteiner.net/tbut-pmf.gif
Since we look at two neutral molecules with a permanent dipole
moment one would expect the PMF to be negative up to a distance
at which the two molecules come into contact, but with "pulldim  Y Y Y"
the PMF is positive, i.e., repulsive, throughout.

With "pulldim  N N Y" I constrain the two molecules in two
dimensions (by freezing the central atom in the x and y dimensions only,
BTW any ideas how else i could achieve that?)
One could argue that a combination of frozen atoms and pull code might be 
problematic, but i freeze and pull in different dimensions, so this should
be OK, and, more importantly: with "pulldim  N N Y" i get results that
look much more reasonable than with  "pulldim  "Y Y Y" and NO additional 
constraints.

With "pulldim  Y Y Y" the RDF (option nolog in g_wahm) looks, in fact, 
like a NON-normalized rdf:
http://www.brunsteiner.net/tbut-rdf.gif
and, interestingly if i normalize it myself (by dividing through z**2 before
taking the logarithm) the resulting PMF looks much more reasonable.

Although what I see suggests that with "pulldim  Y Y Y" the RDF just doesn't
normalized properly, this issue seems to be more involved: 
Looking at the average forces and positions (the content of the f*xvg and x*xvg 
files)
http://www.brunsteiner.net/tbut-f.gif
http://www.brunsteiner.net/tbut-z.gif
suggests that there's already something wrong in the mdrun output
meaning that the problem is further upstream and not connected to anything
done by g_wham.

I repeated the calculations with an even simpler system (2 single atoms
that only interact via a LJ potential) to get qualitatively the same
results:
http://www.brunsteiner.net/2at-pmf.gif
http://www.brunsteiner.net/2at-rdf.gif
http://www.brunsteiner.net/2at-z.gif
http://www.brunsteiner.net/2at-f.gif

A few more remarks:
1) Convergence does not seem to be an issue here as i extended
some of the calculations to include 35+ windows, 2 ns each, and
the PMF remains the same nearly quantitatively.
2) The length of my reaction coordinates is always shorter
than half the box length.
3) I've compared calculations cut-off vs PME, to get very similar
results.
4) If I use "pulldim  Y Y Y" with no additional constraints but with
"comm_mode = Angular" i get results somewhere inbetween the two cases
above.

my specs:
Gromacs version: 4.5.3
Linux 2.6.35-23-generic Ubuntu x86_64
gcc version 4.4.5

I am not sure whether i overlooked something in my input,
or whether there's a problem with code.
I'd be grateful for any ideas/suggestions!

cheers,
Michael



=
mdp file:

integrator  = sd; this is better than md/NHT for systems with very 
few DOFs
tau_t   = 2.0 2.0
ref_t   = 290 290   ; a bit lower than 298 since sd with default 
parameters
; has problems controlling the temperature.
tc_grps = p1 p2 ; also tried System here, no difference
;
dt  = 0.002
tinit   = 0
nsteps  = 10; played with that, results seem to have converged
nstxout = 1000
nstvout = 0
nstfout = 1000
nstenergy   = 100
;
constraint_algorithm = lincs
constraints  = hbonds
continuation = no
;
comm_mode   = Linear
nstcomm = 1
pbc = xyz   ; also tried "pbc = no" same result
nstlist = 10
ns_type = grid 
rlist   = 4.0
rcoulomb= 4.0
rvdw= 4.0
;
coulombtype = Shift
vdwtype = Shift
epsilon_r   = 80
;
pull= umbrella
pull_geometry   = distance
pull_dim= N N Y   ; or "Y Y Y"
pull_start  = yes
pull_ngroups= 1
pull_group0 = p1
pull_group1 = p2
pull_rate1  = 0.0
pull_k1 = 500 ; i let this number vary depending on the pull 
distance;
; also played with the numbers, same result
pull_nstxout= 100
pull_nstfout= 100
pull_pbcatom0   = 1   ; my system geometry is so that the molecules 
never leave the box
pull_pbcatom1   = 16  ; or cross a cell boundary.
;
freezegrps  = c1 c2; with pulldim Y Y Y these line

[gmx-users] Re: PMF from pull code, unexpected results

[Michael Brunsteiner]

chris.neale at utoronto.ca wrote:

> Looking at http://www.brunsteiner.net/tbut-pmf.gif you seem to be getting 
>exactly what I would
> expect. The difference is the entropy term. Note that the spherical shell 
>increases volume as r
> increases for pulldim YYY but this effect is absent for pulldim NNY. This is 
>why you can correct
> as per an RDF.

The "entropy term" (Eqn 6.1 in the manual) has already been the source
of some confusion in this list ... 

If you make a "Gedankenexperiment" and consider the PMF between two atoms in 
vacuum 

that interact ONLY through a simple radial, for example a LJ, potential, and 
then calculate the PMF 

using the pull-code ...  if you don't do umbrella sampling + WHAM but instead 
use constraints 

(pull = constraint) you can do the calculation in your head, to find that the 
resulting PMF 

is, of course, nothing else but the original LJ potential. And you will get the 
SAME result
if you do this calculation in one or in three dimensions (pulldim NNY or 
pulldim 
YYY)
so where does the entropy and the correction term in eqn 6 come in here??

In Section 6.1 of the manual it says: "But when calculating a PMF between two 
solutes in a
solvent, for the purpose of simulating without solvent, the entropic 
contribution should be removed."
I find this a bit confusing, first ... "simulating without solvent" ... why 
would that be
my (or anybodies) purpose? and second:

according to eqn 6.1 this "entropic contribution" is: PMF(r) = - (nc-1)kBT 
log(r)

for this to make sense r would have to be dimensionless wouldn't it?
I could divide r by the distance at which i arbitrarily chose my PMF to vanish,
call it r_c, which would have the advantage that then the correction is zero at 
this
distance ... but then this is just wild guess of mine ... anyway, that would 
mean 

that, if I call the PMF coming out of mdrun/g_wham PMF_wham, 
then the "true" PMF is: PMF(z) = PMF_wham (z) + (nc-1)kBT log(z/r_c)
(the plus comes from the double minus, as in: "removing" something thats 
negative)

is that correct? ... and if so, why does it, seemingly, not apply in the above 
Gedankenexperiment?

Yours truly (and maybe some other people) would really appreciate if somebody 
in 
the know
would clarify that!

thanks in advance!

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

[gmx-users] PMF from pull code, unexpected results

[Michael Brunsteiner]

> Michael: your mdp options indicated using US and my stance is that,   > using 
>US, you are incorrect to say that "And you will get the SAME   > result if you 
>do this calculation in one or in three dimensions   > (pulldim NNY or pulldim 
>YYY)".

this has nothing to do with the mdp file ... i was talking about a 
purely hypothetical case, if you use pull = constraints, rather than
umbrella sampling you don't NEED to run the calculation, since, as i said
you can figure out the result in your head ... with constraints the "average"
force will always be the exact LJ force at that distance, and if you integrate
that you will necessarily end up with the original LJ potential, no
matter whether you do that in one or in three dimensions!
if you only have two spherically symmetric particles, with an interaction
energy that depends only on the distance, and no solvent then the 
"free energy" and the "energy" are identical (apart from the
entropic contribution discussed in the Neumann paper) ...

that means: if you COULD do this calculation with constraints (in 
practice you probably can't because there'd be too few DOFs and the temperature
would be ill defined) you SHOULD get identical results, independent of pulldim
... unless the term in eqn 6.1 is somehow included explicitly, but i can't see
that anywhere in the code ... 
> Nevertheless, I'm a little perturbed by 
> your call out for help from somebody else "in the know", so I'll leave   > 
> you 
>at this point. Good luck.

thanks, doesn't look as if there was much interest, I guess I'll have to go 
back to my statistical mechanics textbooks...

mic


  
-- 
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] topology files for ligands in MD Simulation under OPLS AA force field

[Michael Brunsteiner]

mircial at sjtu.edu.cn mircial at sjtu.edu.cn wrote

> I am using GROMACS package to do molecular dynamics
> simulations under OPLS_AA force field. I encounter
> some problems when preparing the topology files
> of small molecules (ligands).My questions are as follows: 
> 1, how to chose the atom type of each atom from the ligands? 
amber + GAFF would be easier, here the ligand parameterization
is pretty straight forward (there's probably a HOWTO)

with OPLS you want to work with homology, and if you don't need
any high throughput you can do it manually ... look up 
groups in amino acids that are similar to your the components
of your ligands, and choose the atomtypes accordingly.

> 2, how to define the charges of each atom? I know that, when
> the atom type is defined, there will be a corresponding charge
> to this atom type. Does it safe to use this charge in the
> simulations on the ligands (since these charges are designed
> for amino acids)?

yes, OPLS was mostly used for peptides, but AFAIK, the
original parameterization in Jorgensens group started with generic
organic molecules, there's certainly papers from this group that
present parameters for non-peptide molecules, e.g., various
hetero cycles ... (I am not sure whether the gmx topology folder
only contains the part of OPLS relevant for peptides, or all of it,
if not ask someone from the Jorgenson group, they might give you
the original parameter files)
However, if you have some generic organic molecules and assign
OPLS atom types there is no guarantee that the charges add up to
zero (for a neutral molecule) ... what i did before, with
reasonable success, is combine EPSD (SCF or AM1BCC) charges with
OPLS parameters...

In any case, as has been pointed out before many times in this list
mixing different FFs, or using ad-hoc values for some parameters, is 
dangerous, and you really want to test your molecules before using 
them in any extensive calculations. 
A good start might be to compare structures that were minimized
with your FF to the same structure minimized at the DFT or MP2 level
(HF-SCF might do as well depending on the molecule) ... but that's
probably not enough ...

> 3,When the atom type and the charge of atom are defined,
> how to prepare the file in the GROMACS format? Does there easier
> method to prepare such files than manually? 
As for the bonded interactions, in some cases pdb2gmx and
grompp can figure them out by themselves. You'd assign atom- 
(thereby bond-) types, and the bonds, then make a new entry
in the rtp file, including only the atoms and bonds sections. you 
choose a new 3 letter residue name, and starting from a PDB
file of the ligand (with this residue name, and the same atom order
and the same atom names as in the rtp entry) use pdb2gmx/grompp.
However often enough, especially for proper and improper dihedrals
you have to define them manually (again by anlogy) 
for anything but the simplest molecules its a painstaking task ...
so you either take the time, or you use amber (or you talk to 
Schrodinger, they will have some script for assigning OPLS parameters
but that probably wont work with Gromacs)

good luck!
mic


  
-- 
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
>

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

[gmx-users] thermostats, canonical distribution, and stochastic dynamics

[Michael Brunsteiner]

Dear all,

I encountered problems when trying to keep
a simple system at constant temperature, and would
appreciate any help/advice!

I try to calculate a PMF between two relatively simple
particles (about spherical, each 60 atoms) in solution with
umbrella sampling. For now i model the solvent as primitive
model electrolyte (vacuum + epsilonr=80), and what i see is
this:

trying out the following thermostats:
* leapfrog + nose hoover
* velocity verlet + 10 nose hoover chains
* v-rescale
in ALL cases, i get what seems to be
a non-canonical distribution:

Since the distance between the two particles is restrained
via a harmonic potential one would expect (at least for large
distances) to get an approximately Gaussian distribution of
the distance between the two particles, centred near the
equilibrium distance of the harmonic restraint. However, I 
keep seeing a bi-modal distribution (as one would expect 
from a harmonic oscillator with no thermostat) which, I believe,
can be taken as evidence for a non-canonical phase space
coverage (and probably this also screws up the WHAM
analysis)

My system is relatively simple - pathologically so one might
argue - but then, what I recall from articles I read, especially
about Nose Hoover chains, it seems that I could expect getting
a proper canonical distribution even for much simpler systems
with no more than a few degrees of freedom (and my system
has >300 degrees of freedom)

With Langevin dynamics (instead of md) the distance
distribution DOES become Gaussian, but now I encounter
a different problem: More so than with MD convergence
becomes an issue: especially with small tau_t (<5 ps). 
E.g. to cover a distance of about 3 nm it takes at least 44 windows
with >=1 ns each to get a reasonably converged PMF. So I wonder:
i) is there any limit to the size of tau_t I am supposed to use? or 
ii) is there any other way to accelerate convergence?

FTR: I also tried using two tc_grps (one for each particle) 
but, at least qualitatively, i see the same results.

thanks in advance for any help!

regards,
Michael



===
mdp-file:

integrator  = sd   ; or md, or md-vv-avek
tau_t   = 0.2 0.2  ; varies from 0.2 to 10.0 for sd
ref_t   = 300 300
tc_grps = p1 p2; or just System
dt  = 0.001
tinit   = 0
nsteps  = 40   ; or more
nstxtcout   = 0
nstxout = 1000
nstvout = 0
nstfout = 0
nstenergy   = 100
constraints = h-bonds
comm_mode   = Linear
nstcomm = 1
pbc = no
nstlist = 0
ns_type = simple
rlist   = 5.0
rcoulomb= 5.0
rvdw= 5.0
coulombtype = Cut-off
vdwtype = User  ; repulsive only LJ pot
epsilon_r   = 80
pull= umbrella
pull_geometry   = distance
pull_dim= N N Y
pull_start  = no
pull_ngroups= 1
pull_group0 = p1
pull_group1 = p2
pull_rate1  = 0.0
pull_k1 = kkk; varies from 500-5000
pull_nstxout= 100
pull_nstfout= 100
pull_pbcatom0   =  61
pull_pbcatom1   = 122
pull_init1  = ddd ; varies from 1-4
freezegrps  = c1 c2
freezedim   = Y Y N Y Y N ; this is to keep the particles in 1-D
energygrps  = p1 p2


  
-- 
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] IMPLICIT SOLVENT IN AMBER

[Michael Brunsteiner]

I recently came across a paper (J. Phys. Chem. B 2005, 109,
10441-10448) in which  the authors found that using implicit
solvent with explicit ions you easily end up with all the
ions cuddling up in clusters ... which is obviously an artefact.

That is, if you do use this combination you certainly need to
carefully choose your parameters and system set-up, and
run thorough tests before doing any production runs.

... my five cents
mic




Per Larsson wrote:

Hi!

The answer is that I do not now critical it is.
I have seen some papers that seem to hint at it not being overly critical, but, 
again, life does not come at a net charge.
You'd have to see for yourself in your system, given the observables that are 
important I suppose.
Possibly we can implement this in the future, as there indeed seems to be some 
interest in it, but that is not a task I can deal with currently.

/Per


10 mar 2011 kl. 10.52 skrev Yulian Gavrilov:

> Thanks again!
> Don't you know how to make a total charge = 0 in this case, if implicit salt 
>concentration is not implemented currently? Or it is not critically? 
>
> 
> 2011/3/10 Per Larsson 
> Hi!
> 
> Yes, except that in point 2, I'm not sure about the effects of explicit ions 
> in 
>an implicit solvent. 
>
> Do deal with that properly one should use an implicit salt concentration, but 
>that is not implemented currently. 
>
> The choice of water-model with pdb2gmx is not important. You can choose 
> 'None' 
>here.


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

[gmx-users] pull forces

[Michael Brunsteiner]

Dear all,
does anybody know what the forces that are saved
in the f*xvg files from the pull-code actually are? are these:

1) the (sum of the) forces due to the normal non-bonded interactions
in the system acting on the ref and the the full groups.
2) only the forces from the harmonic restraint.
3) the sum of both?

i just did a test writing out the forces from the trajectory
file, and looking at the relevant forces (the ones in the z-dimension
in my case, as i use pulldim N N Y and constrain the groups to stay
put in the x and y dimensions) it seems as if the forces on the
pull group are oscillating as expected, but the forces on the
reference group are close to, but not exactly(!) zero.
i am not sure how to interpret this ...

cheers
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

[gmx-users] mdrun

[michael zhenin]

Dear all,
I am trying to run dynamics (mdrun) on GROMACS 4.5.3 in parallel on 8
processors, but it crashes after a while and refuses to reach to the end.
The error note that pops out is:

1 particles communicated to PME node 4 are more than 2/3 times the cut-off

>> out of the domain decomposition cell of their charge group in dimension
y.

>> This usually means that your system is not well equilibrated


My system contains protein, ATP and Mg. I have equilibirated it for 300 ps
that were divided into 3 parts.
The first 2 parts were the NVT equilibiration (the first one was made in
order to equilibrate only the water) - 100 ps each.


The third was the NPT - 100 ps.
The system wasn't constraint.

My mdp file is :

title = OPLS wt nbd1 MD
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 2000 ; 1 * 2000 = 2 ps, 20 ns


dt = 0.001 ; 1 fs
; Output control
nstxout = 2 ; save coordinates every 2 ps
nstvout = 2 ; save velocities every 2 ps
nstxtcout = 2 ; xtc compressed trajectory output every 2 ps
nstenergy = 2 ; save energies every 2 ps


nstlog = 2 ; update log file every 2 ps
; Bond parameters
continuation = yes ; Restarting after NPT
; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 5 ; 10 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)


rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
ewald_rtol = 1e-05


fourierspacing = 0.12 ; grid spacing for FFT
pme_order = 6 ;
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate


tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors


tau_p = 1.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC

; Dispersion correction

DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
morse = yes


I've searched the net for this error and didn't find any answer. Did anybody
see this error? and what should i do with it..

Thanks,

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

[gmx-users] Re: Is there still interest in rigid-body simulation?

[Michael Hagan]

Ignacio,

I would be very interested in a rigid body capability.

On 3/28/2011 5:45 AM, gmx-users-requ...@gromacs.org wrote:

Is there still interest in rigid-body simulation?

--
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] speed issue, gmx runtime estimate = f(t)

[Michael Brunsteiner]

Dear all, 

I am trying to optimize runtimes/speed for a series of
calculations i plan to do with a rather complex system.
When looking at runtimes I observed one peculiar issue:
the estimated runtime written to stderr by mdrun with the "-v"
option keepsgrowing ... in my experience this estimate normally
converges to a fairly accurate value within a few hundred
steps. not so here ..
for example, starting a job at 12:14 gromacs writes to stderr:

step 100, will finish Wed Apr  6 12:40:35 2011vol 0.97  imb F  1% 
step 200, will finish Wed Apr  6 12:39:42 2011vol 0.95  imb F  2%
step 300, will finish Wed Apr  6 12:39:47 2011vol 0.95  imb F  1%
step 400, will finish Wed Apr  6 12:39:52 2011vol 0.96  imb F  1%
...
step 38400, will finish Wed Apr  6 14:09:11 2011vol 0.91  imb F 29% 

so the estimated run time grows steadily from 25 to 110 minutes 
within about one hour...
this is mdrun_d (gmx-4.5.3) running on a linux PC with 
2 processores/4 nodes/8 threads. There's a couple of other programs (X, editor, 
...)
running at the same time but none of them uses a lot of CPU.

Anyway, what i see lets me expect there might be a way 
to speed up these calculations. I'd appreciate any suggestions!
(below is my mdp file.)

cheers,
Michael


mdp file:

integrator   = md
;
dt   = 0.001
tinit= 0
nsteps   = 10
nstxtcout= 0
nstxout  = 1000
nstvout  = 0
nstfout  = 0
nstenergy= 1000
;
constraint_algorithm = lincs
constraints  = hbonds
;
comm_mode= Linear
nstcomm  = 1
pbc  = xyz
nstlist  = 10
ns_type  = grid 
rlist= 1.4
rcoulomb = 1.4
rvdw = 1.4
;
coulombtype  = PME
vdwtype  = User
fourierspacing   = 0.2
optimize_fft = yes
;
Tcoupl   = nose-hoover
tc_grps  = System
ref_t= 300.0
tau_t= 0.4
gen_vel  = yes
gen_temp = 240.0
gen_seed = 22
;
pull = umbrella
pull_geometry= distance
pull_dim = N N Y
pull_start   = yes
pull_ngroups = 1
pull_group0  = f1
pull_group1  = f2
pull_rate1   = -0.015
pull_k1  = 1000
pull_nstxout = 100
pull_nstfout = 100
pull_pbcatom0= 241
pull_pbcatom1= 482
;
freezegrps   = c1 c2
freezedim= Y Y N Y Y N
;
energygrps   = SOL CC CF
energygrp_table  = CF CF  CF SOL  CF CC
-- 
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] PME

[Michael Brunsteiner]

Elisabeth,

You CAN, in fact calculate the contribution of the reciprocal part
of the PME energy to the binding energy between two components in 
a heterogeneous system, its just quite tedious...
say, your system is molecules A and B for which you want to know
the interaction energy, and the rest of the system, typically
the solvent, we call C.
Now your total Reciprocal Coulomb energy will have six parts: 
ER_tot = ER_AA + ER_BB + ER_CC + ER_AB + ER_AC + ER_BC
but these parts are NOT given in the gromacs output as they
cannot be calculated DIRECTLY, you have to calculate
them by setting the charges on A, B, or C (or combinations thereof)
to zero (there is a tool for setting the charges in a tpr file
to zero) and then do more runs with: "mdrun -rerun" based on the 
original trajectory to get the required contributions.

then E_AB = ER_C0 - ER_A0C0 - ER_B0C0

(or something like it, do double check that formula, i can't be bothered
thinking it through now ... here ER_A0C0, for example,  is the reciprocal
part of the coulomb energy with charges in groups A and C set to zero, etc)

this being said ... it's tedious, time-consuming, and error-prone
(you need to use double precision and save a lot of frames to
get reasonably accurate numbers)
You might be better off using reaction field, or PME and simply
ignore the reciprocal part altogether (if your molecules A, B
are NOT charged and have no permanent and large dipole moment
you might get away with the latter)

What Justin said is correct, PME (or any other Ewald-like
method, PPPM, FMA, etc) is standard these days, and for a good reason.
However, different properties are affected to a different
extent by neglecting the long range interactions, and for 
what you want to calculate it might be OK for getting at least
a qualitative answer, as long as you use PME for the actual MD.
(I'd be VERY surprised if everybody who did LIE in the last 10
years went through the trouble outlined above)

have fun!

mic





Elisabeth wrote:
> Hello Justin,
> 
> Several days ago you answered my question about calculating nonbonded 
> terms:
> 
> Question: If I want to look at nonboded interactions only, do I have to 
> add  Coul. recip.  to [ LJ (SR)  + Coulomb (SR) ] ?
> 
> Answer: The PME-related terms contain both solute-solvent, 
> solvent-solvent, and potentially solute-solute terms (depending on the 
> size and nature of the solute), so trying to interpret this term in some 
> pairwise fashion is an exercise in futility.
> 
> my question is if I want to add up nonbonded related terms to get inter 
> molecular energies, do I have to add Coul. recip. or it is already 
> included in Coulomb (SR)? 
> 

They are separate energy terms.  The PME mesh terms is "Coul. recip." and the 
short-range interactions (contained within rcoulomb, calculated by a modified 
switch potential) are "Coulomb (SR)."

> and also, for a A-B system, I have been using energy groups to extract  
> solute-solvent, solvent-solvent, solute-solute terms. Did you mean that 
> applying doing so with PME as electrostatics treatment is not correct?
> 

PME has been consistently shown to be one of the most accurate long-range 
electrostatics methods and is widely used, but in your case is preventing you 
from extracting the quantity you're after (if it can even be reasonably defined 
at all).  Using energygrps will not resolve the problem I described above.  The 
"Coul. recip." term contains long-range energies between (potentially) A-B, 
A-A, 

and B-B, depending on the nature of what A and B are.  The only terms that are 
decomposed via energygrps are the short-range terms, which are calculated 
pairwise.  Thus, with PME, there is no straightforward way to simply define an 
"intermolecular energy" for a heterogeneous system.  You might be able to 
define 

such a term for a completely homogeneous system (which also assumes that the 
sampling has converged such that the charge densities etc are uniform...but I'm 
sort of thinking out loud on that), but not one that is a mixture.

-Justin

> Thanks for your help!
> Best,
> 
> 
> 

-- 

-- 
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] PME

[Michael Brunsteiner]

Elisabeth wrote:

> The point is in my case there is no option other than ignoring LR since LR
> is not covered by shift or switch functions but at least what PME reports
> for SR is more accurate.

not really ... PME electrostatic energy is a sum of two contributions
LR and SR ... and only the SUM is meaningful, therefore if you use
a switching function you might, in fact, get something thats MORE
accurate than using ONLY the SR part of PME.

> So the decomposed Coulmb. SR terms I am getting
> using energy groups from PME are "reliable ?

see above ... but then again, the PME (or Ewald) SR part
is in fact nothing but the normal 1/r function scaled by an
exponential term so that it goes smoothly to zero at the
cut-off radius, and in that sense you can regard it as the
total Coulomb interaction scaled with yet another tapering
function (like shift or switch)
If you want to make sure to get precise energies, I'd do
the procedure I laid out in my previous mail for a couple
of simple model coumpounds and than compare the energies
with what you get from PME-SR and/or shift/switch for
the same systems (or the same trajectories using
"mdrun -rerun")
And, mind you, here only the RELATIVE energies are
important, not the absolute numbers, as the latter might
differ a lot between different methods, but if your relative
numbers are OK then you should be fine.
by that i mean: if you use two methods X and Y to get
interaction energies (E) for a series of compounds A, B, C, ...
then you should find that (approximately)
E_AB_X/E_CD_X = E_AB_Y/E_CD_Y etc ...
if the two (or more) terms differ a lot than (at least) one
of your methods X or Y gives you the wrong answer.

> BTW: I am dealing with non polar particles i.e alkanes and carbon and
> hydrogen are the only species I have.

does this include the solvent? or ist that water? if so that
won't help, but if you really only have hydrocarbons you will probably get away
with cut-off/shift/switch ...

Anyway, what i would do is this: search for papers of people
who did LIE calculations (linear interaction energy, a method for
approximate binding free energies) and see what they did ... 
Jorgenson, Aqvist, ...

If you are in a rush, just use reaction field electrostatics,
its probably good enough for your purposes.

mic


  
-- 
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] sugar force fields

[Michael Brunsteiner]

Dear all,

I understand that most major force fields (charmm, opls, amber, gromos?...) 
come with parameters for sugars. I wonder:

1)  if any work has been done/published comparing
the accuracy of these different sugar force fields

2) Of the force fields that come with gromacs is there one
for which making sugar topologies is as straight forward as making 
topologies for amino acids, or if there is no such thing for which FF is 
making a sugar topology the least tedious?

BTW: I want to model simple sugar molecules, monomers,
dimers (no polymers)

thanks!

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

[gmx-users] unexpected results with stochastic dynamics

[Michael Brunsteiner]

Hi everybody,

I just ran a number of simulations with two large and roughly spherical 
molecules
in a primitive model electrolyte (epsilon=80), they are constrained to the 
z-axes and
i calculate the PMF using the pull code restraining the distance between the
two molecules. integration is done with langevin dynamcs (gmx-4.5.3)

what i see is this: I start several runs all with the same initial orientations 
of the two
molecules only the distance varies, the dynamics seems to go well, but when i 
compare

the FINAL ORIENTATIONS of the two molecules after multiple runs at different
distances they are NEARLY THE SAME not matter which distance i choose.

below is the relevant part of my mdp file:


integrator   = sd
tau_t    = 2.0
ref_t    = 300
gen_vel  = yes
gen_temp = 240.0
gen_seed = xxx ; this varies
tc_grps  = System
nsteps   = 20
dt   = 0.001

The forces depend on the distance of course (the molecules have
a heterogeneous charge distribution) and this is a highly chaotic system
therefore the probability that the conformations (basically the
rotational DOFs of the two molecules) are the same or very similar after
200 ps at 300K should basically be ZERO ... but this is what i see, 

its re-producible and happens at different distances AND with different
molecular topologies (charge distributions).


a possible explanation for this behaviour might be that of the terms 

in eqn. 3.108 only the velocity dependent and the random term are actually
considered but not force term AND gen_seed doesn't make any difference for 
SD ... then the small differences i see would be only from round-off errors ... 
???

has anybody seen something similar ?

i'd appreciate any suggestions!

thanks
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] implicit solvent LINCS errors

[Michael Daily]

The following mdp file produces a successful dynamics run out to 100K 
steps / 200 ps.  What I discovered is this:


Using the md integrator, it is necessary to turn off pressure coupling.  
However, pressure coupling works with sd (Langevin) integrator.


Mike

---

; title and include files
title= 1EX6-S35P_md1
;cpp  = cpp
;include  = 
-I/home/yoo2/myusr/gromacs-3.3.3/share/gromacs/top/ -I./

define   =
; integrator and input/output setting up
;integrator   = md
integrator = sd

nsteps   = 100 ; 2 ns
dt   = 0.002
nstxout  = 5000
nstvout  = 5000
nstenergy= 500
nstxtcout= 500
nstlog   = 500
xtc_grps = System
energygrps   = System
comm_mode= Linear

;implicit solvent
implicit_solvent = GBSA
gb_algorithm = OBC
gb_saltconc = 0.15 ; note - this feature is not functional yet
rgbradii = 2.0 ; need larger cut-off than atomistic models

; neighbor searching and vdw/pme setting up
nstlist  = 10
ns_type  = grid
pbc  = xyz
;rlist= 1.4
rlist= 2.0

;coulombtype  = pme
coulombtype  = Cut-Off
fourierspacing   = 0.1
pme_order= 6
;rcoulomb = 1.4
rcoulomb = 2.0

;vdwtype  = switch
vdwtype  = Cut-Off
rvdw_switch  = 1.0
;rvdw = 1.2
rvdw = 2.0

; cpt control
tcoupl   = v-rescale
tc-grps  = System
tau_t= 0.4
ref_t= 300.0
Pcoupl   = parrinello-rahman
pcoupltype   = isotropic
tau_p= 1.0
compressibility  = 4.5e-5
ref_p= 1.0

; velocity & temperature control
gen_vel  = yes
gen_temp = 300.0
annealing= no
constraints  = hbonds
constraint_algorithm = lincs
morse= no

---


On 5/30/11 8:45 PM, Michael D. Daily wrote:
Thanks for the recommendations everyone.  I tried all of the mdp 
changes recommended by Justin (increase rlist, rvdw, etc to 2.0; 
change T-coupling to v-scale, and eliminate P-coupling).  When I 
increased the distance cutoffs, it ran about 30 ps then crashed 
instead of crashing immediately.  Only when I also turned off 
P-coupling did it keep running long-term (now it's up to 500K+ steps / 
1 ns), so apparently the problem was a lack of control of system 
volume.  Reducing the timestep from 2 fs to 1fs did not work by itself 
and was not necessary ultimately to get it to run.


Josh, thanks for the physical explanation of why the protein might 
explode.  My starting structure was a minimized xtal structure so I 
wouldn't call it "extended."  I'll try again with 4.5.4 before I scale 
this up any further.


Mike


On 5/30/2011 8:02 PM, Mark Abraham wrote:

On 31/05/2011 10:54 AM, Justin A. Lemkul wrote:



Joshua L. Phillips wrote:
I've found that I often get LINCS warnings like this when starting 
from

highly extended conformations when using implicit solvent. The GBSA
surface tension combined with the lack of viscosity (due to the 
absence
of explicit water) allows the protein to change conformation much 
faster

than LINCS likes by default.

Normally, I just run a short vacuum simulation (keep the same settings
as Justin suggested but set GBSA = No) to let the system relax a 
little

more before starting a GBSA run (EM is often just not enough). Usually
only about 50ps. Also, doing this with position restraints can help 
slow
down the collapse, and usually results in just enough collapse that 
the

following GBSA run will satisfy the default LINCS settings.

Another option might be to run a short simulation with GBSA turned on,
but using Langevin dynamics to include some additional friction that
should slow down the initial collapse.



That's generally a good idea.  We often use the sd integrator when 
doing GBSA calculations.



Also, you could fudge with the environment variables associated with
LINCS, but this seems a little dangerous compared to the above two
suggestions.



I wouldn't say "a little dangerous," I'd say "very dangerous" :)  If 
the constraints are failing, it's not necessarily (or usually) their 
fault.  The system is unstable, and trying to just override this 
will give you a trajectory, but you should be concerned that this 
trajectory is being ruled by spurious forces.



I would be interested in peoples' opinions and suggestions on how to
handle this issue as it is quite common when starting from highly
extended structures using GBSA. Some proteins are more susceptible 
than

others...

Re: [gmx-users] oplsaa vs. charmm

[Michael Daily]

  
  
Do you have some experimental data to compare to your IDP
simulations, like X-ray scattering or some such? I'd imagine that
IDP simulations with either forcefield would only be qualitatively
accurate given that the forcefields are calibrated, as you say, on
rigid proteins and small molecules.

On 6/8/11 8:00 AM, Thomas Evangelidis wrote:
Dear Prof van der Spoel and GROMACS users,
  
  I have read the article "Scrutinizing Molecular Mechanics Force
  Fields..." where it is demonstrated that AMBER99sb yields protein
  conformations that are in better agreement with residual dipolar
  coupling, J-coupling and NOE data, compared with other force
  fields. However, both proteins used for this benchmark study
  (ubiquitin and protein G) are rather rigid, so I was wondering if
  there is a similar analysis for flexible proteins/peptides. I want
  to simulate a few intrinsically disordered proteins/peptides of
  length between 40 and 100 aa and would like to know what would be
  the best choice of force field to use. Any experience or knowledge
  on this matter would be greatly appreciated!
  
  thanks in advance,
  Thomas
  
  
  
  On 27 May 2011 19:01, David van der Spoel
<sp...@xray.bmc.uu.se>
wrote:

  
On 2011-05-27 17.50, simon sham wrote:
  
Hi,
I have recently done two simulations on a protein at
high temperature
near its melting temperature. At the beginning I used
CHARMM forcefield,
and then OPLSAA to double check the results. There are
some differences
in the structures between the forcefield used. I know
different
forcefields can give different results. All parameters
in the
simulations were the same except the forcefield. Is
there anyway I can
tell which forcefield gives more reliable results?

Thanks for the insights,

Simon

  

  
  You might want to check
  
  Oliver Lange, David van der Spoel and Bert de Groot:
  Scrutinizing Molecular Mechanics Force Fields on the
  Microsecond Timescale With NMR Data Biophys. J. 99 pp. 647-655
  (2010)
  
  where we compare a number of FFs to NMR data.
  

-- 
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.se    http://folding.bmc.uu.se
  

  -- 
  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

  

  
  
  
  
  -- 
  ==
  Thomas Evangelidis
  PhD student
  Biomedical Research
Foundation, Academy of Athens
  4 Soranou Ephessiou ,
115 27
Athens, Greece

email: tev...@bioacademy.gr
    teva...@gmail.com
  
website:
https://sites.google.com/site/thomasevangelidishomepage/
  
  
  



-- 
Michael D. Daily, Ph.D.
Postdoctoral Fellow
Qiang Cui group
Department of Chemistry
University of Wisconsin-Madison
  

-- 
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] AFM Simulation

[Michael McGovern]

Hi everyone.  I'm doing a simulation of a system where a peptide is linked to a 
surface and an AFM tip is brought in contact to do force measurements.  

The AFM tip moves very slowly (1nm/ms), so I'm simulating stationary tips at 
various distances.  Does anyone know the most direct measurement from the 
simulations that would correspond to the force measured in the experiments? 
  
Thanks a lot.
-- 
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] Programs to add residues

[Michael Lerner]

On Tue, Jun 21, 2011 at 4:12 PM, Chris Neale wrote:

> Try "Loopy". You can get it to build termini in addition to loops.
>
> http://wiki.c2b2.columbia.edu/**honiglab_public/index.php/**Software:Loopy<http://wiki.c2b2.columbia.edu/honiglab_public/index.php/Software:Loopy>
>
>
I've also seen people use PyMOL's builder to do this. Either way, you'll
need to take (a lot of) extra care to minimize your linker and make sure
that it looks reasonable.

-Michael


> Nevertheless, I'd suggest simply omitting that part of the protein and
> capping your new terminus to remove the charge. You will have more
> difficulties converging the conformation of the unstructured terminus than
> you may expect.


> CHris.
>
> -- original message --
>
> Hi all,
>
> Is there a program that allows the user to add residues to the N and C
> terminus, without using the electron density?  I would like to add a
> short linker to my protein which doesn't exist in the electron
> density.
>
> Thanks a lot,
>
> Sincerely,
> Zack
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
> Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Search<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<http://www.gromacs.org/Support/Mailing_Lists>
>



-- 
Michael Lerner, Ph.D.
IRTA Postdoctoral Fellow
Laboratory of Computational Biology NIH/NHLBI
5635 Fishers Lane, Room T909, MSC 9314
Rockville, MD 20852 (UPS/FedEx/Reality)
Bethesda MD 20892-9314 (USPS)
-- 
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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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] 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]

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]

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] 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] 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] 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

[gmx-users] tabulated potentials and soft core free energy - this should not work ....

[Michael Brunsteiner]

Hi,

i meant to perform a free energy calculation to get the chemical
potential of water, and made a number of input files based on the excellent 
tutorial
by Justin Lemkul.
(Without thinking) I also tried using a tabulated potential for electrostatics
with: coulombtype = User
but it seems to work - all the numbers in dhdl.xvg are finite and non-zero, and 
neither mdrun nor
grompp give me a warning) 
but then it shouldn't work ... how is Gromacs supposed to calculate the 
dU/dlambda term ?
or Is this calculated numerically?
or did i overlook something here?
below my mdp file

thanks for any advice!

michael



mdp:


integrator   = sd
nsteps   = 10
dt   = 0.002
;
pbc  = xyz
nstcomm  = 100
comm_grps    = System
;
nstxtcout    =
0
nstenergy    = 100
nst_log  = 0
nstxout  = 0
nstvout  = 0
;
rlist    = 1.35
;
coulombtype  = User
rcoulomb = 1.2
vdw-type =
switch
rvdw-switch  = 1.0
rvdw = 1.2
;
tc_grps  = system
tau_t    = 1.0
ref_t    = 300
;
Pcoupl   = Berendsen
tau_p    = 0.5
compressibility  =
4.5e-5
ref_p    = 1.0
pcoupltype   = isotropic
;
constraints  = hbonds
lincs_iter   = 8
;
free_energy  = yes
init_lambda  = 0.50
delta_lambda = 0
foreign_lambda   = 0.40 0.60
dhdl_derivatives =
yes
sc_alpha = 0.5
sc_power = 1
sc_sigma = 0.5
nstdhdl  = 100
couple-moltype   = w1
couple-lambda0   = vdw-q
couple-lambda1   = none
couple-intramol  = no



===
Why be happy when you could be normal?
--
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] g_bar ... the right answer for the wrong reason?

[Michael Brunsteiner]

Hi,
i just performed a free energy TI calculation, to get the
free energy of solvation of water in water (the chemical potential of water)
i stuck closely to the templates given in the tutoral by justin lemkul.
The final result i get with g_bar looks good, and the number is within the 
error-bars
of both experiment and literature data - still i keep getting a warning ...
i first turn off the charges and then LJ ... but for the LJ part (only) g_bar 
keeps saying:

WARNING: Some of these results violate the Second Law of Thermodynamics: 
 This is can be the result of severe undersampling, or (more likely)
 there is something wrong with the simulations.

i can't think of anything that's wrong with my simulations, nor do i think
this is undersampling (21 windows each 2 nano seconds)
also the difference between my mdp and justin's are minor (as i believe, see 
below)
why do i get the correct answer but still this warning keeps coming?

at lambda=0 my LJ interactions are fully turned on, and for the DeltaG(0->0.05)
i,e, the first row in the table below i get from g_bar:

 lam_A  lam_B  DG   +/- s_A   +/- s_B   +/-   stdev   +/- 
 0.000  0.050   -2.90  0.03    2.15  0.03 9647088805553410.00 
9570465417587940.00    3.30  0.05
 0.050  0.100   -3.36  0.03    1.48  0.02    1.47  0.03    1.95  0.02
 0.100  0.150   -0.49  0.00    0.03  0.00    0.03  0.00    0.22  0.00
 0.150  0.200    0.12  0.00    0.00  0.00    0.00  0.00    0.05  0.00
etc...

there appears to be a problem with the phase space overlap, but then
all other values, and the end-result, look ok ... can it be that there is an 
issue
with the gmx implementations of soft core potentials and g_bar?


cheers
mic


mdp:
;
integrator   = sd
nsteps   = 100
;dt   = 0.002

pbc  = xyz
nstcomm  = 100
comm_grps    = System
;
nstxtcout    = 0
nstenergy    = 100
nst_log  = 100
nstxout  = 0
nstvout  = 0
;
rlist    = 1.0
;
coulombtype  =
PME
rcoulomb = 1.0
vdw-type = switch
rvdw-switch  = 0.8
rvdw = 0.9
;
DispCorr = EnerPres
;
fourierspacing   = 0.12
pme_order    = 6
ewald_rtol   = 1e-06
epsilon_surface    = 0
optimize_fft = no
;
tc_grps  = system
tau_t    = 1.0
ref_t    = 300
;
Pcoupl   = Berendsen
tau_p    = 0.5
compressibility  = 4.5e-5
ref_p    = 1.0
pcoupltype   = isotropic
;
constraints  = hbonds
lincs_iter   = 8
;
free_energy  = yes
init_lambda  = 0.00
delta_lambda = 0.0
foreign_lambda   = 0.05 0.10
dhdl_derivatives = yes
sc_alpha = 0.5
sc_power = 1
sc_sigma = 0.5   ; i also tried 0.3 same problem
nstdhdl  = 100
couple-moltype   = w1
couple-lambda0   = vdw
couple-lambda1   = none
couple-intramol  = no




===
Why be happy when you could be normal?

--
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] problem with git/4.6

[Michael Brunsteiner]

hi

as suggested in: http://www.gromacs.org/Developer_Zone/Roadmap/GROMACS_4.6

i do:

prompt> git clone git://git.gromacs.org/gromacs.git
Cloning into 'gromacs'...
remote: Counting objects: 119131, done.
remote: Compressing objects: 100% (21428/21428), done.
remote: Total 119131 (delta 100642), reused 115488 (delta 97560)
Receiving objects: 100% (119131/119131), 46.91 MiB | 2.75 MiB/s, done.
Resolving deltas: 100% (100642/100642), done.

and then:

prompt> git checkout --track -b release-4-6 origin/release-4-6
fatal: Not a git repository (or any of the parent directories): .git

if cd to gromacs directory first:

prompt> cd gromacs
prompt> git checkout --track -b release-4-6 origin/release-4-6
Branch release-4-6 set up to track remote branch release-4-6 from origin.
Switched to a new branch 'release-4-6'
prompt> git pull
Already up-to-date.

i am not sure if this is what i want ... the source i have there is still 4.5.5 
or is it not?
... nothing seems to have been downloaded after the first command ...
how do i get the source of the 4.6 development branch?
is there a way to simply get a tar-ball ?

thanks,
mic

-- 
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] problem with git/4.6

[Michael Brunsteiner]

> You should have it.  In CMakeLists.txt, PROJECT_VERSION should be set to 
> "4.6-dev" so you can check that.

ithat what i looks like ... i now get:

prompt> mdrun_d
[...]
:-)  VERSION 4.6-dev-20120629-9c6be1c  (-:
[...]

which gives me:
g_bar_d -f mdv*.xvg -b 100
[...]
total   0.000 -  1.000,   DG -9.00 +/-  0.15

with vanilla 4.5.5 i get:
g_bar_d -f mdv*.xvg -b 100
[...]
total   0.000 -  1.000,   DG -9.00 +/-  0.15

there's two points:

1) this value is fine, as, once combined with the coulomb part, it is very 
close to literature and
exptl data. however, in both cases  get this warning ...  Second Law of 
Thermodynamics etc...
and riduculously high s_B values at lambda=0.

2) the reason i saw different DeltaG values earlier was that i had included (by 
mistake)
the dhdl file from a single additional window at lambda=0 that was produced 
with exactly the same
input apart from the value of 0.3 (instead of 0.5) for sc_sigma ... as far as i 
understand things (not very
well perhaps) the resulting free energy difference should NOT depend on this 
parameter ... so neither should
the result of g_bar ... but maybe it does when using foreign lambdas as i did 
here ? if this is so
it might be a good idea making g_bar read tpr files to give a warning in such a 
case ...

i guess for now i'll wait for a stable 4.6 - any ideas when this will be there?

thanks!
michael
--
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] BAR gives different result than TI

[Michael Brunsteiner]

Hi David,


i saw your recent post in the gmx mailing list about "BAR gives different 
result than TI" 
i wonder did you get any good answers to this queston so far!?

i was also trying to do some free energy calculations, with BAR i seem to get 
the correct
answer (compared to expt and literature) but there still seems to be some 
issues -
see: http://lists.gromacs.org/pipermail/gmx-users/2012-June/072867.html


also, i assume you used soft core potentials ... i remember quite a while ago i 
compared
results for free energies of solvation using  either sc_power=1 or 2 ... and 
the results
of (as i am sure, converged) calculations differed significantly. but this was 
a while ago
and in the meantime things might have changed.

another point: you write that your "numbers seem to have converged pretty good" 
...

the molecules you use are obviously flexible and i recall talking once to a 
this guy
from some software company, they had done extensive FE calculations and he told
me that their main problem was getting converged results for the gas-phase leg 
of
the simulations, in your mail you only gave the net-results for the two 
molecules,
how do the two separate calculations (mutation in vacuum and mutation in 
solvent)
compare?

cheers
michael


===
Why be happy when you could be normal?

--
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] water with atom-based cut-off ...

[Michael Brunsteiner]

hi,

i'd like to do a simulation of a solute in water - trying to reproduce some 
literature data
i try to be consistent and use an atom-based cut-off. for this each atom needs 
to be
its own charge-group ...
so when i change spc.itp
from:
 1  opls_116   1    SOL OW  1  -0.82
 2  opls_117   1    SOL    HW1  1   0.41
 3  opls_117   1    SOL    HW2  1   0.41
to:
    1  opls_116   1    SOL OW  1  -0.82
    2  opls_117   1    SOL    HW1  2   0.41
    3  opls_117   1    SOL    HW2  3   0.41

grompp complains that for settle all atoms need to be in the same charge group.
so i removed the settle and exclusion statements from spc.itp and instead
put the commands


constraints   = h-angles
constraint_algorithm  = LINCS
lincs_order    = 8
lincs_iter   = 2


in the mdp file ... but now constraints are not applied to waters at all (O and 
H look as if they move independently)
it also didn't help renaming the water molecule, residue and atom names ...

any suggestions on how i can apply atom-based cut-off with lincs or shake to 
water?

thanks,

michael


===
Why be happy when you could be normal?
--
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] tricky business

[Michael Brunsteiner]

In June 2005 Berk Hess wrote (in this list):

> Determining ionic solvation free energies is very tricky business,
> experimentally as well as theoretically.
> Last week I heard a talk by P. H. Hünenberger who explained
> all the proper corrections you need to make for different electrostatics
> treatments to obtain the correct energy. For PME several corrections
> are required (most of which are not in the code). 

are they NOW? if so what are they? ... and which are the corrections that were 
already implemented
in 2005? (i didn't find anything in the manual, neither in terms of methods nor 
references)

I don't know what P. H. Hünenberger said in his talk, but i am mostly concerned
about correction terms discussed here: J. Phys. Chem. 1996, 100, 1206-1215

thanks,
michael



===
Why be happy when you could be normal?
--
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] 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

[gmx-users] NMA and g_nmeig segmentation fault

[Michael Howard]

Hello All -- I'm trying to do a normal modes analysis of a fairly
large crystal system (~20,000 atoms) in double-precision GROMACS. The
system is xy periodic with 2 9-3 walls, and has periodic molecules. I
first minimized the structure with L-BFGS to 1e-5 tolerance with
switched vdW interactions and PME electrostatics (3dc geometry).

When I tried to run g_nmeig on the matrix, I received the following
error message from the shell

/var/spool/pbs/mom_priv/jobs/2480011.lionxj.rcc.psu.edu.SC: line 7:
3205 Segmentation fault  g_nmeig_d -f wall3_nm -s wall3_nm -T 300
-first 1 -last 2 -of wall3_efreq -ol wall3_eigval -v wall3_eigvec
>&nm.g_nmeig

The last thing written to the g_nmeig output is
Reading file wall3_nm.tpr, VERSION 4.5.5 (double precision)
Reading file wall3_nm.tpr, VERSION 4.5.5 (double precision)
Reading double precision matrix generated by Gromacs VERSION 4.5.5

I was able to successfully run NMA using the same code on a similar
system that was only 1/3 of the current size (reaction-field-zero
electrostatics) and also on a system 1/3 of the size with xyz
periodicity and PME electrostatics, and received no segmentation fault
error. I'm a little bit baffled about the source of the segmentation
error, because I would have expected it to crash during the Hessian
calculation if there were a problem with allocating system memory or
essentially singular forces. Has anyone else encountered a similar
problem?

Thanks in advance for your help.

-- Mike Howard
-- 
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] 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] 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

[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: 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

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: 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:
>>> >>

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] 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

[gmx-users] Charge groups in Charmm for lipids

[Michael McGovern]

Hi, I'm making a lipid topology for DPPG in the charmm forcefield in gromacs 
4.5.  I'm putting it together from pieces of what already exists as suggested 
in the charmm files.  I noticed that there are charge groups in the original 
charmm files, and these were reflected in the gromacs topology in the 4.5 beta 
release, but in 4.5.3 and 4.5.4, every atom has its own charge group.  

So for instace, POPC in the beta release had a topology that began:
[ POPC ]
 [ atoms ]
        N       NTL     -0.60   0
        C11     CTL2    -0.10   0
        C12     CTL5    -0.35   0
        C13     CTL5    -0.35   0
        C14     CTL5    -0.35   0
        H11     HL      0.25    0
        H12     HL      0.25    0
        H21     HL      0.25    0
        H22     HL      0.25    0
        H23     HL      0.25    0
        H31     HL      0.25    0
        H32     HL      0.25    0
        H33     HL      0.25    0
        H41     HL      0.25    0
        H42     HL      0.25    0
        H43     HL      0.25    0
        C15     CTL2    -0.08   1
        H51     HAL2    0.09    1
        H52     HAL2    0.09    1
        P1      PL      1.50    2
        O3      O2L     -0.78   2
        O4      O2L     -0.78   2
        O1      OSL     -0.57   2
        O2      OSL     -0.57   2
        C1      CTL2    -0.08   3
...
While in 4.5.4 this is
[ POPC ]
 [ atoms ]
        N       NTL     -0.60   0
        C11     CTL2    -0.10   1
        C12     CTL5    -0.35   2
        C13     CTL5    -0.35   3
        C14     CTL5    -0.35   4
        H11     HL      0.25    5
        H12     HL      0.25    6
        H21     HL      0.25    7
        H22     HL      0.25    8
        H23     HL      0.25    9
        H31     HL      0.25    10
        H32     HL      0.25    11
        H33     HL      0.25    12
        H41     HL      0.25    13
        H42     HL      0.25    14
        H43     HL      0.25    15
        C15     CTL2    -0.08   16
        H51     HAL2    0.09    17
...
The charge groups do not have integer charge and are all one atom.  The beta 
version has groups that agree with the charmm groups and have integer charge.  
Was this changed on purpose? The beta version seems correct to me, but am I 
missing something as to why this change would be made?  
Thanks-- 
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] Charge groups in Charmm for lipids

[Michael McGovern]

Thanks a lot, I will check that out.



From: Thomas Piggot 
To: Discussion list for GROMACS users 
Sent: Friday, August 5, 2011 6:48 PM
Subject: Re: [gmx-users] Charge groups in Charmm for lipids


 Hi,

Yes, this was done on purpose to reflect the different use of charge
groups in CHARMM and GROMACS. Check out the mailing list archives
for more details about this. 

On another note you might like to take a look at the contributed
CHARMM36 force field (available to download on the GROMACS website).
There is a DPPG entry in the lipids.rtp file that should help in
making a CHARMM27 DPPG rtp entry (if you wish to use the CHARMM27
lipids rather than the CHARMM36 ones).

Cheers

Tom 

On 06/08/11 00:22, Michael McGovern wrote: 
Hi, I'm making a lipid topology for DPPG in the charmm forcefield in gromacs 
4.5.  I'm putting it together from pieces of what already exists as suggested 
in the charmm files.  I noticed that there are charge groups in the original 
charmm files, and these were reflected in the gromacs topology in the 4.5 beta 
release, but in 4.5.3 and 4.5.4, every atom has its own charge group.  
>
>
>So for instace, POPC in the beta release had a topology that began:
>[ POPC ]
> [ atoms ]
>        N       NTL     -0.60   0
>        C11     CTL2    -0.10   0
>        C12     CTL5    -0.35   0
>        C13     CTL5    -0.35   0
>        C14     CTL5    -0.35   0
>        H11     HL      0.25    0
>        H12     HL      0.25    0
>        H21     HL      0.25    0
>        H22     HL      0.25    0
>        H23     HL      0.25    0
>        H31     HL      0.25    0
>        H32     HL      0.25    0
>        H33     HL      0.25    0
>        H41     HL      0.25    0
>        H42     HL      0.25    0
>        H43     HL      0.25    0
>        C15     CTL2    -0.08   1
>        H51     HAL2    0.09    1
>        H52     HAL2    0.09    1
>        P1      PL      1.50    2
>        O3      O2L     -0.78   2
>        O4      O2L     -0.78   2
>        O1      OSL     -0.57   2
>        O2      OSL     -0.57   2
>        C1      CTL2    -0.08   3
>...
>While in 4.5.4 this is
>[ POPC ]
> [ atoms ]
>        N       NTL     -0.60   0
>        C11     CTL2    -0.10   1
>        C12     CTL5    -0.35   2
>        C13     CTL5    -0.35   3
>        C14     CTL5    -0.35   4
>        H11     HL      0.25    5
>        H12     HL      0.25    6
>        H21     HL      0.25    7
>        H22     HL      0.25    8
>        H23     HL      0.25    9
>        H31     HL      0.25    10
>        H32     HL      0.25    11
>        H33     HL      0.25    12
>        H41     HL      0.25    13
>        H42     HL      0.25    14
>        H43     HL      0.25    15
>        C15     CTL2    -0.08   16
>        H51     HAL2    0.09    17
>...
>The charge groups do not have integer charge and are all one atom.  The beta 
>version has groups that agree with the charmm groups and have integer charge.  
>Was this changed on purpose? The beta version seems correct to me, but am I 
>missing something as to why this change would be made?  
>Thanks

-- 
Dr Thomas Piggot
University of Southampton, UK. 
-- 
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] 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]

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

[gmx-users] Using genbox to REMOVE solvent

[Michael Daily]

Hi all,

I am looking for a relatively easy way to make minor modifications to a
solvated system in gromacs without having to replace the whole solvent
layer.  Specifically, I'm swapping out some molecules (technically, just
swapping residue and atom names) in a MARTINI model and I think this causes
clashes with the first solvation shell if I increase the size of the
headgroup (e.g. POPE -> POPC).  So what I want to do is simply find any
solvent molecules that have severe clashes after this modification and
remove them - I know genbox checks for overlaps in this way when initially
solvating a system, so is it possible to make it check for and remove such
overlaps in an EXISTING system?

Thanks,
Mike
-- 
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] Using genbox to REMOVE solvent

[Michael Daily]

I tried an experiment similar to yours - after converting the POPE's to
POPC's when I remove all waters that are near the (lipid-emulsified)
nanoparticle based on a radial criterion and then run genbox again, the
total number of waters decreases by about 1100.  Compared to a control case
(ran genbox on the starting structure) that added 1600 waters, this means a
net loss of about 2700 waters.

In any case this is an important discussion for any of you running MARTINI
models - yes it's straightforward to "convert" lipid molecules but the
result can be a terrible starting structure - fortunately the system can be
regenerated from scratch with a new lipid config overnight.

On Wed, Aug 17, 2011 at 11:54 PM, Mark Abraham wrote:

> On 18/08/2011 4:32 PM, Michael Daily wrote:
>
>> Hi all,
>>
>> I am looking for a relatively easy way to make minor modifications to a
>> solvated system in gromacs without having to replace the whole solvent
>> layer.  Specifically, I'm swapping out some molecules (technically, just
>> swapping residue and atom names) in a MARTINI model and I think this causes
>> clashes with the first solvation shell if I increase the size of the
>> headgroup (e.g. POPE -> POPC).  So what I want to do is simply find any
>> solvent molecules that have severe clashes after this modification and
>> remove them - I know genbox checks for overlaps in this way when initially
>> solvating a system, so is it possible to make it check for and remove such
>> overlaps in an EXISTING system?
>>
>
> No.
>
> One way to test your hypothesis is to strip the waters from POPC, resolvate
> it with the right atom size, change it to POPE, solvate the
> previously-solvated structure and see if more waters get added. Now they'll
> be in a contiguous lump at the end of the .gro file.
>
> Mark
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
> Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Search<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<http://www.gromacs.org/Support/Mailing_Lists>
>



-- 

Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] Using genbox to REMOVE solvent

[Michael Daily]

And I am happy to report that this method produces a viable starting
structure :)

On Thu, Aug 18, 2011 at 7:02 PM, Michael Daily wrote:

> I tried an experiment similar to yours - after converting the POPE's to
> POPC's when I remove all waters that are near the (lipid-emulsified)
> nanoparticle based on a radial criterion and then run genbox again, the
> total number of waters decreases by about 1100.  Compared to a control case
> (ran genbox on the starting structure) that added 1600 waters, this means a
> net loss of about 2700 waters.
>
> In any case this is an important discussion for any of you running MARTINI
> models - yes it's straightforward to "convert" lipid molecules but the
> result can be a terrible starting structure - fortunately the system can be
> regenerated from scratch with a new lipid config overnight.
>
>
> On Wed, Aug 17, 2011 at 11:54 PM, Mark Abraham wrote:
>
>> On 18/08/2011 4:32 PM, Michael Daily wrote:
>>
>>> Hi all,
>>>
>>> I am looking for a relatively easy way to make minor modifications to a
>>> solvated system in gromacs without having to replace the whole solvent
>>> layer.  Specifically, I'm swapping out some molecules (technically, just
>>> swapping residue and atom names) in a MARTINI model and I think this causes
>>> clashes with the first solvation shell if I increase the size of the
>>> headgroup (e.g. POPE -> POPC).  So what I want to do is simply find any
>>> solvent molecules that have severe clashes after this modification and
>>> remove them - I know genbox checks for overlaps in this way when initially
>>> solvating a system, so is it possible to make it check for and remove such
>>> overlaps in an EXISTING system?
>>>
>>
>> No.
>>
>> One way to test your hypothesis is to strip the waters from POPC,
>> resolvate it with the right atom size, change it to POPE, solvate the
>> previously-solvated structure and see if more waters get added. Now they'll
>> be in a contiguous lump at the end of the .gro file.
>>
>> Mark
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
>> Please search the archive at http://www.gromacs.org/**
>> Support/Mailing_Lists/Search<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<http://www.gromacs.org/Support/Mailing_Lists>
>>
>
>
>
> --
> 
> Michael D. Daily
> Postdoctoral research associate
> Pacific Northwest National Lab (PNNL)
> 509-375-4581
> (formerly Qiang Cui group, University of Wisconsin-Madison)
>



-- 

Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] surface area calculations for MARTINI

[Michael Daily]

Hi,

I'm trying to do some surface area calculations on MARTINI models using
APBS, which requires supplying vdw radii via a pqr file.  Given how MARTINI
treats LJ interactions (0.47 nm radius for all particle types), how do I
calculate reasonable effective vdw radii? Obviously the PC headgroup (NC3)
will have a larger vdw radius than the PE headgroup (NH3), and g_sas seems
to account for this (reported ASA for -NC3 group is always higher than for
-NH3 group) but I can't figure out how.  Indeed, using "editconf -mead" to
produce a pqr file from a MARTINI structure shows 2.35 A as the vdw radius
for all lipid groups, which I fear will cause inaccurate calculations in
APBS.  Has anyone had this problem before?

Mike

-- 
====
Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] MARTINI / all-atom mapping

[Michael Daily]

Hi all,

I am trying to reverse-map some martini lipids to united atom.  In order to
do this, I'd prefer to have an EXACT definition of the aa-to-cg mapping.  I
cannot find this, only an imprecise graphic, in the MARTINI paper; the
martini.itp file doesn't appear to list which heavy atoms are represented by
each CG bead either.  For example, I'm looking for something like:

'NC3' : ['C1', 'C2', 'C3', 'N4', 'C5', 'C6'],
'PO4' : ['O7', 'P8', 'O9', 'O10', 'O11','C12'],
'GL1' : ['C13', 'O14', 'C15', 'O16'],
'GL2' : ['C32', 'O33', 'C34', 'O35'],

etc.

For some atoms it's obvious which MARTINI groups they belong in, but others
on the borderline are not obvious.  For example, does C12 belong in PO4 or
GL1?

Anybody have a master list like this?

Thanks,
Mike

-- 

Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] MARTINI / all-atom mapping

[Michael Daily]

Xavier,

I did find the atom2cg.awk script on the downloads-> tools of the martini
website, but there is no corresponding one for lipids.  In any case I can
probably figure out the mapping by trial and error, just based on inter-bead
distances, but it would be nice to have it officially documented.

Mike

On Tue, Aug 30, 2011 at 3:06 AM, XAvier Periole  wrote:

>
> it must be some example of mapping lipids on the website: cgmartini.nl
>
> On Aug 30, 2011, at 3:55 AM, Michael Daily wrote:
>
> Hi all,
>
> I am trying to reverse-map some martini lipids to united atom.  In order to
> do this, I'd prefer to have an EXACT definition of the aa-to-cg mapping.  I
> cannot find this, only an imprecise graphic, in the MARTINI paper; the
> martini.itp file doesn't appear to list which heavy atoms are represented by
> each CG bead either.  For example, I'm looking for something like:
>
> 'NC3' : ['C1', 'C2', 'C3', 'N4', 'C5', 'C6'],
> 'PO4' : ['O7', 'P8', 'O9', 'O10', 'O11','C12'],
> 'GL1' : ['C13', 'O14', 'C15', 'O16'],
> 'GL2' : ['C32', 'O33', 'C34', 'O35'],
>
> etc.
>
> For some atoms it's obvious which MARTINI groups they belong in, but others
> on the borderline are not obvious.  For example, does C12 belong in PO4 or
> GL1?
>
> Anybody have a master list like this?
>
> Thanks,
> Mike
>
> --
> 
> Michael D. Daily
> Postdoctoral research associate
> Pacific Northwest National Lab (PNNL)
> 509-375-4581
> (formerly Qiang Cui group, University of Wisconsin-Madison)
>  --
> 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 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
>



-- 

Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] Questions regarding REMD simulation

[michael zhenin]

Hi All,

I have few questions regarding REMD simulation.

I assume REMD works in the following way: Each replica starts with a given
temperature and the temperature of each replica changes along the simulation
(as described in Sugita and Okamoto publication, 1999). Then demux.pl and
trjcat collect the trajectories of the lowest temperature from each replica.



My problem is that while analyzing the results, I found in the edr
files that the temperatures along each replica are relatively

constant (+- 2K from the starting temperature) and the average temperature
is equal to the initial temperature.


How can this be explained? Does the temperature of each replica change
along the simulation or maybe the temperature of each

simulation stays close to its initial value and the coordinates are changed
between the different replica? (if so, why do i need to use

the post-processing tools such as the demux.pl script?)


Thanks,


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] Membrane protein simulation

[Michael Daily]

Unfortunately genbox will put waters anywhere there is a space, including
inside the membrane.  This can easily be fixed by making a script to remove
waters that are z +/- ~2 nm from the membrane center (you should run
g_density on the system to figure out the optimal distance filter).  You can
run this after running genbox.

On Sun, Sep 18, 2011 at 5:12 PM, Mark Abraham wrote:

>  On 19/09/2011 9:42 AM, Sweta Iyer wrote:
>
> Hi,
> I embedded my protein of interest into a DMPC membrane by the g_membed tool
> with the following command:
>  g membed -f input.tpr -p system.top  -n index.ndx -xyinit 0.1 -xyend 1.0
> -nxy 1000 -zinit 1.1 -zend 1.0 -nz 100
>
>  I then energy minimized the resultant structure for 1 ns before the
> position restraint dynamics and the productive run.
> My structure looks fine after the em. However, when I do the PR and look at
> the structure it looks weird in that half of the protein is hanging out of
> the membrane and there seems to be a patch of water molecules that seem to
> have entered the membrane.
>
>  I have no clue what must be possibly going wrong here. Should I have
> equilibrated  the system for longer than 1ns or is it something wrong with
> the membrane insertion.
>
>
> Sounds like the advice here might be useful.
> http://www.gromacs.org/Documentation/Terminology/Periodic_Boundary_Conditions
>
> 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
>



-- 

Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] opls parameters fir gramicidin terminal groups

[Michael Daily]

Does anyone know where I can find opls-aa parameters for strange peptide
terminal groups, such as the formyl N-terminus and ethanolamine C-terminus
of gramicidin A? I know for charmm you can auto-generate these using
swissparam but I don't know of any equivalent for opls.

Thanks,

-- 

Michael D. Daily
Postdoctoral research associate
Pacific Northwest National Lab (PNNL)
509-375-4581
(formerly Qiang Cui group, University of Wisconsin-Madison)
-- 
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] TI/FEP, BAR, and topologies

[Michael Brunsteiner]

Hi,

I'd like to perform TI calculations as described in section 3.12.2 of the 
(version 4.5.4) manual.
my questions are:

1) i understand in Gromacs TI/FEP is implemented as a single topology, and not 
dual topolgy
algorithm, is that correct?


2)  to successfully analyze the results with BAR what is a good frequency at 
which i would

save energies and dU/dlambda values?
3) how do i set up topologies in practice?

If for example i wanted to mutate an ASN to a PHE residue in solution I seem to 
have at least two options:

  a) a direct mutation, where i'd mutate, e.g., CD1  and CD2 in PHE to OD1 and 
ND2 respectively in ASN,

and so forth, including a transformation of a few PHE atoms to dummies, and a 
transformation
of some bonded energy terms.


  b) construct a bastard residue from a PHE and a ASN residue, where the two 
residues share

the backbone atoms, including C-alpha, but where there are two full and 
different copies of the side chain -
and then mutate from state A with all sidechain atoms of, say, ASN set to 
dummies and the PHE atoms
fully turned on, to state B (vice versa) ... this would be like emulating a of 
dual topology approach,
wouldn't it?


is there any good reason for using either a or b? To me option a seems to be 
the more tedious one
as i'd have to include the changes in bonded energy terms by hand into the 
topology file, while with
option b pdb2gmx and grompp would do most of that work for me, correct? Also if 
i used option a
then at the ASN state would I not get some spurious effects from the ring 
topology of PHE, from the PHE
atoms that have their non-bonded interactions turned off but their bonded 
interactions are still there??

(i seem to recall that totally turning off any bonded interactions opens a can 
of worms)
thanks in advance for any help!

cheers
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] 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] 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] 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] 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] 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

[gmx-users] Charmm Sugars/carbohydrates

[Michael McGovern]

Hi everyone. I'm doing a simulation in gromacs using the charmm36 parameters 
from the gromacs website. The parameters don't seem to have carbohydrates, 
which are part of the charmm36 force field. In particular, I need parameters 
for trehalose. Is there anywhere I can find these parameters in a 
format suitable for gromacs?-- 
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] shell molecular dynamics

[michael zhenin]

Hello,

I would like to run a simulation using the shell option.

I tried to restraint solvent molecules which were found within the range
of 5-10 angstroms from the box center. I couldn't find any way to do
so. I didn't find any information about this topic anywhere else
i.e. GROMACS manual, online tutorials etc.

Can anybody help me with this problem or guide me to appropriate tutorial.

Thanks in advance.

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] 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: 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] 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

[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]

> 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

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] 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] 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] 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] 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]

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] 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: 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