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

[chris . neale]

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.


Please note that I didn;t check everything thoroughly or look at the  
other images so there could still be something weird going on, but I  
disagree that http://www.brunsteiner.net/tbut-pmf.gif shows anything  
strange.


Chris.

-- original message --


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

[gmx-users] Can g_wham support using different temperature for different windows?

[chris . neale]

Sounds like unconverged sampling. You would be astounded how long  
systems like this can take to converge. An all-atom simulation like  
this can easily require >10 us (microseconds!) per umbrella. I don't  
know about martini, probably a lot less.


Chris.

-- original message --

Thanks Justin.
I tried your suggestions by either increase more windows and change the force
constant, but it seems the samplings are still bad in some windows. When I did
pulling in (0 0 1) direction and a reverse pulling in (0 0 -1)  
direction, I got

different configurations at certain reaction coordinates. And the windowed
umbrella sampling seems depends strongly on the initial configurations in that
window. Therefore I got different PMFs using pulling in (0 0 1) direction and
reverse pulling in ?0 0 -1) direction.


In my simulation, I exert constraints on phosphate atoms in z direction, so
there is no lipid flip-flop and the membrane will be stable at high
temperatures. Then I am thinking of increasing temperature in those  
bad windows

to enhance sampling...

best regards,
Jianguo





--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] test on g_density

[chris . neale]

Could it be this:

http://www.mail-archive.com/gmx-users@gromacs.org/msg35107.html

I've posted a few times on the fact that this tool is broken for  
constant pressure simulations.


I'm not sure why both even and odd would be under the overall average,  
but it seems possible.


Chris.

-- original message --


Hi all.

I've noticed an unexpected behaviour using g_density. I have a
trajectory (with time step 2) and it is split into two subsets (using
trjconv) with, lets say, even and odd steps respectively. According
the the algorithm I expected that the density obtained with the
original trajectory would be the average of the density with each
generated trajectory. However, that was not the case, and for each
slice, while density calculated from the subsets was practically the
same, the density calculated with the original was different (greater
or lower). My system is a membrane, the original trajectory is 50ns
long (25001 steps) and I calculate the density along the bilayer
normal (z). I am using GROMACS4.5.3 and a I've issued the following
commands:

*To get the subsets:
trjconv -f bilayer.xtc -s bilayer.trp -o bilayer_odd -dt 4
trjconv -f bilayer.xtc -s bilayer.trp -o bilayer_even -dt 4 -b 2

The generated trajectories have 12501 and 12500 steps respectively.

*To calculate the densities:
g_density -f bilayer.xtc -s bilayer.tpr -o dens_all
g_density -f bilayer_odd.xtc -s bilayer.tpr -o dens_odd
g_density -f bilayer_even.xtc -s bilayer.tpr -o dens_even

The result is attached. dens_odd and dens_even are very similar, but
surprisingly (at least for me), dens_all is different (it actually
integrates more atoms than the other two).

I haven't seen anything in the code that could explain that result
(maybe when removing pbc??). Is it a bug or a mistake from me?

Thank you!

Javier
-- next part --


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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

[chris . neale]

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)". 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.


Chris.

-- original message --

Re: PMF from pull code, unexpected results
Michael Brunsteiner mbx0009 at yahoo.com
Wed Feb 23 23:16:07 CET 2011

* Previous message: [gmx-users] widom insertion
* Next message: [gmx-users] trr files transferability
* Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

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] Questions on combination of Berger's united-atom force field for lipid and OPLS_AA force field for protein

[chris . neale]

Your first question is addressed in the following document:

http://www.pomeslab.com/files/lipidCombinationRules.pdf

Your second and third question are addressed (in general form) in the  
gromacs manual.


If you want a topology file, send me a personal email indicating the  
lipid that you are using (no other questions by personal email though  
please).


Chris.

-- original message --

Dear All:

I am using GROMACS package to simulate membrane proteins. I plan to =20
use Berger's united atom force field (Berger's UA FF) for lipids and =20
OPLS_AA force field (OPLS_AA FF) for proteins.

I have read the guide post kindly by Prof. Dr. Chris Neale in previous =20
mailing lists very carefully. According to his guide, the following =20
steps are needed:

1). Added [atomtypes] from lipid.itp to ffoplsaanb.itp -- after =20
changing c6/c12
to sigma/epsilon. Also added atomtype H from olsa_369 to match H expected by
pope.itp
   - sigma   =3D (c12/c6)^1/6
   - epsilon =3D c6/(4*sigma^6)

2). Added [pairtypes] from lipid.itp to ffoplsaanb.itp -- after =20
changing c6/c12
to sigma/epsilon. (gives effective fudgeLJ of 0.125). Also changed all =20
reference
to OW to opls_116 (opls spc water oxygen) and simply removed any with =20
reference
to HW as it will be zero regardless.

3). Added [dihedraltypes] from lipid.itp to ffoplsaabon.itp.
   - Prior to running ensure that the non-RB dihedral does not exist for the=
se
groups.

4). make a topology file like this:

My questions are as follows:

1, Berger's United atom force field scale LJ1-4 interaction by 0.125 =20
and scale Coulombic1-4 interaction by 1.0. OPLS_AA FF scale both of =20
them by 0.5, right? From the above changes, I find that the changes of =20
the sigma or epsilon (or c6 and c12) is only associate with the LJ1-4 =20
interaction. So, how the Coulombic1-4 interaction is scaled properly? =20
Do I need additional changes?

2, the non-bond interactions between the atoms from the lipid and the =20
atoms from the protein are need to be calculated, right?  My question =20
is how these interactions are calculated? Does the changes in the =20
[atomtypes]affect these interactions?

3, How the non-bond interactions between the atoms from the lipid are =20
calculated (e.g., the LJ and coulombic interactions)? Does the changes =20
in the [atomtypes]affect these interactions?

Thank you very much for your time and your kindness. I really =20
appreciate your help.

Regards

Ruo-Xu Gu


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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 one thermostat for entire system ??

[chris . neale]

I'm not sure what the reviewer has in mind. Therefore I would split  
the response into a number of possibilities.


1. If you simulated a box of those solute chemicals in the absence of  
a large solute, and if those solutes do not phase separate, then the  
collisions will equipartition the kinetic energy. The solutes are all  
small and chemically similar insofar as our forcefields treat them.


2. If one then added a large solute, one may be concerned that the  
thermal fluctuations of the solute will partition separately from the  
solute. This is a real concern and this is why I use Langevin dynamics  
(the sd integrator). But you seem to be not very concerned about the  
thermal fluctuations of the solute, so this, on its own, seems to be  
not a concern.


3. The solvent beside a frozen solute: this seems to be a serious  
possibility of a problem. Did you freeze it entirely? If so, does  
gromacs discount this from the kinetic energy? If so, then all seems  
fine. But if you used harmonic restraints then does that somehow suck  
energy out of nearby solvent molecules? I don't know but I might be  
concerned about that if I was trying to draw kinetic conclusions about  
the system.


4. What are your conclusions? Do they even depend on the equipartition  
of energy? If so, how?


So my base answer is that yes, you do possibly have problems and these  
can all be solved with Langevin dynamics (but then kinetics are  
irrelevant and conclusions can not be drawn). But more applicably,  
might these problems affect your conclusions? and here I honestly  
doubt that you have any problem by coupling solvents together. For  
example, if you coupled a box of water together,

 would you expect problems with energy partitioning? I think not.

Finally, to the reviewer's probable question, you have a frozen or  
restrained solute so likely you are not possibly experiencing a frozen  
ice cube type problem. If the reviewer has a particular problem in  
mind, perhaps they would be kind enough to share that with you so that  
you could address it in specific.


Chris.

-- original message --

Thanks Mark for the advice.

I have just rerun a test simulation with each of my gas species
coupled separately to a thermostat and have got similar values for the
quantities of interest (diffusion coefficients). However I am not sure
that this will satisfy the reviewer without a bit more justification
of why I chose to use a single thermostat for my system.

I remembered reading some gromacs information on the application of
thermostats (on which I based my decision to apply one thermostat for
the entire system). I have just managed to locate it on the FAQ page:

http://www.gromacs.org/Documentation/Terminology/Thermostats

What Not To Do
Some hints on practices that generally not a good idea to use:
?   Do not use separate thermostats for every component of your system.
Some molecular dynamics thermostats only work well in the
thermodynamic limit.  A group must be of sufficient size to justify
its own thermostat.  If you use one thermostat for, say, a small
molecule, another for protein, and another for water, you are likely
introducing errors and artifacts that are hard to predict. In
particular, do not couple ions in aqueous solvent in a separate group
from that solvent.

My system consists of
55 molecules of CO2
38 molecules of N2
3 molecules of O2
One large molecule of MCM-41 (mostly frozen but with mobile surface
groups) consisting of 4284 atoms.

My take on this was that the gas molecules were in small supply
compared to the other part of the system so I should avoid using
separate thermostats. Was I mistaken in making this assumption?

Is there any feeling for what "sufficient size" as referred to above
is? Does anyone know of any papers I could reference in my explanation?

Any comments welcome

Thanks

Jenny






Quoting Mark Abraham :


On 8/03/2011 2:50 AM, Jennifer Williams wrote:

Hi,

I simulated the diffusion of small gases (CO2, N2) in a framework   
structure which was mostly frozen with some mobile surface groups.   
I applied a temperature thermostat to the entire system (i.e I   
didn't couple the gas molecules and framework separately). I have   
now been asked the following by a reviewer:


"Please comment on any artefacts that might arise as a result of   
non-equipartition of energies.  For example, what is the calculated  
 temperature of each of the gas species and mobile species?"


I have tried to explore the temperature of each species separately   
but g_energy will only give me the energy of the system as a whole.  
 Are there any other tools within gromacs to look at the  
temperature  of each species in turn from the md runs I already have?


Does anyone have a suggestion as to how I can address the   
reviewer's comment and proove that using one thermostat for the   
whole system is OK (that is what I am hoping!).  Is there anything   
in particular that I should be looking closely at

[gmx-users] Zero Potential of Mean Force with g_wham

[chris . neale]

g_wham is not the only version of wham. Try using alan grossfield's  
version. Too often, I am afraid, gromacs accessory programs get broken  
in an update (not sure what the general solution is here beyond  
renewed calls for a proper test suite. Perhaps having 20+ programs is  
not ideal for a single software suite where the real focus is only on  
mdrun and grompp?)


Chris.

-- original message --

Hi,

I ran  g_wham 4.5.2 and did get a non-zero PMF curve.
I assume that there is something going on with g_wham on version 4.5.1.

Thank you for your help.

Susana

On Wed, Mar 9, 2011 at 3:00 PM, Mark Abraham anu.edu.au>wrote:



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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_wham PMF profile

[chris . neale]

There is no such thing as over-sampling. There is no problem induced  
by having more sampling in one region except in relation to the fact  
that you now have less sampling in some other region. That is to say,  
if you have a totally converged PMF and then sample the first half of  
the umbrellas for another 100x as long, the PMF will not become  
incorrect.


Please note, though, that there are some possible problems here, but  
they involve one's conclusions and not the sampling or the wham. If,  
for example, you concluded that the error was less or the PMF had more  
features in a region that you had sampled for much longer times, then  
you can see how this conclusion could be an artefact of your method.


Generally, you need to have stronger umbrellas and thus closer  
umbrellas to obtain sufficient overlap in places where your PMF is  
concave down. The more strongly concave down the PMF is, the more you  
require umbrellas with strong force constants in that locale.  
Otherwise, you can introduce massive global PMF errors because the two  
umbrellas on either side of the local maximum fall to either side and  
the overlap becomes very poor.


Chris.

-- original message --

Thank you!

and, is there a problem if I'm over-sampling some regions in the case
where in the other regions I have got a good (not a under-one)
sampling?
This over-sampled region can invalidate my DeltaG value?

best,

Anna


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Pulling

[chris . neale]

1. yes. it is acceptable. It is different, but neither method is de  
facto better.


2. to enhance convergence by limiting the amount of phase space that  
must be sampled. Changing the restraints can change the profile, but  
if you care only about the integrated standard binding free energy  
then it does not change the converged result. See, for example, D. L.  
Mobley, J. D. Chodera, K. A. Dill. "On the use of orientational  
restraints and symmetry number corrections in alchemical free energy  
calculations", ...


Chris.

-- original message --

Dear All
Afew question about Pulling in Umberella Sampling

1-the goal of pulling is making some primary structures (in different
distances) to do umberella sampling for each one of them.
 I can make these states by transporting my ligands along a vector to
prepare these primary structures.Is this correct?Now I can do US for each
one!without any need to doing pulling.

2-Why do we keep fix the relative orientation of Protein-ligand during the
pulling ? I think changing the orientation of ligand during the
pulling(suppose the protein is restrain) can chang our result?
  Because our umbrella sampling maintain this orientation too.am I right?

Thanks in advance
-- next part --


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Umberella sampling

[chris . neale]

1. Depends on how you set them in your .mdp file. It could be either.

2. There is no general method. Use trial and error. Also, your  
question is flawed, K defines X. Unless by "length of one windows" you  
meant the distance between neighbouring centers of restraint  
(umbrellas).


3. you need overlapping distributions. Let's leave it at that. I have  
not seen any general treatment of what K is required and I am rather  
sure that it is impossible to predetermine the necessary force  
constants (K) unless you know the PMF. If you need to know the PMF  
beforehand then this is not a general solution and also probably  
useless since the whole purpose is to get the PMF.


Chris.

-- original message --

Dear All
afew question in umberella sampling tutorial:
1-We do umberella sampling for each of 25 simulation windows,while using a
spring(harmonic potential),Are these springs 1 or 3 dimensional?

2-Suppose the length of one windows is X nm,what is the  approperiate K
(spring constant) for this window?Is there a general way to determine this
value?

3-I think the K must be such that the oscilation amplitude be  a few larger
than X/2, because we need overlapping of density distribution for analysing
with wham method, am I right?

Thanks in advance



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Pulling

[chris . neale]

2-When you limit your sampling phase space,it can make some Errors,
Because you may ignore some important regions(where you don't know  
them and you can't predict there).


I agree with the idea that underlies your point, but it is not a  
problem. You don't need to sample all intervening phase space, just  
converge one well defined part of it. Read that paper again, and some  
others.


... snip ...


Is it possible to do Umbrella Sampling while the drug be free to rotate
along with oscilation in windows?


sure, but you will have more convergence problems.

Chris.

On Tue, Mar 15, 2011 at 9:49 AM,  wrote:


1. yes. it is acceptable. It is different, but neither method is de facto
better.

2. to enhance convergence by limiting the amount of phase space that must
be sampled. Changing the restraints can change the profile, but if you care
only about the integrated standard binding free energy then it does not
change the converged result. See, for example, D. L. Mobley, J. D. Chodera,
K. A. Dill. "On the use of orientational restraints and symmetry number
corrections in alchemical free energy calculations", ...

Chris.

-- original message --

Dear All
Afew question about Pulling in Umberella Sampling

1-the goal of pulling is making some primary structures (in different
distances) to do umberella sampling for each one of them.
 I can make these states by transporting my ligands along a vector to
prepare these primary structures.Is this correct?Now I can do US for each
one!without any need to doing pulling.

2-Why do we keep fix the relative orientation of Protein-ligand during the
pulling ? I think changing the orientation of ligand during the
pulling(suppose the protein is restrain) can chang our result?
 Because our umbrella sampling maintain this orientation too.am I right?

Thanks in advance
-- next part --



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Umberella sampling

[chris . neale]

Mohsen:

To be honest, this just sounds like debating. I have offered you my  
perspective already.


Chris.

-- original message --

Dear chris
Thanks for your reply

1-I am not sure,Since we need the know the variations of free energy along a
specific degree of freedom of system(for example z axis),so  the springs
must be 1 dimensional to allow drug to oscilate only in one dimension(z
axis).
Let me say my question in other words:
Do we need to do sampling just in phase space related to z axis(1
dimensional) or we can do sampling in other places(if it be 3
dimensional),Although we still want the variation of free energy along z
axis?

2-evry windows have a width(for example 0.02 nm),Besides according to Hock's
law knowing the oscilation amplitude is possible when you have the K and the
mass and total energy:1/2 k(x/2)^2=1/2 m v(max)^2

3-my question in other words:
Since we need overlapping the distribution densities of each windows,Does
the drug need to enter to the neighbour windows a little ?
Thanks in advance for your reply
On Tue, Mar 15, 2011 at 9:54 AM,  wrote:


1. Depends on how you set them in your .mdp file. It could be either.

2. There is no general method. Use trial and error. Also, your question is
flawed, K defines X. Unless by "length of one windows" you meant the
distance between neighbouring centers of restraint (umbrellas).

3. you need overlapping distributions. Let's leave it at that. I have not
seen any general treatment of what K is required and I am rather sure that
it is impossible to predetermine the necessary force constants (K) unless
you know the PMF. If you need to know the PMF beforehand then this is not a
general solution and also probably useless since the whole purpose is to get
the PMF.

Chris.

-- original message --

Dear All



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

[gmx-users] pick in pull.xvg

[chris . neale]

I did a pulling simulation for protein-drug system.


This could mean almost anything. Please provide details.

Chris.

-- original message --

Dear All

I did a pulling simulation for protein-drug system.

My output (pull.xvg) has a pick in 100 ps (equivalent to 1 nm),
the rest of plot is flat(of course with oscilation around the average),of
course with NOT ZERO average in the flat region.

Can I result ?
" the interaction between protein and drug has finished after 1 nm"

If yes, why the average of flat region is not zero(after the pick)?

Thanks in advance
-- next part --
An HTML attachment was scrubbed...
URL:  
http://lists.gromacs.org/pipermail/gmx-users/attachments/20110319/6447ad26/attachment.html



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] plot of force via time

[chris . neale]

It depends on what you are trying to do.

Also, I'd suggest that you'll get more assistance if you show us that  
you have done some reading and some work to try to do something prior  
to mailing the list.


Chris.

-- original message --

plot of force via time
mohsen ramezanpour ramezanpour.mohsen at gmail.com
Sat Mar 19 14:49:40 CET 2011

* Previous message: [gmx-users] g_sham
* Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

Dear All

I have a plot,Force via time.(pull.xvg in pulling)
what is important,the value of force in each time or its average? because it
is oscilationg around a line.

I am in doubt , because the important in tempreture and pressure coupling
were thier averages not their values in time.

Thanks ina dvance


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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_sham

[chris . neale]

I didn't look, but I would be amazed if this information is not  
availablein the manual or in g_sham -h


It's looking a lot like you are stabbing around in the dark and hoping  
that the list will fill in all the details for you. In fact, often  
people will. But it is better for you if we suggest to you to try  
figuring this out on your own.


Or perhaps I misunderstand your question? f so, then please rephrase  
to make it a lot more specific and to show that you tried to find the  
answer but could not.


Chris.

-- original message --

Dear All

I have a trajectory file for 500 ps simulation
I used g_sham tool,the output was an ener.xvg file (0-40 on y axis and 0-33
on x axis ??!!)

But I don't know what quantities are ploted via eachother?
Please let me know
Thanks in advance


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] compressing a box of water droplets into a homogeneous solution of liquid water

[chris . neale]

Dear users:

I recently came across a system that was composed of tip4p water vapor  
droplets separated by vacuum. This system is what you might get if you  
did a NVT simulation of water with a box that was 10 times too large  
for the number of water molecules.


I was surprised to see that this system did not collapse to any  
significant extent during 200 ps of NPT equilibration at 1 atm using  
the Berendsen thermostat with tau_p=1.0 and the sd integrator and a  
colombic cut-off. (We also tried a number of other integrator/pressure  
coupling combinations with the same results).


I had assumed that such collapse would occur quite rapidly. This does  
not seem to be the case (no noticeable contraction within 200 ps).


Has anybody else done anything like this? Can anybody comment on their  
expectations/experience of collapse from the gas state to the liquid  
state under standard NPT conditions?


Thank you,
Chris.

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

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

[gmx-users] Simulation of slow folding proteins

[chris . neale]

Dear Bharat:

I hope that this doesn't impede others giving you advice about how to  
go about doing what you want to do, but here's my two cents for what  
it's worth:


My suggestion is to forget about trying to do that. 3-ns simulations  
are not going to give you equilibrium populations of conformational  
basins (and neither would 3-us simulations, I suspect). Not every  
question is amenable to MD on the currently available timescales.


If you really want to try to do it, I would find out what  
conformations lead to flourescence and see if you get more drift away  
from those conformations in some models with longer loops. But again,  
I think it's probably going to be a waste of time.


Chris.

-- original message --

Hi,

I simulated a protein (GFP) with one of its loop replaced with a longer loop
. After simulating for some 5 to 10 models with different loop sequences , I
decided to carry out wet-lab experiments for one model on the basis of
simulation result. The analysis of simulation was done in the following way
:-

1) Comparison of RMSD values of backbone - for both GFP wild type and loop
replaced (mutated GFP)
2) Comparison of RMSF of the residues common to both the proteins.
3) Visual Inspection of entry of water into GFP.

The simulation was carried out for 3ns for both the proteins.

After that I did the wet-lab experiment for the modeled structure and I
found the fluorescence to be 50% of the wild type.

Now I want to investigate the reason for this reduction in fluorescence
??... How can this be quantified using MD.


--
Bharat
Ph.D. Candidate
Room No. : 7202A, 2nd Floor
Biomolecular Engineering Laboratory
Division of Chemical Engineering and Polymer Science
Pusan National University
Busan -609735
South Korea
Lab phone no. - +82-51-510-3680, +82-51-583-8343
Mobile no. - 010-5818-3680
E-mail : monu46010 at yahoo.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] compressing a box of water droplets into a homogeneous solution of liquid water

[chris . neale]

Thanks Patrick and Andre!

We repeated this with a few box sizes just to get a quick handle on  
it. The equilibrium volume is about 64 nm^3. If we start with a volume  
of 1000 nm^3 then the overall box does not collapse at all within 200  
ps of NPT Langevin dynamics at 1 atm.If we start with a volume of 200  
nm^3, then it does collapse to approximately 64 nm^3 within 200 ps of  
such simulation.


My best guess is that the rapid collapse is driven by nonbonded  
interactions and thus the rapid collapse does not occur when the  
system is so large with such low density that water forms isolated  
vapour droplets that do not interact with each other by LJ  
interactions. Sure, it is expected to collapse eventually from the 1  
atm pressure coupling, and we have also observed that high pressure  
works, but at 1 atm it might take a very long time to reach equilibrium.


I agree with Andre that none of this matters to regular simulations as  
there is no good reason to go through this type of state when one  
wants to simulate dense liquids. I just found it curious that  
Berendsen pressure coupling at 1 atm was not sufficient to quickly  
equilibrate the volume in a system where the vacuum regions are large  
in comparison to the LJ cutoffs.


Chris.

-- original message --

Hi Chris,
I experienced the same kind of thing. In the process of building a box
of liquid (organic solvent), at some point I wanted to get rid of a
layer of vacuum around my system. So for shrinking the box I used
similar settings as you and found also that the collapse was going
slower than I'd have expected.
One solution to accelerate this (if your goal is to shrink the box) is
to increase the pressure (to say 100 atm). But it's important to stop
the simulation in time (i.e. once the layer of vacuum has disapeared)
otherwise the system shrinks too much and density is off.
So to come back to your system which has a very big layer of vacuum
around, and according to my experience, the volume is probably
decreasing but too slowly to see anything signigicant (compared to the
initial value) in 200 ps .
Ciao,

Patrick

Le 21/03/2011 16:53, chris.ne...@utoronto.ca a écrit :

[Hide Quoted Text]
Dear users:

I recently came across a system that was composed of tip4p water vapor
droplets separated by vacuum. This system is what you might get if you
did a NVT simulation of water with a box that was 10 times too large for
the number of water molecules.

I was surprised to see that this system did not collapse to any
significant extent during 200 ps of NPT equilibration at 1 atm using the
Berendsen thermostat with tau_p=1.0 and the sd integrator and a colombic
cut-off. (We also tried a number of other integrator/pressure coupling
combinations with the same results).

I had assumed that such collapse would occur quite rapidly. This does
not seem to be the case (no noticeable contraction within 200 ps).

Has anybody else done anything like this? Can anybody comment on their
expectations/experience of collapse from the gas state to the liquid
state under standard NPT conditions?

Thank you,
Chris.

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

[gmx-users] Umbrella sampling

[chris . neale]

can you show us your mdp and the pmf and the histograms for the data  
that you put into wham? It's a lot easier to diagnose with the real  
data.


In the general case where umbrellas are spaced equally along your  
reaction coordinate, sampling overlap between umbrellas will always  
decrease anywhere the PMF is concave down, such as your barrier  
region. I would suggest adding windows at every 0.005 nm spanning the  
region that you consider to have a sampling problem (e.g. 0.76 to 0.83  
?) and use a much stronger force constant here (e.g. 2-3 times as  
strong as you used for umbrellas with 0.02 nm spacing). You have  
identified a problem, so my suggestion is to not bother fiddling with  
adding one extra replica but sample that region to death with strong  
force constants. I presume that this will not impact your overall  
efficiency too much if the reaction coordinate is long overall.


Finally, you will need to evaluate the convergence of your PMF overall  
and perhaps of this region in particular, especially if you want to  
know the dG to cross it or between two low energy states on either  
side of it.


Chris.

-- original message --


Sorry I am not sure that I follow. Will the window with r0 =0.80 giving
the distribution centred around 0.78nm not drive my free energy profile
up. If I remove this window prior to running g_wham the free energy goes
down. Should I increase the force constant so that the mean of the
window is 0.80nm (bearing in mind that this is near the barrier region).

Gavin
XAvier Periole wrote:


You can present the data differently:
you have two windows at 0.78 nm giving different distribution.

That indicates these windows are not converged. Does not mean
that the others (0.80 nm) are converged :))

On Mar 31, 2011, at 12:20 PM, Gavin Melaugh wrote:


Hi Xavier

Thanks for the reply. With respect to your answer of my first query.
What if you had two windows practically on top of each other, but one
was not supposed to be there. e.g A window with r0 of 0.80 nm and
centred at 0.78 nm and a window with r0 of 0.78 nm centred at 0.78nm.

Gavin

XAvier Periole wrote:


On Mar 31, 2011, at 11:53 AM, Gavin Melaugh wrote:


Hi All

I have generated several PMF curves for the one system using umbrella
sampling. In the first part of the curve (barrier region) I use a high
force constant with small intervals between the windows. The latter
part
of the curve I use a lower force constant with larger window spacing.
Anyway I have a few issues that I need clarifying:
1 - Can you have too much overlap between windows?

no, there no such a thing of too much overlap :)) You could even put
two identical windows with same 100% overlap ... no problem.

2 - Does the distribution at each window have to centered around the
desired r0? (If not does this affect the free energy?)

The deviation of the distribution from the r0 is what dictates the
profile. The more away from the disired r0 the higher the free energy
of the system.

3- If you over sample one particular window, will it affect the curve?

There is no such a thing of over sampling ... the only thing you can
have is not enough sampling.


Many thanks

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




--
gmx-users mailing listgmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at 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] Umbrella Sampling

[chris . neale]

looks fine to me, no need to do that extra sampling that I suggested  
since it appears that you already did this -- benefits of seeing real  
data ;). If you want to understand why your histograms are not always  
centered at r0 (note that this is just fine) then you should read more  
about US, WHAM, and how to bias/debias the data for US (I am sure that  
there are textbooks around that explain this). The only case in which  
all of your histograms will be centered at their respective r0 is when  
the underlying PMF is exactly flat.


Chris.

-- original message --

Hi Chris many thanks again for the advise. I have, or at least I thought
have sampled my barrier region to death, but as I say some histograms
may not be centred around r0. I will proceed with what you suggest.
Please find attached a picture of the histograms, the corresponding
profile, and a sample mdp file that I use.

-- next part --
A non-text attachment was scrubbed...
Name: groprofile.agr.gz
Type: application/x-gzip
Size: 3805 bytes
Desc: not available
Url :  
http://lists.gromacs.org/pipermail/gmx-users/attachments/20110331/2a5b0bd6/groprofile.agr.gz

-- next part --
A non-text attachment was scrubbed...
Name: hist.agr.gz
Type: application/x-gzip
Size: 25250 bytes
Desc: not available
Url :  
http://lists.gromacs.org/pipermail/gmx-users/attachments/20110331/2a5b0bd6/hist.agr.gz


--

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Umbrella Sampling

[chris . neale]

your comment:

which should be centred around 0.80nm

is flawed, as i mentioned earlier. also, it is not g_wham that is  
sensitive but the convergence and sampling of phase space that is  
sensitive. don`t remove any data. do evaluate your convergence.  
without convergence measures, a pmf is worse than useless.


chris.

-- original message --

Cheers Chris

If I remove the red histogram (the first of the wider distributions),
which should be centred around 0.80nm but is actually centred around
0.78 nm; and add in some more histograms with higher force constants the
profile changes slightly. It seems that  g_wham is very sensitive to
these subtleties. How do I know which curve is correct? I have about six
such curves that differ slightly in this manner.

Gavin

chris.ne...@utoronto.ca wrote:

[Hide Quoted Text]
looks fine to me, no need to do that extra sampling that I suggested
since it appears that you already did this -- benefits of seeing real
data ;). If you want to understand why your histograms are not always
centered at r0 (note that this is just fine) then you should read more
about US, WHAM, and how to bias/debias the data for US (I am sure that
there are textbooks around that explain this). The only case in which
all of your histograms will be centered at their respective r0 is when
the underlying PMF is exactly flat.

Chris.

-- original message --

Hi Chris many thanks again for the advise. I have, or at least I thought
have sampled my barrier region to death, but as I say some histograms
may not be centred around r0. I will proceed with what you suggest.
Please find attached a picture of the histograms, the corresponding
profile, and a sample mdp file that I use.


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Umbrella Sampling

[chris . neale]

1. how do you know that data is "not good"? It is not good science to  
remove data just because it doesn't match with your expected results.  
For now, there is no compelling reason to remove that data.


2. My first step in checking the convergence is to block average the  
time. If you have 10 ns total per umbrella, then extract the first 0-1  
ns  and make a PMF based on that. Then extract the 1-2 ns data and  
make a PMF based on that. Continue until you make a PMF based on the  
9-10 ns data. Now check and see if this PMF is drifting over time. You  
can do this by looking at the PMFs or by integrating the bound state  
and plotting the binding free energy as a function of equilibration  
time. If the binding free energy is still drifiting as a function of  
increasing equilibration time then you are certainly unconverged. Note  
that there will always be noise in this measure, but what you are  
looking for is unidirectional drift in the value. -- One could write a  
textbook on this topic. i suggest that you seek out alternative  
resources to understand how to check convergence properly.


Chris.


Cheers Chris

Is the best way to check for convergence; to keep adding in more
histograms until the curves converge. Also your comment  'don't remove
any data', do you mean to keep histograms that are not so good.

Gavin
chris.neale at utoronto.ca wrote:

your comment:

which should be centred around 0.80nm

is flawed, as i mentioned earlier. also, it is not g_wham that is
sensitive but the convergence and sampling of phase space that is
sensitive. don`t remove any data. do evaluate your convergence.
without convergence measures, a pmf is worse than useless.

chris.

-- original message --

Cheers Chris

If I remove the red histogram (the first of the wider distributions),
which should be centred around 0.80nm but is actually centred around
0.78 nm; and add in some more histograms with higher force constants the
profile changes slightly. It seems that  g_wham is very sensitive to
these subtleties. How do I know which curve is correct? I have about six
such curves that differ slightly in this manner.

Gavin

chris.neale at utoronto.ca wrote:

[Hide Quoted Text]
looks fine to me, no need to do that extra sampling that I suggested
since it appears that you already did this -- benefits of seeing real
data ;). If you want to understand why your histograms are not always
centered at r0 (note that this is just fine) then you should read more
about US, WHAM, and how to bias/debias the data for US (I am sure that
there are textbooks around that explain this). The only case in which
all of your histograms will be centered at their respective r0 is when
the underlying PMF is exactly flat.

Chris.

-- original message --

Hi Chris many thanks again for the advise. I have, or at least I thought
have sampled my barrier region to death, but as I say some histograms
may not be centred around r0. I will proceed with what you suggest.
Please find attached a picture of the histograms, the corresponding
profile, and a sample mdp file that I use.



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] FEP and loss of performance

[chris . neale]

If we accept your text at face value, then the simulation slowed down  
by a factor of 1500%, certainly not the 16% of the load balancing.


Please let us know what version of gromacs and cut and paste your  
cammands that you used to run gromacs (so we can verify that you ran  
on the same number of processors) and cut and paste a diff of the .mdp  
files (so that we can verify that you ran for the same number of steps).


You might be correct about the slowdown, but let's rule out some other  
more obvious problems first.


Chris.

-- original message --


Dear all,
when I run a single free energy simulation
i noticed that there is a loss of performace with respect to
the normal MD

free_energy= yes
init_lambda= 0.9
delta_lambda   = 0.0
couple-moltype = Protein_Chain_P
couple-lambda0 = vdw-q
couple-lambda0 = none
couple-intramol= yes

   Average load imbalance: 16.3 %
   Part of the total run time spent waiting due to load imbalance: 12.2 %
   Steps where the load balancing was limited by -rdd, -rcon and/or -dds: X0 %
   Time:   1852.712   1852.712100.0

free_energy= no
   Average load imbalance: 2.7 %
   Part of the total run time spent waiting due to load imbalance: 1.7 %
   Time:127.394127.394100.0

It seems that the loss of performace is due in part to in the load imbalance
in the domain decomposition, however I tried to change
these keywords without benefit
Any comment is welcome.

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

[gmx-users] minimization and simulation problems

[chris . neale]

What are your initial box dimensions prior to em? Also, please copy  
and paste your .mdp options. Also, what happens when you run the same  
post-em simulation with nsteps=1 ?


-- original message --


Dear all,
I'm trying to run simulation of 30 proteins in water using the Martini
force field. I used water.gro file in order to solvate the proteins.
For minimization I used the em.mdp file published at Martini site
(http://md.chem.rug.nl/cgmartini/index.php/home). When I set the emtol
parameter to 10 the system can't converge. So I used emtol 100 and
then the system converged. I use it as an input for the simulation.
The file can't be attached as it is too big nut I can send it if needed.
However, the sumulation crushes when I'm trying to run MD using md.mdp
also from the Martini site. I'm getting the following warnings and
errors:
Warning: Only triclinic boxes with the first vector parallel to the
x-axis and the second vector in the xy-plane are supported.
  Box (3x3):
 Box[0]={ nan,  nan,  nan}
 Box[1]={ nan,  nan,  nan}
 Box[2]={ nan,  nan,  nan}
  Can not fix pbc.
Warning: Only triclinic boxes with the first vector parallel to the
x-axis and the second vector in the xy-plane are supported.
  Box (3x3):
 Box[0]={ nan,  nan,  nan}
 Box[1]={ nan,  nan,  nan}
 Box[2]={ nan,  nan,  nan}
  Can not fix pbc.
Warning: Only triclinic boxes with the first vector parallel to the
x-axis and the second vector in the xy-plane are supported.
  Box (3x3):
 Box[0]={ nan,  nan,  nan}
 Box[1]={ nan,  nan,  nan}
 Box[2]={ nan,  nan,  nan}
  Can not fix pbc.
Warning: Only triclinic boxes with the first vector parallel to the
x-axis and the second vector in the xy-plane are supported.
  Box (3x3):
 Box[0]={ nan,  nan,  nan}
 Box[1]={ nan,  nan,  nan}
 Box[2]={ nan,  nan,  nan}
  Can not fix pbc.
Warning: Only triclinic boxes with the first vector parallel to the
x-axis and the second vector in the xy-plane are supported.
  Box (3x3):
 Box[0]={ nan,  nan,  nan}
 Box[1]={ nan,  nan,  nan}
 Box[2]={ nan,  nan,  nan}
  Can not fix pbc.
Warning: Only triclinic boxes with the first vector parallel to the
x-axis and the second vector in the xy-plane are supported.
  Box (3x3):
 Box[0]={ nan,  nan,  nan}
 Box[1]={ nan,  nan,  nan}
 Box[2]={ nan,  nan,  nan}
  Can not fix pbc.

---
Program mdrun_mpi, VERSION 4.0.3
Source code file: nsgrid.c, line: 348

Fatal error:
Number of grid cells is zero. Probably the system and box collapsed.

---

"It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

Error on node 0, will try to stop all the nodes
Halting parallel program mdrun_mpi on CPU 0 out of 8

gcq#166: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0[cli_0]:
aborting job:
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0

---
Program mdrun_mpi, VERSION 4.0.3
Source code file: nsgrid.c, line: 348

Fatal error:
Number of grid cells is zero. Probably the system and box collapsed.

---
Error on node 1, will try to stop all the nodes
Halting parallel program mdrun_mpi on CPU 1 out of 8

gcq#166: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

application called MPI_Abort(MPI_COMM_WORLD, -1) - process 3[cli_3]:
aborting job:
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 3

gcq#166: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

application called MPI_Abort(MPI_COMM_WORLD, -1) - process 5[cli_5]:
aborting job:
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 5

gcq#166: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

application called MPI_Abort(MPI_COMM_WORLD, -1) - process 4[cli_4]:
aborting job:
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 4

gcq#166: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

application called MPI_Abort(MPI_COMM_WORLD, -1) - process 6[cli_6]:
aborting job:
application called MPI_Abort(MPI_COMM_WORLD, -1)

[gmx-users] FEP and loss of performance

[Chris Neale]

Load balancing problems I can understand, but why would it take longer 
in absolute time? I would have thought that some nodes would simple be 
sitting idle, but this should not cause an increase in the overall 
simulation time (15x at that!).


There must be some extra communication?

I agree with Justin that this seems like a strange thing to do, but 
still I think that there must be some underlying coding issue (probably 
one that only exists because of a reasonable assumption that nobody 
would annihilate the largest part of their system).


Chris.


Luca Bellucci wrote:

/  Hi Chris,

/>/  thank for the suggestions,
/>/  in the previous mail there is a mistake because
/>/  couple-moltype = SOL (for solvent) and not "Protein_chaim_P".
/>/  Now the problem of the load balance seems reasonable, because
/>/  the water box is large ~9.0 nm.
/
Now your outcome makes a lot more sense.  You're decoupling all of the solvent?
  I don't see how that is going to be physically stable or terribly meaningful,
but it explains your performance loss.  You're annihilating a significant number
of interactions (probably the vast majority of all the nonbonded interactions in
the system), which I would expect would cause continuous load balancing issues.

-Justin


/  However the problem exist and the performance loss is very high, so I have

/>/  redone calculations with this command:
/>/
/>/  grompp -f
/>/  md.mdp -c ../Run-02/confout.gro -t ../Run-02/state.cpt -p ../topo.top -n 
../index.ndx -o
/>/  md.tpr -maxwarn 1
/>/
/>/  mdrun -s md.tpr -o md
/>/
/>/  this is part of the md.mdp file:
/>/
/>/  ; Run parameters
/>/  ; define  = -DPOSRES
/>/  integrator  = md;
/>/  nsteps  = 1000  ;
/>/  dt  = 0.002 ;
/>/  [..]
/>/  free_energy= yes ; /no
/>/  init_lambda= 0.9
/>/  delta_lambda   = 0.0
/>/  couple-moltype = SOL; solvent water
/>/  couple-lambda0 = vdw-q
/>/  couple-lambda1 = none
/>/  couple-intramol= yes
/>/
/>/  Result for free energy calculation
/>/   Computing: Nodes Number G-CyclesSeconds %
/>/  ---
/>/   Domain decomp.   8126   22.0508.3 0.1
/>/   DD comm. load  8 150.0090.0 0.0
/>/   DD comm. bounds 8 120.0310.0 0.0
/>/   Comm. coord.8   1001   17.3196.5 0.0
/>/   Neighbor search8127  436.569  163.7 1.1
/>/   Force   8   100134241.57612840.9
87.8
/>/   Wait + Comm. F8   1001   19.4867.3 0.0
/>/   PME mesh  8   1001 4190.758 1571.610.7
/>/   Write traj.  8  71.8270.7 0.0
/>/   Update  8   1001   12.5574.7 0.0
/>/   Constraints   8   1001   26.4969.9 0.1
/>/   Comm. energies  8   1002   10.7104.0 0.0
/>/   Rest   8  25.1429.4 0.1
/>/  ---
/>/   Total  8   39004.53114627.1   100.0
/>/  ---
/>/  ---
/>/   PME redist. X/F  8   3003 3479.771 1304.9 8.9
/>/   PME spread/gather   8   4004  277.574  104.1 0.7
/>/   PME 3D-FFT   8   4004  378.090  141.8 1.0
/>/   PME solve  8   2002   55.033   20.6 0.1
/>/  ---
/>/  Parallel run - timing based on wallclock.
/>/
/>/ NODE (s)   Real (s)  (%)
/>/ Time:   1828.385   1828.385100.0
/>/ 30:28
/>/   (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
/>/  Performance:  3.115  3.223  0.095253.689
/>/
/>/   I Switched off only the free_energy keyword and I redone the calculation
/>/  I have:
/>/   Computing: Nodes Number G-CyclesSeconds %
/>/  ---
/>/   Domain decomp.  8 77   10.9754.1 0.6
/>/   DD comm. load 8  10.0010.0 0.0
/>/   Comm. coord.   8   1001   14.4805.4 0.8
/>/   Neighbor search   8 78  136.479   51.2 7.3
/>/   Force 8   1001 1141.115  427.961.3
/>/   Wait + Comm. F  8   1001   17.8456.7 1.0
/>/   PME mesh8   1001  484.581  181.726.0
/>/   Write traj.   8  51.22

[gmx-users] minimization and simulation problems

[Chris Neale]

Can you please redo the md part with gen_vel=yes and see if that makes 
any difference?


Generally, you need to narrow down the problem for us. Does it crash in 
serial as well as parallel? How many steps does it go before the crash? 
what happens to the system volume as a function of time for the duration 
of the simulation prior to the crash.


Chris.

Quoting politr at fh.huji.ac.il 
:


Dear gromacs users,

/  my box dimensions are 368A and when I run the simulation with

/>/  nsteps=1 it works fine. The mdp files used for minimization and
/>/  post-em simulation are attached.
/Thanks again for your help.
Regina

/

/>/
/>/  Quotingchris.neale at utoronto.ca  
:
/>/
/>>/  What are your initial box dimensions prior to em? Also, please copy
/>>/   and paste your .mdp options. Also, what happens when you run the
/>>/  same post-em simulation with nsteps=1 ?
/>>/
/>>/  -- original message --
/>>/
/>>/
/>>/  Dear all,
/>>/  I'm trying to run simulation of 30 proteins in water using the Martini
/>>/  force field. I used water.gro file in order to solvate the proteins.
/>>/  For minimization I used the em.mdp file published at Martini site
/>>/  (http://md.chem.rug.nl/cgmartini/index.php/home). When I set the emtol
/>>/  parameter to 10 the system can't converge. So I used emtol 100 and
/>>/  then the system converged. I use it as an input for the simulation.
/>>/  The file can't be attached as it is too big nut I can send it if needed.
/>>/  However, the sumulation crushes when I'm trying to run MD using md.mdp
/>>/  also from the Martini site. I'm getting the following warnings and
/>>/  errors:
/>>/  Warning: Only triclinic boxes with the first vector parallel to the
/>>/  x-axis and the second vector in the xy-plane are supported.
/>>/   Box (3x3):
/>>/  Box[0]={ nan,  nan,  nan}
/>>/  Box[1]={ nan,  nan,  nan}
/>>/  Box[2]={ nan,  nan,  nan}
/>>/   Can not fix pbc.
/>>/  Warning: Only triclinic boxes with the first vector parallel to the
/>>/  x-axis and the second vector in the xy-plane are supported.
/>>/   Box (3x3):
/>>/  Box[0]={ nan,  nan,  nan}
/>>/  Box[1]={ nan,  nan,  nan}
/>>/  Box[2]={ nan,  nan,  nan}
/>>/   Can not fix pbc.
/>>/  Warning: Only triclinic boxes with the first vector parallel to the
/>>/  x-axis and the second vector in the xy-plane are supported.
/>>/   Box (3x3):
/>>/  Box[0]={ nan,  nan,  nan}
/>>/  Box[1]={ nan,  nan,  nan}
/>>/  Box[2]={ nan,  nan,  nan}
/>>/   Can not fix pbc.
/>>/  Warning: Only triclinic boxes with the first vector parallel to the
/>>/  x-axis and the second vector in the xy-plane are supported.
/>>/   Box (3x3):
/>>/  Box[0]={ nan,  nan,  nan}
/>>/  Box[1]={ nan,  nan,  nan}
/>>/  Box[2]={ nan,  nan,  nan}
/>>/   Can not fix pbc.
/>>/  Warning: Only triclinic boxes with the first vector parallel to the
/>>/  x-axis and the second vector in the xy-plane are supported.
/>>/   Box (3x3):
/>>/  Box[0]={ nan,  nan,  nan}
/>>/  Box[1]={ nan,  nan,  nan}
/>>/  Box[2]={ nan,  nan,  nan}
/>>/   Can not fix pbc.
/>>/  Warning: Only triclinic boxes with the first vector parallel to the
/>>/  x-axis and the second vector in the xy-plane are supported.
/>>/   Box (3x3):
/>>/  Box[0]={ nan,  nan,  nan}
/>>/  Box[1]={ nan,  nan,  nan}
/>>/  Box[2]={ nan,  nan,  nan}
/>>/   Can not fix pbc.
/>>/
/>>/  ---
/>>/  Program mdrun_mpi, VERSION 4.0.3
/>>/  Source code file: nsgrid.c, line: 348
/>>/
/>>/  Fatal error:
/>>/  Number of grid cells is zero. Probably the system and box collapsed.
/>>/
/>>/  ---
/>>/
/>>/  "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)
/>>/
/>>/  Error on node 0, will try to stop all the nodes
/>>/  Halting parallel program mdrun_mpi on CPU 0 out of 8
/>>/
/>>/  gcq#166: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)
/>>/
/>>/  application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0[cli_0]:
/>>/  aborting job:
/>>/  application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0
/>>/
/>>/  -

[gmx-users] FEP and loss of performance

[Chris Neale]

>> Dear Chris and Justin


/  Thank you for your precious suggestions

/>>/  This is a test that i perform in a single machine with 8 cores
/>>/  and gromacs 4.5.4.
/>>/
/>>/  I am trying  to enhance the  sampling of a protein using the decoupling 
scheme
/>>/  of the free energy module of gromacs.  However when i decouple only the
/>>/  protein, the protein collapsed. Because i simulated in NVT i thought that
/>>/  this was an effect of the solvent. I was trying to decouple also the 
solvent
/>>/  to understand the system behavior.
/>>/
/>

Rather than suspect that the solvent is the problem, it's more likely that
decoupling an entire protein simply isn't stable.  I have never tried anything
that enormous, but the volume change in the system could be unstable, along with
any number of factors, depending on how you approach it.

If you're looking for better sampling, REMD is a much more robust approach than
trying to manipulate the interactions of huge parts of your system using the
free energy code.


Presumably Luca is interested in some type of hamiltonian exchange where lambda 
represents the interactions between the protein and the solvent?
This can actually be a useful method for enhancing sampling. I think it's dangerous if we 
rely to heavily on "try something else". I still see no methodological
reason a priori why there should be any actual slowdown, so that makes me think 
that it's an implementation thing, and there is at least the possibility that 
this is
something that could be fixed as an enhancement.

Chris.


-Justin


/   I expected a loss of performance, but not so drastic.

/>/  Luca
/>/
/>>/  Load balancing problems I can understand, but why would it take longer
/>>/  in absolute time? I would have thought that some nodes would simple be
/>>/  sitting idle, but this should not cause an increase in the overall
/>>/  simulation time (15x at that!).
/>>/
/>>/  There must be some extra communication?
/>>/
/>>/  I agree with Justin that this seems like a strange thing to do, but
/>>/  still I think that there must be some underlying coding issue (probably
/>>/  one that only exists because of a reasonable assumption that nobody
/>>/  would annihilate the largest part of their system).
/>>/
/>>/  Chris.
/>>/
/>>/  Luca Bellucci wrote:
/>>>/  /  Hi Chris,
/>>/  />/  thank for the suggestions,
/>>/  />/  in the previous mail there is a mistake because
/>>/  />/  couple-moltype = SOL (for solvent) and not "Protein_chaim_P".
/>>/  />/  Now the problem of the load balance seems reasonable, because
/>>/  />/  the water box is large ~9.0 nm.
/>>/  /
/>>/  Now your outcome makes a lot more sense.  You're decoupling all of the
/>>/  solvent? I don't see how that is going to be physically stable or terribly
/

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] FEP and loss of performance

[chris . neale]

I don't know if it is possible or not. I think that you can enhance  
your chances of developer attention if you develop a small and simple  
test system that reproduces the slowdown and very explicitly state  
your case for why you can't use some other method. I would suggest  
posting that to the mailing list and, if you don't get any response,  
post it as an enhancement request on the redmine page (or whatever has  
taken over from bugzilla).


Good luck,
Chris.

-- original message --


Yes i am testing the possibility to perform an Hamiltonian-REMD
Energy barriers can be overcome  increasing the temperature system or scaling
potential energy  with a lambda value, these methods are "equivalent".
Both have advantages and disavantages, at this stage it is not the right place
to debate on it. The main problem seems to be how to overcome to the the loss
of gromacs performance in such calculation.  At this moment it seems an
intrinsic code problem.
Is it possible?


 >> Dear Chris and Justin
>>
>>/  Thank you for your precious suggestions

/>>/  This is a test that i perform in a single machine with 8 cores
/>>/  and gromacs 4.5.4.
/>>/
/>>/  I am trying  to enhance the  sampling of a protein using the
decoupling scheme />>/  of the free energy module of gromacs.  However when
i decouple only the />>/  protein, the protein collapsed. Because i
simulated in NVT i thought that />>/  this was an effect of the solvent. I
was trying to decouple also the solvent />>/  to understand the system
behavior.
/>>/
/>

>Rather than suspect that the solvent is the problem, it's more likely that
>decoupling an entire protein simply isn't stable.  I have never tried
> anything that enormous, but the volume change in the system could be
> unstable, along with any number of factors, depending on how you approach
> it.
>
>If you're looking for better sampling, REMD is a much more robust approach
> than trying to manipulate the interactions of huge parts of your system
> using the free energy code.

Presumably Luca is interested in some type of hamiltonian exchange where
lambda represents the interactions between the protein and the solvent?
This can actually be a useful method for enhancing sampling. I think it's
dangerous if we rely to heavily on "try something else". I still see no
methodological reason a priori why there should be any actual slowdown, so
that makes me think that it's an implementation thing, and there is at
least the possibility that this is something that could be fixed as an
enhancement.

Chris.


-Justin

>/   I expected a loss of performance, but not so drastic.

/>/  Luca
/>/
/>>/  Load balancing problems I can understand, but why would it take
longer />>/  in absolute time? I would have thought that some nodes would
simple be />>/  sitting idle, but this should not cause an increase in the
overall />>/  simulation time (15x at that!).
/>>/
/>>/  There must be some extra communication?
/>>/
/>>/  I agree with Justin that this seems like a strange thing to do, but
/>>/  still I think that there must be some underlying coding issue
(probably />>/  one that only exists because of a reasonable assumption
that nobody />>/  would annihilate the largest part of their system).
/>>/
/>>/  Chris.
/>>/
/>>/  Luca Bellucci wrote:
/>>>/  /  Hi Chris,
/>>/  />/  thank for the suggestions,
/>>/  />/  in the previous mail there is a mistake because
/>>/  />/  couple-moltype = SOL (for solvent) and not "Protein_chaim_P".
/>>/  />/  Now the problem of the load balance seems reasonable, because
/>>/  />/  the water box is large ~9.0 nm.
/>>/  /
/>>/  Now your outcome makes a lot more sense.  You're decoupling all of
the />>/  solvent? I don't see how that is going to be physically stable or
terribly /


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] converging value

[chris . neale]

You need to look at the dependence of the statistical value (in this  
case the average pressure) with (a) increasing equilibration time that  
you discard prior to collecting pressure data, and (b) increasing  
production simulation time that you include in your average.


Create the first plot, and discard any data in the early part of the  
simulation for which there appears to be a drift (the discarded data  
will not be used in the second plot). Then create the second plot,  
where you should see large oscillations about an average value and  
these oscillations should get smaller as the simulation time is  
increased. You can gauge the remaining error by block averaging or  
(much better) comparing with another simulation that you ran. There  
are also other statistical methods available, but I prefer the  
simplicity of these.


If your question is simply that you did these things but the  
oscillations remain large and you want to know how to get a more  
precise value, then the answer is simply to simulate longer.


Chris.

-- original message --

Dear All

How can I determine the converged value of a simulation?
Because the pressure has big oscilations around it's converging value but it
is difficult to determine that value.

can everybody guide me?
Thanks in advance
-- next part --
An HTML attachment was scrubbed...
URL: http://lists.gromacs.org/pipermail/gm

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] rationale behind tcoupling ligand with protein instead of with SOL

[chris . neale]

You didn't state your usage, but if you're doing US or decoupling, etc  
(some method where you lose the dynamics anyway) I suggest that you  
use Langevin dynamics. You will get the correct ensemble. Separate  
temperature coupling groups is a trick that helps in some cases, but  
it still does not give you the correct ensemble.


Chris.

-- original message --

Any rationale behind the thermostat coupling of a ligand with the protein
instead of the ligand with the solvent (as shown in Justin's T4 Lysozyme
binding example)? Especially with small drug-type molecules as generally the
ligand might/would take the usual place of solvent within a binding region
or other solvent accessible surface and it might be more realistic to
simulate the ligand as part of the solvent's ensemble (one might run two
simulations in order to compare the ligand-protein interaction to the
water-protein interaction at the interaction interface...)

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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 arising from T- and P-coupling methods

[chris . neale]

I can't tell you if there is a problem or not.

The only intended difference that I can see is that for 4.5.X:

# grompp by default sets the new nsttcouple parameter equal to  
nstlist, this means T-coupling is done less frequently; grompp checks  
if tau_t is large enough
# grompp by default sets the new nstpcouple parameter equal to  
nstlist, this means P-coupling is done less frequently; grompp checks  
if tau_p is large enough


(see http://www.gromacs.org/About_Gromacs/Release_Notes/Versions_4.5.x )

Note that there were some fixes to P-R scaling in the 4.0.X series --  
se http://www.gromacs.org/About_Gromacs/Release_Notes/Revisions_in_4.0


To test this further, I suggest that you need to define a good  
procedure to ensure that what you are seeing is (b) from a well  
equilibrated system and is also (b) statistically significant.


I suggest that one way to do this is to cycle through your T- and P-  
coupling options at a given temperature. For example, at T=500 K, run  
X ns of Berendsen, then X ns of Parrinello-Rahman, then back to  
Berendsen, and so on for a few cycles. Each time you switch coupling  
method, be sure to use the structure output from the previous run.  
Also be sure that X ns is as long as you can afford. This way, you  
will be able to pick out systematic changes as temperature/density  
oscillations with periods that are related to your changes of  
algorithm. This also ensures that what you are seeing is not simply an  
artifact of having a poorly equilibrated density in your initial  
structure.


You could also run the same cycle with Berendsen and V-rescale.

Chris.

-- original message --

Dear gmxers,
   According to my recent practice, we find that the Berensen methods  
for T- and P- coupling can yield reasonable averaged density as a  
function of temperature, but when the v-rescale method and the  
Parrinello-Rahman method are employed for T- and P- coupling, somewhat  
unexpected results (i.e. density at higher temperature is bigger than  
that at lower tempearature) are obtained. Generally, the latter setup  
is considered to be prefered to the former one in simulating realistic  
ensemble. I am using gmx-4.5.3, and previously I have also performed  
one similar work using 4.0 which can generate expected results using  
the latter setup. I wonder if this version 4.5.3 has some bugs in  
calculating T and P, and are they dealt with in 4.5.4? Please give me  
some hints.


 Yours sincerely,
 Chaofu Wu, Dr.
  --
  Department of Chemistry and Materials Science, Hunan University of  
Humanities, Science and Technology, Loudi 417000, the People?s  
Republic of China (P.R. China)

-- next part --
An HTML attachment was scrubbed...
URL:  
http://lists.gromacs.org/pipermail/gmx-users/attachments/20110412/ceadf085/attachment.html



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Different TI free energy values in 4.0.7 and 4.5.3

[chris . neale]

Dear Yan:

If you want to get to the bottom of this then (a) please report your  
mdp files and (b) show some evidence that your values are converged  
by, for instance, block averaging the run-time.


Chris.

-- original message --

Dear Gromacs users,

  I got different free energy values by Thermodynamic Integration (TI) in
Gromacs 4.0.7 and 4.5.3.

 When I calculated the hydration free energy of Chloride ion by using
Thermodynamic Integration method, the results for the free energy of
un-charging the Chloride ion from q=-e to 0 in SPC water are 359+/-6 kJ/mol
in Gromacs 4.0.7 and 385 +/-6 kJ/mol in Gromacs 4.5.3, respectively.

  The literature value is 392 kJ/mol (J. Phys. Chem. 1996, 100, 1206-1215,
Table 6), which is close to the result in Gromacs 4.5.3 but not in 4.0.7.

  As a check, Gromacs 3.3.1 gave the result as 381+/-3 kJ/mol.

  I searched the mail list and the bugzilla of Gromacs but didn't find
explanation for this issue of version.

  Does anyone know the reason for the different TI free energy values in
those different version of 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] Different TI free energy values in 4.0.7 and 4.5.3

[chris . neale]

Based on your use of:

sc-sigma = 0.3

I wonder what happens with 4.5.3 when you set the env.var.  
GMX_SCSIGMA_MIN to 0


quoting http://www.gromacs.org/About_Gromacs/Release_Notes/Versions_4.5.x

"for free-energy calculations sc-sigma now also sets the minimum  
soft-core sigma (old tpr files retain the old behavior, which can be  
enforced by setting the env.var. GMX_SCSIGMA_MIN to 0)"


-- original message --

Dear Gromacs users,

   One example of my mdp files with lambda=0.5 is like the following:

integrator   = md
tinit= 0
dt   = 0.004
nsteps   = 25
simulation_part  = 1
comm_mode= Linear
nstcomm  = 10


nstlist  = 2
ns_type  = grid
pbc  = xyz
periodic_molecules   = no
rlist= 1.4


coulombtype  = PME
rcoulomb-switch  = 0
rcoulomb = 1.4
epsilon_r= 1
epsilon_rf   = 1
vdw_type = cut-off
rvdw = 1.4
fourierspacing   = 0.168
pmeorder = 4
ewald_rtol   = 1.0e-5
ewald_geometry   = 3d
optimize_fft = yes


tcoupl   = v-rescale
tc-grps  = ioq !ioq
tau_t= 0.1 0.1
ref_t= 298.15 298.15
pcoupl   = Berendsen
Pcoupltype   = Isotropic
tau_p= 1.0
compressibility  = 4.5e-5
ref_p= 1.0


constraints  = all-bonds
constraint_algorithm = lincs
lincs_order  = 4
lincs_iter   = 1

free_energy  = yes
init_lambda  = 0.5
delta_lambda = 0
foreign_lambda   =
sc-alpha = 0.7
sc-power = 1
sc-sigma = 0.3
nstdhdl  = 10
separate-dhdl-file   = yes
dhdl-derivatives = yes
dh_hist_size = 0
dh_hist_spacing  = 0.1
couple-moltype   =
couple-lambda0   = vdw-q
couple-lambda1   = vdw-q
couple-intramol  = no

The system has 1 Chloride ion with ffG43a1 force field parameter and SPC
water molecules with heavy hydrogen in a cubic box of size 3 nm. I have
sampled 21 lambda windows ranging from 0 to 1.

   The values are converged by block averaging.

   Regards,
   Yan

On Thu, Apr 14, 2011 at 3:35 PM,  wrote:


Dear Yan:

If you want to get to the bottom of this then (a) please report your mdp
files and (b) show some evidence that your values are converged by, for
instance, block averaging the run-time.

Chris.

-- original message --

Dear Gromacs users,

 I got different free energy values by Thermodynamic Integration (TI) in
Gromacs 4.0.7 and 4.5.3.

 When I calculated the hydration free energy of Chloride ion by using
Thermodynamic Integration method, the results for the free energy of
un-charging the Chloride ion from q=-e to 0 in SPC water are 359+/-6 kJ/mol
in Gromacs 4.0.7 and 385 +/-6 kJ/mol in Gromacs 4.5.3, respectively.

 The literature value is 392 kJ/mol (J. Phys. Chem. 1996, 100, 1206-1215,
Table 6), which is close to the result in Gromacs 4.5.3 but not in 4.0.7.

 As a check, Gromacs 3.3.1 gave the result as 381+/-3 kJ/mol.

 I searched the mail list and the bugzilla of Gromacs but didn't find
explanation for this issue of version.

 Does anyone know the reason for the different TI free energy values in
those different version of 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] Cylinder Pulling output file pullx.xvg wron

[chris . neale]

That part of the pull code has a different printf statement in version  
4.5.3. Can you test 4.0.7 output vs. 4.5.3 output and see if 4.5.3 is  
correct? You could simple use nsteps=0 genvel=no and load a .gro with  
no velocities in your mdrun.


Chris.

-- original message --

I am trying cylinder pulling as an option of umbrella sampling using gromacs
4.0.7. Things seem to work but I noticed that in the output file pullx.xvg,
the referece group's postion is always 0. I don't think it is right. When I
use position pulling, the reference position is never 0.

Thanks for you help.
Kun

Below is my .mdp file
pull= umbrella
pull_geometry   = cylinder
pull_r1 = 1.0
Pull_r0 = 1.0
pull_dim= N N Y
pull_start  = yes
pull_init1  = 0.0
pull_ngroups= 1
pull_group0 = group0
pull_group1 = group1
pull_rate1  = 0.01
pull_vec1   = 0.0 0.0 -1.0
pull_k1 = 3000  ; kJ mol^-1 nm^-2
pull_nstxout= 1000
pull_nstfout= 1000

Here is the output:

@title "Pull COM"
@xaxis  label "Time (ps)"
@yaxis  label "Position (nm)"
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "0 Z"
@ s1 legend "1 dZ"
0.  0   5.69818
2.  0   5.67685
4.  0   5.68564
6.  0   5.65461
8.  0   5.63482
10. 0   5.61658
12. 0   5.61679
14. 0   5.54434
16. 0   5.55352
18. 0   5.4964
-- next part -

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Use different fudgeLJ and fudgeQQ values in simulations with the AMBER force field

[chris . neale]

Stephane:

The problem is the coulombic 1-4 interactions. There is not (and never  
has been as far as I can tell) a way to get different 1-4 coulombic  
scaling for different pairs. There is a trick that works when  
combining one ff that uses fudgeQQ 1.0 and another that uses fudgeQQ  
0.5:  
http://oldwww.gromacs.org/pipermail/gmx-users/2010-March/049428.html .


This trick can be applied to obtain micelle fudgeLJ 1.0 and fudgeQQ  
1.0 and, simultaneously, peptide fudgeLJ 0.5 and fudgeQQ = 5/6 (close  
but not exactly 0.8333 due to rounding) as follows:


1. set the fudgeLJ to 1.0 and fudgeQQ to 0.16

2. calculate the micelle [ pairtypes ] interactions according to the  
combination rules of the ff and then divide the values by 6. You can  
include these "values/6" in the pairtypes section or directly in the  
pairs section.


3. calculate the protein [ pairtypes ] interactions according to the  
combination rules of the ff and then divide the values by 10. You can  
include these "values/12" in the pairtypes section or directly in the  
pairs section.


4. Include all of the micelle [ pairs ] 6 times and all of the peptide  
[ pairs ] 5 times in their respective topologies


The energy should be nearly correct. The only problems will be (a)  
different internal rounding due to different order of operations (this  
will be very small and will get even smaller in double precision) and  
(b) based on the fact that 0.8333 is not exactly 5/6 in the first  
place (you're stuck with this one).


You will need to test that I got the values correctly above. To do  
that, you can compute the pair energy of micelle and protein systems  
individually using the standard ff files and then again using the  
modified version that is detailed above and that should work for both.


If you use this, an appropriate reference would be as an extension of  
the half-epsilon double-pairlist method:

Biophys J. 2010 Mar 3;98(5):784-792. "An Iris-Like
Mechanism of Pore Dilation in the CorA Magnesium Transport System."
Chakrabarti N, Neale C, Payandeh J, Pai EF, Pomès R.

Chris.

-- original message --

Dear gmx users,

Recently I have parametrized a new force field for glycolipid molecules
based on the GLYCAM parameters. I would like to use it with the
AMBER99SB-ILDN force field in gromacs for simulations of peptides and
glycolipds. As you know GLYCAM uses two different values for fudgeLJ and
fudgeQQ (e.g. set to 1.0) instead of 0.5 and 0.8333 values. So it is
possible to define two different fudgeLJ and fudgeQQ values in GROMACS
(as in AMBER11, for example) for the micelle and the peptides.

I have found an old response (2007) in the list archive for a old
version of gromacs

http://lists.gromacs.org/pipermail/gmx-users/2007-February/025907.html

Since this response have been posted 4 years ago, i am not sure that
this approach can be also use with the current version of gromacs ?

Thank you in advance for your response


Stephane





--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] lipid ff decisions

[chris . neale]

Has anybody done or read a good comparative study on what lipid ff  
works best in combination with proteins? There are so many options and  
I've realized that when anybody asks me what ff combination to use,  
all I can do is to tell them what combination I use and that so far it  
seems to work fine for me...


Thanks,
Chris.

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

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

[gmx-users] g_traj -ox -com slight inconsistency

[Chris Neale]

Dear Users:

I have noticed an inconsistency with g_traj -ox -com output when using 
-f a.gro and either -s a.gro or -s a.tpr where a.tpr was constructed 
from a.gro.


There does not appear to be any difference when not using the -com flag.

The difference is on the order of <=0.008 nm so I suppose that it could 
be related to rounding and order of operations or precision.


The same results exist for at least 4.0.5 and 4.5.3.

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

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

[gmx-users] g_select is pretty fantastic

[chris . neale]

Dear g_select developers:

This is a shout out of thanks for this new tool, which I find very useful.

For users who have not played around with this, I'd suggest trying it  
(especially if you are already familiar with VMD-style selection  
syntax).


I want to make special thanks for the VMD-style selection language,  
which avoids us all having to learn yet another protocol. I would  
suggest direct interaction between the g_select developers and the VMD  
developers to bring this into lock-step over the coming years (perhaps  
a universal selection syntax?). Another suggestion is that I only  
realized that the syntax is essentially a VMD-syntax by trying it out  
(there is no documentation that mentions this). I'm thus not sure how  
overlapping the syntaxes really are, but would request that they  
become fully-overlapping over time and that the current possibilities  
are mentioned in the help information.


Another important note is that you need to supply a .ndx file from a  
simple run through make_ndx to get things to work like "group  
Protein". I think it was a good decision to wall-off the creation of  
standard selection groups (make_ndx) from the g_select program, but I  
think that this fact would be useful in the documentation.


To be more explicit about the possibilities, one can do things like:

GRO=indo.1.2.6.5us.gro;
echo -e "aP8\nq\n" | make_ndx -f ${GRO} -o initial.ndx;
g_select -f ${GRO} -n initial.ndx -on select.ndx -select "same residue  
as group P8 and within 1 of group Protein" -s full.tpr;


Pretty cool.

Thanks again,
Chris.

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

[gmx-users] free energy

[chris . neale]

Your windows are restrained. The PMF that you get out of WHAM is a  
representation of the relative sampling after removing the umbrella  
biases. Sounds like you are saying that you look at the still-biased  
trajectories and you see different a different distribution of states  
than you do in the de-biased PMF. not sure what the problem is  
here. Perhaps go back to some review literature on US.


Now, if you saw more bound than unbound in unrestrained simulations,  
then that's a different story, but that doesn't appear to be what you  
are talking about.


Chris.

-- original message --

Hi all

I have performed a PMF simulation of taking part of a molecule out of
the cavity of a host using umbrella sampling. The free energy curve
suggests that the guest prefers to be outside the host as the bound
state is higher in energy, or the free energy difference to go in is
positive. However when I view the trajectories for each window it
appears that there is always more bound states than unbound. This leads
me to doubt my free energy profile?

Cheers

Gavin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] free energy

[chris . neale]

yes it will, and so will it affect the profile if water molecules or  
ions go in when both chains are absent.


You'll need to determine what question you are trying to answer and  
also think pretty hard about what your PMF really means in the context  
of this system.


Chris.

-- original message --


Hi Chris

My windows are restrained obviously using the force constant in the mdp
file. The trajectories that I have viewed are those of the individual
biased sampling windows. I have not put on the unbiased simulations yet.
There is also the following issue: The simulations involve two identical
molecules containing hydrocarbon chains. I calculate the PMF to take a
specific hydrocarbon chain of one molecule out of a specific site on the
neighbouring molecule. When this guest chain goes out, other chains can
go in (intramolecular or intermolecular). Will this affect the profile?

Gavin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] free energy

[chris . neale]

Dear Gavin:

I'm not seeing a lot of effort on your part to answer these questions.  
I suspect that you can answer some of this. Good luck!


Chris.

-- original message --

The current simulations are currently in vacuum. Does the following mdp
file seem ok. Note that this is a production run using the final
configuration after a lengthy equilibration. Also I was thinking about
trying to prevent the other chains from entering the specific site; is
there a a way to implement this in gromacs.

Gavin

chris.neale at utoronto.ca wrote:

yes it will, and so will it affect the profile if water molecules or
ions go in when both chains are absent.

You'll need to determine what question you are trying to answer and
also think pretty hard about what your PMF really means in the context
of this system.

Chris.

-- original message --


Hi Chris

My windows are restrained obviously using the force constant in the mdp
file. The trajectories that I have viewed are those of the individual
biased sampling windows. I have not put on the unbiased simulations yet.
There is also the following issue: The simulations involve two identical
molecules containing hydrocarbon chains. I calculate the PMF to take a
specific hydrocarbon chain of one molecule out of a specific site on the
neighbouring molecule. When this guest chain goes out, other chains can
go in (intramolecular or intermolecular). Will this affect the profile?

Gavin



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] difference of protein dynamics between NPT and NVE ensemble?

[chris . neale]

I'm not sure about the V vs. P part, but the E vs. T part depends on  
the thermostat that you use. For some thermostats, the problem can be  
large: http://www.ncbi.nlm.nih.gov/pubmed/20046980


-- original message --

Does anyone have experience of doing NVE simulation? I wish to know, for a
system with box size 6.1*6.1*6.1nm, whether the microcanonical ensemble for
a NVE simulation will converge with the canonical ensemble from a NPT
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] Count mismatch for state entry SDx, code count is 754728, file count is 0

[Chris Neale]

Dear Users:

Using gromacs 4.0.5, I find that there are at least some cases where 
some type of disk error can get propagated through both my.tpr and 
my_prev.tpr, complicating restarts. This used to be a bigger problem in 
gromacs 3, and I don't recall ever seeing it in gromacs 4 so I thought I 
would post a notification.


I'm just going to extract some coordinates and restart, but ideally this 
wouldn't happen. A google search for the relevant error "Count mismatch 
for state entry" only turns up some online source code.


I don't know if this error occurs in 4.5.3, and it's not binary 
reproducible so that would be difficult to check. Still, the error 
checking that regularly occurs prior to overwriting the previous (and 
without error) _prev.cpt file with a new (and with error) _prev.cpt file 
seemed to not catch this problem, at least with gromacs 4.0.5.


The run that wrote out the .tpr finished normally due to -maxh, with a 
stderr that looked like this:


... < snip > ...
starting mdrun 'Generated by genbox'
1000 steps,  2.0 ps (continuing from step 3769350,   7538.7 ps).
[gpc-f138n034:06165] 15 more processes have sent help message 
help-mpi-btl-base.txt / btl:no-nics
[gpc-f138n034:06165] Set MCA parameter "orte_base_help_aggregate" to 0 
to see all help / error messages


Step 5036590: Run time exceeded 47.322 hours, will terminate the run

Step 5036600: Run time exceeded 47.322 hours, will terminate the run

 Average load imbalance: 0.2 %
 Part of the total run time spent waiting due to load imbalance: 0.2 %
 Steps where the load balancing was limited by -rdd, -rcon and/or -dds: 
X 0 % Z 0 %

 Average PME mesh/force load: 0.745
 Part of the total run time spent waiting due to PP/PME imbalance: 4.9 %


Parallel run - timing based on wallclock.

   NODE (s)   Real (s)  (%)
   Time: 170485.000 170485.000100.0
   1d23h21:25
   (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:625.583 31.889  1.284 18.685

gcq#165: "I'm a Jerk" (F. Black)


gcq#165: "I'm a Jerk" (F. Black)

#

And then when I gmxcheck both of the .cpt files I get the exact same 
error, although the files do differ:


$ diff md1.cpt md1_prev.cpt
Binary files md1.cpt and md1_prev.cpt differ


$ gmxcheck  -f md1.cpt
 :-)  G  R  O  M  A  C  S  (-:

  S  C  A  M  O  R  G

:-)  VERSION 4.0.5  (-:


  Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
   Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 Copyright (c) 2001-2008, The GROMACS development team,
check out http://www.gromacs.org for more information.

 This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
 as published by the Free Software Foundation; either version 2
 of the License, or (at your option) any later version.

   :-)  gmxcheck  (-:

Option Filename  Type Description

  -fmd1.cpt  Input, Opt!  Trajectory: xtc trr trj gro g96 pdb cpt
 -f2   traj.xtc  Input, Opt.  Trajectory: xtc trr trj gro g96 pdb cpt
 -s1   top1.tpr  Input, Opt.  Run input file: tpr tpb tpa
 -s2   top2.tpr  Input, Opt.  Run input file: tpr tpb tpa
  -c  topol.tpr  Input, Opt.  Structure+mass(db): tpr tpb tpa gro 
g96 pdb

  -e   ener.edr  Input, Opt.  Energy file: edr ene
 -e2  ener2.edr  Input, Opt.  Energy file: edr ene
  -n  index.ndx  Input, Opt.  Index file
  -mdoc.tex  Output, Opt. LaTeX file

Option   Type   Value   Description
--
-[no]h   bool   no  Print help info and quit
-niceint0   Set the nicelevel
-vdwfac  real   0.8 Fraction of sum of VdW radii used as warning
cutoff
-bonlo   real   0.4 Min. fract. of sum of VdW radii for bonded atoms
-bonhi   real   0.7 Max. fract. of sum of VdW radii for bonded atoms
-tol real   0.001   Relative tolerance for comparing real values
defined as 2*(a-b)/(|a|+|b|)
-[no]ab  bool   no  Compare the A and B topology from one file
-lastenerstring Last energy term to compare (if not given 
all are

tested). It makes sense to go up until the
Pressure.

Checking file md1.cpt

---
Program gmxcheck, VERSION 4.0.5
Source code file: checkpoint.c, line: 186

Fatal error:
Count mismatch for state entry SDx, code count is 754728, file count is 0

---

"Confirmed" (Star Trek)

##

[gmx-users] Re: Count mismatch for state entry SDx, code count is 754728, file count is 0

[Chris Neale]

Apologies: my.tpr and my_prev.tpr should have read my.cpt and my_prev.cpt.

On 11-05-05 12:36 PM, Chris Neale wrote:

Dear Users:

Using gromacs 4.0.5, I find that there are at least some cases where 
some type of disk error can get propagated through both my.tpr and 
my_prev.tpr, complicating restarts. This used to be a bigger problem 
in gromacs 3, and I don't recall ever seeing it in gromacs 4 so I 
thought I would post a notification.


I'm just going to extract some coordinates and restart, but ideally 
this wouldn't happen. A google search for the relevant error "Count 
mismatch for state entry" only turns up some online source code.


I don't know if this error occurs in 4.5.3, and it's not binary 
reproducible so that would be difficult to check. Still, the error 
checking that regularly occurs prior to overwriting the previous (and 
without error) _prev.cpt file with a new (and with error) _prev.cpt 
file seemed to not catch this problem, at least with gromacs 4.0.5.


The run that wrote out the .tpr finished normally due to -maxh, with a 
stderr that looked like this:


... < snip > ...
starting mdrun 'Generated by genbox'
1000 steps,  2.0 ps (continuing from step 3769350,   7538.7 ps).
[gpc-f138n034:06165] 15 more processes have sent help message 
help-mpi-btl-base.txt / btl:no-nics
[gpc-f138n034:06165] Set MCA parameter "orte_base_help_aggregate" to 0 
to see all help / error messages


Step 5036590: Run time exceeded 47.322 hours, will terminate the run

Step 5036600: Run time exceeded 47.322 hours, will terminate the run

 Average load imbalance: 0.2 %
 Part of the total run time spent waiting due to load imbalance: 0.2 %
 Steps where the load balancing was limited by -rdd, -rcon and/or 
-dds: X 0 % Z 0 %

 Average PME mesh/force load: 0.745
 Part of the total run time spent waiting due to PP/PME imbalance: 4.9 %


Parallel run - timing based on wallclock.

   NODE (s)   Real (s)  (%)
   Time: 170485.000 170485.000100.0
   1d23h21:25
   (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:625.583 31.889  1.284 18.685

gcq#165: "I'm a Jerk" (F. Black)


gcq#165: "I'm a Jerk" (F. Black)

#

And then when I gmxcheck both of the .cpt files I get the exact same 
error, although the files do differ:


$ diff md1.cpt md1_prev.cpt
Binary files md1.cpt and md1_prev.cpt differ


$ gmxcheck  -f md1.cpt
 :-)  G  R  O  M  A  C  S  (-:

  S  C  A  M  O  R  G

:-)  VERSION 4.0.5  (-:


  Written by David van der Spoel, Erik Lindahl, Berk Hess, and 
others.

   Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 Copyright (c) 2001-2008, The GROMACS development team,
check out http://www.gromacs.org for more information.

 This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
 as published by the Free Software Foundation; either version 2
 of the License, or (at your option) any later version.

   :-)  gmxcheck  (-:

Option Filename  Type Description

  -fmd1.cpt  Input, Opt!  Trajectory: xtc trr trj gro g96 pdb cpt
 -f2   traj.xtc  Input, Opt.  Trajectory: xtc trr trj gro g96 pdb cpt
 -s1   top1.tpr  Input, Opt.  Run input file: tpr tpb tpa
 -s2   top2.tpr  Input, Opt.  Run input file: tpr tpb tpa
  -c  topol.tpr  Input, Opt.  Structure+mass(db): tpr tpb tpa gro 
g96 pdb

  -e   ener.edr  Input, Opt.  Energy file: edr ene
 -e2  ener2.edr  Input, Opt.  Energy file: edr ene
  -n  index.ndx  Input, Opt.  Index file
  -mdoc.tex  Output, Opt. LaTeX file

Option   Type   Value   Description
--
-[no]h   bool   no  Print help info and quit
-niceint0   Set the nicelevel
-vdwfac  real   0.8 Fraction of sum of VdW radii used as warning
cutoff
-bonlo   real   0.4 Min. fract. of sum of VdW radii for bonded 
atoms
-bonhi   real   0.7 Max. fract. of sum of VdW radii for bonded 
atoms

-tol real   0.001   Relative tolerance for comparing real values
defined as 2*(a-b)/(|a|+|b|)
-[no]ab  bool   no  Compare the A and B topology from one file
-lastenerstring Last energy term to compare (if not given 
all are

tested). It makes sense to go up until the
Pressure.

Checking file md1.cpt

---
Program gmxcheck, VERSION 4.0.5
Source code file: che

[gmx-users] constraints

[Chris Neale]

generally, look at mdout.mdp from grompp

-- original message --

In my input file if I don't specify "constraints" then

What is default

constraints=none

constraints=all bonds


NIlesh




--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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 potential energy=nan for r=0 nm and charged atoms

[chris . neale]

Dear Users:

I am trying to work with tabulated potentials for the first time. I  
would like to allow 2 charges to be at the same spot. Imagine taking a  
sodium and chloride ion and removing the LJ entirely and I want the  
system to evaluate energies correctly and be stable during the  
simulation. It was my intention to add a sort-of soft-core to the  
charge-charge interactions at very close distances so that the system  
did not crash. However, my testing didn't even get that far. I find  
that the Coulomb(SR) energy = nan even when using tables that suggest  
it should be zero.


I suppose that there is some other place in the code where this  
problem develops?


While looking at share/gromacs/top/table6-12.xvg it seems that this  
should not lead to energies of nan, but it does.


$ head share/gromacs/top/table6-9.xvg
#
# Table AB1: ndisp=6 nrep=9
#
0.00e+00   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
2.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
4.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
6.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
8.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
1.00e-02   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
1.20e-02   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00


and therefore, according to manual 4.5.4 page 151 equation 6.23, the  
[(qi*qj)/(4*pi*epsilon)]*f(rij) coulomb SR value should equal zero. I  
tried a few things, like modifying the values so that they are  
constant from 0 to 0.04 nm in the table.xvg file, but this did not help.


 here is the .top file that I used along with  
ffoplsaa to test out this idea:


#include "ffoplsaa.itp"
#include "na.itp"

[ system ]
; name
tables

[ molecules ]
; name  number
MOL_A  1
MOL_B  1

 and here is the .itp file that I used along with  
ffoplsaa to test out this idea:


[ atomtypes ]
 aaa   AA  8 15.999400.0   A0.0  0.0

[ moleculetype ]
; molname   nrexcl
MOL_A 1

[ atoms ]
; idat type   res nr  residu name at name cg nr   charge
1   aaa   1   MLA AA  11.0

[ moleculetype ]
; molname   nrexcl
MOL_B 1

[ atoms ]
; idat type   res nr  residu name at name cg nr   charge
1   aaa   1   MLB BB  1-1.0

 And the .gro file:

title
2
1MLA AA1   0.000   0.000   0.000
1MLB BB1   0.000   0.000   0.000
  5.0  5.0  5.0

### And my .mdp file is

nsteps = 0
ns_type = grid
rlist = 1
rcoulomb = 1
rvdw = 1
coulombtype = user
tc_grps =  System
tau_t   =  1.0
ld_seed =  -1
ref_t = 300
gen_temp = 300
gen_vel = yes
unconstrained_start = no
gen_seed = -1
Pcoupl = berendsen
pcoupltype = isotropic
tau_p = 4
compressibility = 4.5e-5
ref_p = 1.0



Still, when I run a 0-step mdrun (mdrun -deffnm md1 -table  
mod2.table.xvg), I get:


   Energies (kJ/mol)
LJ (SR)   Coulomb (SR)  PotentialKinetic En.   Total Energy
0.0e+00nannannannan
Temperature Pressure (bar)
nannan

   # I get the same thing if the two charges are like or dislike...

But if I then modity my topologies to turn the +1 charge to zero, I  
get sensible output:


LJ (SR)   Coulomb (SR)  PotentialKinetic En.   Total Energy
0.0e+000.0e+000.0e+005.27174e+005.27174e+00
Temperature Pressure (bar)
4.22694e+024.66876e-01

Also, if I leave the charges on and separate the two ions by 0.001 nm,  
I also don't get the nan values.


   Energies (kJ/mol)
LJ (SR)   Coulomb (SR)  PotentialKinetic En.   Total Energy
0.0e+000.0e+000.0e+006.64240e-016.64240e-01
Temperature Pressure (bar)
5.32595e+015.88265e-02


Thank you for your assistance,
Chris.



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

[gmx-users] beware gen-vel=no

[chris . neale]

Regarding gen-vel=no as discussed here:
http://lists.gromacs.org/pipermail/gmx-users/2011-May/061151.html

I would caution against the general use of gen-vel=no and then  
coupling to a temperature of 300 K with the Berendsen thermostat. One  
problem that can arise is concerted unfolding. Temperature is random  
velocity, but with gen-vel=no and a structure that is just slightly  
too compact, you will have a slight net force that gets massively  
scaled up by the coupling algorithm and results, not in true  
temperature, but in large scale coordinated motion.


It's just something to look out for. gen-vel=no should not be a  
problem with the sd algorithm or with some type of position restrained  
equilibration procedure or, of course, if loading in velocities from  
some cpt or trr file.


Chris.

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

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

[gmx-users] tabulated potential energy=nan for r=0 nm and charged atoms

[chris . neale]

Dear users:

with assistance from Berk on the developers list, I am posting a  
work-around to this issue.


the nonbonded interactions calculate the nonbonded distance, r, using  
the gmx_invsqrt function. Thus, the optimized kernel from 4.5.4 and  
4.0.7 both give nan when two charges occupy the space, even when using  
tables that specify the force and potential energy to be zero.


I had expected that I would need to modify the generic kernel source,  
but that turned out to be unnecessary -- just using it was enough.


For some reason, using "export GMX_NOOPTIMIZEDKERNELS=1" prior to  
running mdrun solves the issue for 4.5.4 but not 4.0.7 (on my  
architecture). To solve the problem in v4.0.7 (and this also works for  
v4.5.4) I must "export GMX_NB_GENERIC=1; export GMX_NO_SOLV_OPT=1"  
prior to running mdrun. Note that the generic nonbonded kernel runs  
1.1x slower then the unoptimized kernel, which itself runs 2.4x slower  
than the optimized kernel (on my architecture for this system).


Thank you,
Chris.

-- original message --

Dear Users:

I am trying to work with tabulated potentials for the first time. I
would like to allow 2 charges to be at the same spot. Imagine taking a
sodium and chloride ion and removing the LJ entirely and I want the
system to evaluate energies correctly and be stable during the
simulation. It was my intention to add a sort-of soft-core to the
charge-charge interactions at very close distances so that the system
did not crash. However, my testing didn't even get that far. I find
that the Coulomb(SR) energy = nan even when using tables that suggest
it should be zero.

I suppose that there is some other place in the code where this
problem develops?

While looking at share/gromacs/top/table6-12.xvg it seems that this
should not lead to energies of nan, but it does.

$ head share/gromacs/top/table6-9.xvg
#
# Table AB1: ndisp=6 nrep=9
#
0.00e+00   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
2.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
4.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
6.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
8.00e-03   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
1.00e-02   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00
1.20e-02   0.00e+00 0.00e+00
0.00e+00 0.00e+00   0.00e+00 0.00e+00

and therefore, according to manual 4.5.4 page 151 equation 6.23, the
[(qi*qj)/(4*pi*epsilon)]*f(rij) coulomb SR value should equal zero. I
tried a few things, like modifying the values so that they are
constant from 0 to 0.04 nm in the table.xvg file, but this did not help.

 here is the .top file that I used along with
ffoplsaa to test out this idea:

#include "ffoplsaa.itp"
#include "na.itp"

[ system ]
; name
tables

[ molecules ]
; name  number
MOL_A  1
MOL_B  1

 and here is the .itp file that I used along with
ffoplsaa to test out this idea:

[ atomtypes ]
  aaa   AA  8 15.999400.0   A0.0  0.0

[ moleculetype ]
; molname   nrexcl
MOL_A 1

[ atoms ]
; idat type   res nr  residu name at name cg nr   charge
1   aaa   1   MLA AA  11.0

[ moleculetype ]
; molname   nrexcl
MOL_B 1

[ atoms ]
; idat type   res nr  residu name at name cg nr   charge
1   aaa   1   MLB BB  1-1.0

 And the .gro file:

title
 2
 1MLA AA1   0.000   0.000   0.000
 1MLB BB1   0.000   0.000   0.000
   5.0  5.0  5.0

### And my .mdp file is

nsteps = 0
ns_type = grid
rlist = 1
rcoulomb = 1
rvdw = 1
coulombtype = user
tc_grps =  System
tau_t   =  1.0
ld_seed =  -1
ref_t = 300
gen_temp = 300
gen_vel = yes
unconstrained_start = no
gen_seed = -1
Pcoupl = berendsen
pcoupltype = isotropic
tau_p = 4
compressibility = 4.5e-5
ref_p = 1.0



Still, when I run a 0-step mdrun (mdrun -deffnm md1 -table
mod2.table.xvg), I get:

Energies (kJ/mol)
 LJ (SR)   Coulomb (SR)  PotentialKinetic En.   Total Energy
 0.0e+00nannannannan
 Temperature Pressure (bar)
 nannan

# I get the same thing if the two charges are like or dislike...

But if I then modity my topologies to turn the +1 charge to zero, I
get sensible output

[gmx-users] Rigid structure and flexible structure present in sing system

[chris . neale]

If I read between the lines correctly, you know how to do this in  
gromacs, but you wish that you got a big speedup from freezing most of  
the atoms in your system. If that is the case, then I think that  
gromacs can not help you in its current form. Therefore, I suggest  
that you try the Charmm software. It may be slower overall, but it  
gives you the speedup you expect when you freeze 95% of the atoms in  
your system so it may be much faster for your usage. Other software  
may also do this, but I have no experience with it. If you have  
explicit water, this won't be of much use since you're still only  
freezing <10% of the atoms in your system.


If you simply want to know how to freeze one structure, then use  
freeze groups and energy exclusions without pressure coupling, or use  
position restraints and refcoordscaling=com with position restraints,  
or create an elastic network of restraints.


Chris.

-- original message --

Hi,

For a while I've been looking for a way to include a rigid
(macromolecular) structure and a flexible (small-protein) structure in
a single system and have not had much success.

I looked for a way to apply a different set of constraints to each
structure, but couldn't find one. I looked into modifying the topology
file, but the only method I could find was to artificially increase
the mass of every atom in the rigid structure and this does not save
any computing time.

Does anyone know how I can fix the relative co-ordinates of one
structure and only calculate the MD of another?

Any constructive criticism (even if to inform me that I'm doing
something ridiculous) is very gratefully received.

Kind regards,
Luke



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] How to get angle distribution between two tyrosine stacking residues

[chris . neale]

1. Determine some mathematical calculation that you want to apply to  
some atoms.

2. make a .ndx group that includes all of those atoms
3. Run g_traj -ox to output the coordinates of those atoms
4. Apply your mathematical calculation using awk.

-- original message --

Can anyone tell how to calculate the normal-normal angle between two
stacking tyrosine residues as a function of time.


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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 during energy minimization

[chris . neale]

1. you're not doing energy minimization as your title suggests, please  
use a better title:

   integrator = md

2. generally, when you have a problem it is useful to try to simplify  
the system (e.g. get rid of those tables) and try to reproduce the  
problem with the simplest system possible before posting.


3. I know that you commented it out, but the following command looks  
like a typo or a misunderstanding... perhaps you meant -DFLEXIBLE ?

   ;define= -DEFLEXIBLE

4. Most relevant to you: This is a little ambitious, don't you think?
   dt = 0.02 ; ps !

Chris.

-- original message --



Dear Gmx users,
  I am trying to do  minimization of my system .i
have no problem wehen i grompp it but when i do the mdrun its giving me
some  segmentation error.I had attached the output of grompp and  mdrun
below.Any suggestions please.Thanks in advance
*OUTPUT OF GROMPP*

Ignoring obsolete mdp entry 'title'
Ignoring obsolete mdp entry 'cpp'

Back Off! I just backed up mdout.mdp to ./#mdout.mdp.2#
Generated 332520 of the 332520 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 332520 of the 332520 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'Ion'
Excluding 2 bonded neighbours molecule type 'SOL'
Analysing residue names:
There are: 2Ion residues
There are:   875  Water residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting
into groups...
Number of degrees of freedom in T-Coupling group rest is 5253.00
Largest charge group radii for Van der Waals: 0.039, 0.039 nm
Largest charge group radii for Coulomb:   0.085, 0.085 nm
This run will generate roughly 7 Mb of data

Back Off! I just backed up em.tpr to ./#em.tpr.2#





*OUTPUT OF MDRUN*

Reading file em.tpr, VERSION 4.5.3 (single precision)
Starting 2 threads

WARNING: For the 1498 non-zero entries for table 2 in table_Na_Cl.xvg the
forces deviate on average -2147483648% from minus the numerical derivative
of the potential


WARNING: For the 1498 non-zero entries for table 2 in table_Na_Cl.xvg the
forces deviate on average -2147483648% from minus the numerical derivative
of the potential

Making 1D domain decomposition 2 x 1 x 1
starting mdrun 'NA SODIUM ION in water'
1 steps,200.0 ps.

step 0: Water molecule starting at atom 2553 can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.

step 0: Water molecule starting at atom 2595 can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
Wrote pdb files with previous and current coordinates
Segmentation fault

*EM.MDP*
title = nacl
cpp = /usr/bin/cpp ; the c preprocessor
;define= -DEFLEXIBLE
integrator = md
dt = 0.02 ; ps !
nsteps = 1
nstlist = -1
coulombtype = user
energygrps = Na Cl Sol
energygrp_table = Na Cl
vdwtype= user
ns_type = grid
rlist = 1.4
rcoulomb = 1.0
rvdw = 1.0
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
pbc=xyz;
; Energy minimizing stuff
emtol = 1000.0
emstep = 0.01

regards,
sree
-- next part --


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] binaries for mac os x 10.6

[chris . neale]

Dear Users:

Does somebody have a set of gromacs binaries for mac os X 10.6.7 that  
they can post somewhere (perhaps in  
http://www.gromacs.org/Downloads/User_contributions/Gromacs_binaries)?


Macbook Air 3.1
Version 10.6.7
1.4 GHz Intel Core 2 Duo
2 GB 1067 MHz DDR3
SMC version 1.67f3


I am trying to install it on a macbook air and I can't imagine using  
up 4 GB of the 64 GB hard drive to install xcode and do an install  
from source. Even tools like fink seem to require xcode.


I'm very new to the mac (it's not even my mac) so if anybody has  
better ideas for getting a gromacs installation then that would be  
great.


Please note that I did try to search the archives, but for some reason  
a search for the word "gromacs" turned up many hits but a search for  
"mac" or "max os X 10.6" or variation therein turned up no hits...  
possibly because the words are so short? A google serch turned up a  
bit, but not what I am looking for.


Thank you very much for your time,
Chris.

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

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

[gmx-users] Why does the -append option exist?

[chris . neale]

Dear Dimitar:

to me, your posts have indeed appeared very aggressive. So much so  
that there are probably many people who can help who have decided not  
to bother due to what they perceive your tone to be.


On to the problem at hand: I agree with the previous suggestion that  
you may have somehow had two runs going in the same directory and  
competing to write to the same file.


Can you test this by submitting one script, then waiting for an hour,  
then submitting the same script again in the same directory? If you  
can reproduce the problem with this test then it would suggest that  
this is indeed one possible source of the problem, however it arrived.


I probably won't respond on this topic again, as I don't want to find  
myself involved in an argument on the list. I simply wanted to point  
out that there has been at least one valid suggestion that you are in  
a good position to test.


Very generally, when you mail the list, people are going to try to  
suggest possible sources of the problem and then, ideally, you can  
test those and report back.


Chris.


 Here is an example:

 
 [dpachov]$ ll -rth run1*  \#run1*
-rw-rw-r-- 1 dpachov dpachov  11K May  2 02:59 run1.po.mdp
-rw-rw-r-- 1 dpachov dpachov 4.6K May  2 02:59 run1.grompp.out
-rw-rw-r-- 1 dpachov dpachov 3.5M May 13 19:09 run1.gro
-rw-rw-r-- 1 dpachov dpachov 2.3M May 14 00:40 run1.tpr
-rw-rw-r-- 1 dpachov dpachov 2.3M May 14 00:40 run1-i.tpr
-rw-rw-r-- 1 dpachov dpachov0 May 29 21:53 run1.trr
-rw-rw-r-- 1 dpachov dpachov 1.2M May 31 10:45 run1.cpt
-rw-rw-r-- 1 dpachov dpachov 1.2M May 31 10:45 run1_prev.cpt
-rw-rw-r-- 1 dpachov dpachov0 Jun  3 14:03 run1.xtc
-rw-rw-r-- 1 dpachov dpachov0 Jun  3 14:03 run1.edr
-rw-rw-r-- 1 dpachov dpachov  15M Jun  3 17:03 run1.log
 

 Submitted by:

ii=1
ifmpi="mpirun -np $NSLOTS"

   if [ ! -f run${ii}-i.tpr ];then
   cp run${ii}.tpr run${ii}-i.tpr
  tpbconv -s run${ii}-i.tpr -until 20 -o run${ii}.tpr
   fi

k=`ls md-${ii}*.out | wc -l`
   outfile="md-${ii}-$k.out"
   if [[ -f run${ii}.cpt ]]; then

   $ifmpi `which mdrun` -s run${ii}.tpr -cpi run${ii}.cpt -v -deffnm
run${ii} -npme 0 > $outfile  2>&1

fi
 =


This script is not using mdrun -append.



-append is the default, it doesn't need to be explicitly listed.



Your original post suggested the use of -append was a problem. Why aren't
we seeing a script with mdrun -append? Also, please provide the full script
- it looks like there might be a loop around your tpbconv-then-mdrun
fragment.



There is no loop; this is a job script with PBS directives. The header of it
looks like:
===
#!/bin/bash
#$ -S /bin/bash
#$ -pe mpich 8
#$ -ckpt reloc
#$ -l mem_total=6G
===

as usual submitted by:

qsub -N  myjob.q





Note that a useful trouble-shooting technique can be to construct your
command line in a shell variable, echo it to stdout (redirected as suitable)
and then execute the contents of the variable. Now, nobody has to parse a
shell script to know what command line generated what output, and it can be
co-located with the command's stdout.



I somewhat understand your point, but could give an example if you think it
is really necessary?

  


 Writing checkpoint, step 51879590 at Tue May 31 10:45:22 2011
Energies (kJ/mol)
U-BProper Dih.  Improper Dih.  CMAP Dih.  LJ-14
8.33208e+034.72300e+035.31983e+02   -1.21532e+032.89586e+03
 Coulomb-14LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.
 3.00900e+049.31785e+04   -3.87790e+03   -7.40841e+05
-8.36838e+04
  PotentialKinetic En.   Total EnergyTemperature Pres. DC (bar)
   -6.89867e+051.28721e+05   -5.61146e+052.99472e+02   -1.24229e+02
 Pressure (bar)   Constr. rmsd
   -1.03491e+022.99840e-05
 


So the -append restart looks like it did fine here.


 Last output files from restarts:

 [dpachov]$ ll -rth md-1-*out | tail -10
-rw-rw-r-- 1 dpachov dpachov 6.1K Jun  3 16:40 md-1-2428.out
 -rw-rw-r-- 1 dpachov dpachov 6.2K Jun  3 16:44 md-1-2429.out
-rw-rw-r-- 1 dpachov dpachov 5.9K Jun  3 16:46 md-1-2430.out
-rw-rw-r-- 1 dpachov dpachov 5.9K Jun  3 16:48 md-1-2431.out
-rw-rw-r-- 1 dpachov dpachov 6.1K Jun  3 16:50 md-1-2432.out
-rw-rw-r-- 1 dpachov dpachov0 Jun  3 16:52 md-1-2433.out
-rw-rw-r-- 1 dpachov dpachov 6.2K Jun  3 16:55 md-1-2434.out
-rw-rw-r-- 1 dpachov dpachov 6.2K Jun  3 16:58 md-1-2435.out
-rw-rw-r-- 1 dpachov dpachov 5.9K Jun  3 17:03 md-1-2436.out
*-rw-rw-r-- 1 dpachov dpachov 5.8K Jun  3 17:04 md-1-2437.out*
 
+ around the time when the run1.xtc file seems to have been saved:

 [dpachov]$ ll -rth md-1-23[5-6][0-9]*out
-rw-rw-r-- 1 dpachov dpachov 6.2K J

[gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[chris . neale]

Dear Thomas:

I agree with your conclusions about g_order.

My usage of the Berger POPC lipids in combination with the g_order  
command (as per the instructions in  
http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order )  
yields order parameters of 0.193 and 0.051 for the first and second
carbons in the double bond at 300K. This is similar to but, as  
expected, slightly larger than the values described in  
http://www.sciencedirect.com/science/article/pii/S0014579305014365  
where they simulated at 310K and also used g_order. I suspect that one  
could find tens of publications that use g_order and have published  
similar values. The Poger paper (On the Validation of Molecular  
Dynamics Simulations of Saturated and cis-Monounsaturated  
Phosphatidylcholine Lipid Bilayers: A Comparison with Experiment)  
gives clearly different values for the first carbon in the double bond  
than is reported by g_order. One probable reason that this has gone  
unnoticed so long is that there appears to be no experimental oleoyl  
order parameters -- as far as I know -- if somebody finds them, please  
post them.


I took a look at the source and it is clearly incorrect for the i+1  
unsaturated partner in a double bond over i->i+1. I have developed a  
new version of g_order to fix this, but I'd like to test it first. If  
we could agree on a single .xtc file to analyze and compare our values  
then that would be great. I can provide you with one of mine (Berger  
lipids) or you can provide me with one of yours (off-list).


Thanks,
Chris.


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

[gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[chris . neale]

I have patched g_order version 4.5.4 to compute the order parameters  
in a different way. Despite trying, I don't understand the original  
g_order implementation. In the version that I am providing, the  
hydrogens are explicitly built and the coordinates of these hydrogens  
are output to stdout.


Compared to the original g_order: there are large differences for the  
2 carbons in the double bond and all single bonded carbons are the  
same within what appears to be rounding error.


Compared to the VMD script that Thomas Piggot mentioned, the results  
differ more, but only because that VMD script is using the real  
hydrogen positions and not some with idealized geometry. There is a  
larger difference for the double-bonded carbons, but if I use the  
idealized hydrogen positions (output by my new version of g_order)  
with the VMD .tcl script then I get the same values.


notes on the patch:
1. it is not complete for all options. Use it only with g_order -od  
and with or without -unsat.
2. There are extraneous print statement that will provide you with the  
coordinates of the atoms that are built. This is for debugging. Either  
run with g_order >/dev/null or comment out the 3 fprintf lines near  
the comments containing the word LOOK in capitals.

3. I did my best, but it could always be in error. Please use with caution.

To apply the patch (version 4.5.4)
cd src/tools
patch < g_order.patch

If it will not compile, then the spacing/hard-returns of the patch  
must have been mangled in the process of posting. Contact me off-list  
and I will provide you with a copy.


Here is the patch:

--- gmx_order.c 2011-03-04 06:10:44.0 -0500
+++ gmx_order.c 2011-06-09 17:57:11.586591000 -0400
@@ -379,6 +379,13 @@
real arcdist;
gmx_rmpbc_t  gpbc=NULL;

+// BEGIN CN MOD  


+// variables for modification by CN
+rvec ab,ac,bc,e,ce,eg,ef,g,f,cg,cf;
+float dab,def,deg,gdot,fdot,edot;
+// END CN MOD  


+
+
   /* PBC added for center-of-mass vector*/
   /* Initiate the pbc structure */
   memset(&pbc,0,sizeof(pbc));
@@ -501,98 +508,66 @@
direction[0],direction[1],direction[2]);*/
  }

-   if (bUnsat) {
- /* Using convention for unsaturated carbons */
- /* first get Sz, the vector from Cn to Cn+1 */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
- length = norm(dist);
- check_length(length, a[index[i]+j], a[index[i+1]+j]);
- svmul(1/length, dist, Sz);
-
- /* this is actually the cosine of the angle between the double bond
-and axis, because Sz is normalized and the two other components of
-the axis on the bilayer are zero */
- if (use_unitvector)
- {
-	sdbangle += acos(iprod(direction,Sz));  /*this can probably be  
optimized*/

- }
- else
-   sdbangle += acos(Sz[axis]);
-   } else {
- /* get vector dist(Cn-1,Cn+1) for tail atoms */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
- length = norm(dist);  /* determine distance between two atoms */
- check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
-
- svmul(1/length, dist, Sz);
- /* Sz is now the molecular axis Sz, normalized and all that */
-   }
-
-   /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
-  we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
-   rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
-   rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
-   cprod(tmp1, tmp2, Sx);
-   svmul(1/norm(Sx), Sx, Sx);
-
-   /* now we can get Sy from the outer product of Sx and Sz   */
-   cprod(Sz, Sx, Sy);
-   svmul(1/norm(Sy), Sy, Sy);
-
-   /* the square of cosine of the angle between dist and the axis.
-  Using the innerproduct, but two of the three elements are zero
-  Determine the sum of the orderparameter of all atoms in group
-  */
-   if (use_unitvector)
-   {
-	cossum[XX] = sqr(iprod(Sx,direction)); /* this is allowed, since Sa  
is normalized */

-   cossum[YY] = sqr(iprod(Sy,direction));
-   cossum[ZZ] = sqr(iprod(Sz,direction));
-   }
-   else
-   {
-   cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized 
*/
-   cossum[YY] = sqr(Sy[axis]);
-   cossum[ZZ] = sqr(Sz[axis]);
+// BEGIN CN MOD  


+  if (bUnsat) {
+//exactly as for the case without unsat, but here we select the  
ce vector to represent the hydrogen
+// place the 2 imaginary hydrogen atoms by enscribing a  
tetrahedron in a cube
+// a,b,c are the 3 backbone atoms, where a and b are in two cube  
corn

[gmx-users] gromacs 2.1

[chris . neale]

Dear users:

Does anybody have the source code for gromacs version 2.1? I would  
like to check the original source of g_order, but only versions 3.0  
and above are available on the gromacs website.


Thank you,
Chris.

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

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

[gmx-users] Gradually increasing the force constant during restrained MD runs

[chris . neale]

Dear Aditi:

you can do this by scripting gromacs to run in short segments,  
modifying the force constant for each run. e.g.:


#!/bin/bash
for((i=1;i<=100; i++)); do
  let "j=i-1"
  grompp -f ${i}.mdp -p ${i}.top -f ${j}.gro -o ${i}.tpr
  mdrun -deffnm ${i}
done

You only need to use some more scripting to create the ${i}.top and  
${i}.mdp files, although I think that orientational restraints are  
going to be part of the topology so the .mdp file can even remain  
constant.


You can then trjcat -settime once you are done.

-- original message --

I am running an MD simulation on RNA with orientation restraints.  
Instead of using a single force constant for the whole trajectory, I  
want to start with a very low force constant and increase it gradually  
and regularly up to a final suitably high value and then continue the  
simulation at this final value of force constant.


Can you suggest if there is a way to do so in GROMACS?

Thanks!

Aditi


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

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

[gmx-users] Question about umbrella sampling through a nanotube

[chris . neale]

Unfortunately, I don't think that there is any way to use this data  
(and only this data) to derive the PMF along z. The way that you did  
your US, the PMF along z is convoluted with xy motion. You can get the  
PMF along the reaction coordinate that you actually used (XYZ) using  
standard WHAM and try to interpret it, I suppose. In your special  
case, there is probably a mathematical conversion between the XYZ PMF  
and the Z-only PMF due to the symmetry of the system.


What one would desire is two separate restraints, one operating in Z  
and one operating in XY. That is, unfortunately, not possible in the  
standard version of gromacs.


-- original message --

Dear GMXers,

I'm trying to compute the PMF of a molecule along the centerline of a
nanotube (the axial direction of the nanotube is parallel to the z axis).
The nanotube is used as the reference group and the molecule as the pulling
group.

pull_geometry   = position
pull_dim= Y Y Y
pull_start  = no
pull_init1  = 0 0 1.2
pull_rate1  = 0.00
pull_k1 = 1000
pull_nstxout= 1000
pull_nstfout= 1000


"pull_dim=Y Y Y" is used to make sure that the molecule is always close to
the centerline of the nanotube.

From the output of "pullx.xvg", the COM positions in all three directions

(i.e. x,y,z) are included.
How to get the PMF only along z axis from here? Is it OK if I just take the
corresponding z column and do g_wham directly?

Thank you.

Best,
Yanbin
-- next part --


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] centering a bilayer at a cartesian coordinate

[chris . neale]

Dear users:

I am trying to process a trajectory so that a group of molecules has  
their center of mass at a constant position in Cartesian space. I have  
not been able to figure out how to do this.


The reason that I would like this is that I have conducted umbrella  
sampling of a solute along the normal to a lipid bilayer and I would  
like to create a movie of a false-trajectory of the immersion profile.  
Such processing is also important for generating SDFs and for use with  
g_density.


I expect that the problems that I am having are due to volume  
fluctuations in the NPT ensemble.


The best idea that I have tried was (gromacs 4.5.3):

trjconv -center -box 10 10 15

where the original dimensions are much less that 10 10 15

But still, when I run g_traj on the processed .gro file I do not get  
similar values for the COM.


In more detail:

echo 0 | trjconv -s md1.tpr -f md1_50ps_md1.part0001.xtc -o tmp.xtc -e  
10 -pbc mol

# select POPC for centering and System for output and use a .tpr with pbc=no
echo -e "13\n0\n" | trjconv -s empty.tpr -f tmp.xtc -box 10 10 15  
-center -pbc none -o tmp2.xtc

# select POPC for coordinate output
echo 13 | g_traj  -s empty.tpr -f tmp2.xtc -ox -com -nopbc

 99550  4.8424  5.00025 7.49519
 99600  4.79691 4.94304 7.39868
 99650  4.75789 4.85189 7.41581
 99700  4.74974 4.80423 7.48257
 99750  5.07858 5.04788 7.51654
 99800  4.99091 4.78976 7.47778
 99850  5.11406 4.84656 7.48769
 99900  5.09739 5.12395 7.47889
 99950  5.0411  4.88659 7.52737
 10 4.91271 5.07101 7.50397

For Z, the value is 7.48037+/-0.072562 (MIN=7.3 and MAX=7.8)

If I look at the frames from the min and max z values, they do clearly  
differ in the placement of the bilayer along z.


###

To simplify things, I redid this while only using the phosphorous  
atoms (1 per lipid) in the centering and coordinate output.


Here, Z=7.49323+/-0.0658727 (MIN=7.3 and MAX=7.7)

To verify this and since all the atoms have the same mass, I checked  
it with awk:


echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 47550 -b 4 -o look_low.gro
echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 80750 -b 8 -o  
look_high.gro

cat look_low.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.39197
cat look_high.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.62207

###

Finally, I rechecked with a single atom and here it worked exactly as  
expected. The atoms is placed at 5 5 7.5 in every single frame.


###

Then I checked using 2 phosphorous atoms from the same (or,  
separately, different) leaflets, and the result was the same:


$ tail coord.xvg
 99550  5   5   7.5
 99600  5   5.0005  7.5
 99650  5.0005  5.0005  7.5005
 99700  5   5.0005  7.5
 99750  5.0005  5   7.5005
 99800  5   5.0005  7.5005
 99850  5   5   7.5005
 99900  5.0005  5   7.5005
 99950  5.0005  5   7.4995
 10 5   5   7.5


Then using 3 phosphorous atoms from (containing atoms from both  
leaflets), and now there is more deviation:


$ tail coord.xvg
 99550  5.48833 5.30233 6.972
 99600  5.48733 5.29233 6.97167
 99650  5.50033 5.29567 6.93467
 99700  5.52133 5.27433 6.93967
 99750  5.489   5.253   6.99133
 99800  5.481   5.319   6.91533
 99850  5.47867 5.309   6.827
 99900  5.48733 5.292   6.83133
 99950  5.52767 5.276.86867
 10 5.518   5.319   6.89167

(MIN=6.7 and MAX=7.1)

I originally had problems with pbc during trjconv in which the -center  
option would place the bilayer at the peripheries of the unit cell in  
some frames (COM should still be at the center). This is what lead me  
to the route that I have used above (a .tpr that was generated with  
pbc=no and trjconv -pbc none -center). Nevertheless, I think that  
perhaps trjconv is using pbc information even in this case where ti  
should not.


##

If anybody can see what I am doing wrong or has a method to center a  
bilayer at a cartesian coordinate, then I am very interested.


Thank you for your time,
Chris.


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

[gmx-users] centering a bilayer at a cartesian coordinate

[chris . neale]

I solved it. It was my own error to assume that trjconv -center  
centers the COM. Actually, it centers the value of (min+max)/2 in each  
dimension.


I can modify trjconv locally to do what I want.

Thank you,
Chris.

static void center_x(int ecenter,rvec x[],matrix box,int n,int  
nc,atom_id ci[])

{
int i,m,ai;
rvec cmin,cmax,box_center,dx;

if (nc > 0) {
copy_rvec(x[ci[0]],cmin);
copy_rvec(x[ci[0]],cmax);
for(i=0; i cmax[m])
cmax[m]=x[ai][m];
}
}
calc_box_center(ecenter,box,box_center);
for(m=0; m
Dear users:

I am trying to process a trajectory so that a group of molecules has
their center of mass at a constant position in Cartesian space. I have
not been able to figure out how to do this.

The reason that I would like this is that I have conducted umbrella
sampling of a solute along the normal to a lipid bilayer and I would
like to create a movie of a false-trajectory of the immersion profile.
Such processing is also important for generating SDFs and for use with
g_density.

I expect that the problems that I am having are due to volume
fluctuations in the NPT ensemble.

The best idea that I have tried was (gromacs 4.5.3):

trjconv -center -box 10 10 15

where the original dimensions are much less that 10 10 15

But still, when I run g_traj on the processed .gro file I do not get
similar values for the COM.

In more detail:

echo 0 | trjconv -s md1.tpr -f md1_50ps_md1.part0001.xtc -o tmp.xtc -e
10 -pbc mol
# select POPC for centering and System for output and use a .tpr with pbc=no
echo -e "13\n0\n" | trjconv -s empty.tpr -f tmp.xtc -box 10 10 15
-center -pbc none -o tmp2.xtc
# select POPC for coordinate output
echo 13 | g_traj  -s empty.tpr -f tmp2.xtc -ox -com -nopbc

  99550 4.8424  5.00025 7.49519
  99600 4.79691 4.94304 7.39868
  99650 4.75789 4.85189 7.41581
  99700 4.74974 4.80423 7.48257
  99750 5.07858 5.04788 7.51654
  99800 4.99091 4.78976 7.47778
  99850 5.11406 4.84656 7.48769
  99900 5.09739 5.12395 7.47889
  99950 5.0411  4.88659 7.52737
  104.91271 5.07101 7.50397

For Z, the value is 7.48037+/-0.072562 (MIN=7.3 and MAX=7.8)

If I look at the frames from the min and max z values, they do clearly
differ in the placement of the bilayer along z.

###

To simplify things, I redid this while only using the phosphorous
atoms (1 per lipid) in the centering and coordinate output.

Here, Z=7.49323+/-0.0658727 (MIN=7.3 and MAX=7.7)

To verify this and since all the atoms have the same mass, I checked
it with awk:

echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 47550 -b 4 -o   
look_low.gro

echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 80750 -b 8 -o
look_high.gro
cat look_low.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.39197
cat look_high.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.62207

###

Finally, I rechecked with a single atom and here it worked exactly as
expected. The atoms is placed at 5 5 7.5 in every single frame.

###

Then I checked using 2 phosphorous atoms from the same (or,
separately, different) leaflets, and the result was the same:

$ tail coord.xvg
  99550 5   5   7.5
  99600 5   5.0005  7.5
  99650 5.0005  5.0005  7.5005
  99700 5   5.0005  7.5
  99750 5.0005  5   7.5005
  99800 5   5.0005  7.5005
  99850 5   5   7.5005
  99900 5.0005  5   7.5005
  99950 5.0005  5   7.4995
  105   5   7.5


Then using 3 phosphorous atoms from (containing atoms from both
leaflets), and now there is more deviation:

$ tail coord.xvg
  99550 5.48833 5.30233 6.972
  99600 5.48733 5.29233 6.97167
  99650 5.50033 5.29567 6.93467
  99700 5.52133 5.27433 6.93967
  99750 5.489   5.253   6.99133
  99800 5.481   5.319   6.91533
  99850 5.47867 5.309   6.827
  99900 5.48733 5.292   6.83133
  99950 5.52767 5.276.86867
  105.518   5.319   6.89167

(MIN=6.7 and MAX=7.1)

I originally had problems with pbc during trjconv in which the -center
option would place the bilayer at the peripheries of the unit cell in
some frames (COM should still be at the center). This is what lead me
to the route that I have used above (a .tpr that was generated with
pbc=no and trjconv -pbc none -center). Nevertheless, I think that
perhaps trjconv is using pbc information even in this case where ti
should not.

##

If anybody can see what I am doing wrong or has a method to center a
bilayer at a cartesian coordinate, then I am very interested.

Thank you for your time,
Chris.







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

[gmx-users] g_density

[chris . neale]

Dear Matthias:

Did you run trjconv -center -pbc mol at some point before running  
g_density? With pbc, there are multiple ways to put the center of a  
system at the center of the box. This can lead to having the water  
slab in the center of the unit cell and the bilayer on the top and  
bottom z.


Without seeing the commands and the output, it is hard to say any more.

Chris.

-- original message --

Hi,

I have aligned a trajectory to a reference structure and have now run
g_density on both the original and the aligned trajectory in order to
find the bilayer headgroup position and the bilayer thickness.

g_density works fine with the original trajectory but there seems to
be a bug when running it to the aligned structure (Alignment done
using trjconv -fit progressive).

The density is calculated along the z direction (-d Z). The density
for positive z is calculated correctly, but the density for negative z
is shifted about 15nm to to the top end of the box size, as if the
density distribution along the z-axis were discontinuous.

Other then stitching the data together by hands, are there any obvious
better solutions?

Thanks,

Matthias

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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_dist

[chris . neale]

Dear Nilesh:

You don't seem to have made any progress since Mark gave you some  
reading hints. Your best chance to get a useful response follows the  
general form:


1. I want to do A.
2. I tried method B, here are the exact commands that I used (copy and paste)
3. With that method, I obtain C
4. But I am confused about why C does not look like D.
5. What is the problem?

Basically, it helps to show that you are doing your own work too.

For example, try a command on a single frame and then look at that  
frame with VMD and count for yourself and verify the answer that way.


Chris.

-- original message --

I have a system with 128 emi (cations) and 128 Cl (anions).

I want to calculate how many CL atoms are in cutoff distance relative to
hydorgen aotm of cation.  I considered all CL atoms are distinguishable.

Basically I want to calcualte the distance between each CL atom and
correponing hydorgen and registered those CL atoms which are whthin
cutoff.

How can I do ?

I am using Gromacs 4.0.7 version.

Thanks

Nilesh










--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Can't unfold the protein

[chris . neale]

Dear Hsin-Lin:

The first thing that comes to mind is using very high temperatures with NVT.
Obtaining the desired denatured state is a complex challenge, but you  
might try this paper:


Phys. Rev. Lett. 93, 238105 (2004)
Reversible Temperature and Pressure Denaturation of a Protein  
Fragment: A Replica Exchange Molecular Dynamics Simulation Study


Still, your post seems to be about simply being unable to unfold at  
all. I have unfolded a very stable beta-barrel by either simulating at  
3000 K or using the pull code to pull the termini apart. This was not  
a real study -- I only used this to make a reverse folding movie for  
non-scientists -- but I can verify that it is possible to unfold even  
very stable proteins by these methods.


Chris.

-- original message --

Hi,

I have question about unfolding.
I use three different ways respectively but all failed.
The three way is,
1. 600K in 20ns, then the system explode.
2. 400K in 10ns, the protein is very stable and the value is almost  
the same in radius of gyration.
3. heating up from 300K to 400K in 2ns and the other 2ns for 400K  
only, then the protein is still stable and the value is almost the  
same in radius of gyration.

Below is the mdp file of the heating.
What's wrong with my system?



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Can't unfold the protein

[chris . neale]

Dear Hsin-Lin:

I am no expert in this area. I am just saying that if you get  
densities that are way too low in NPT, then you might alleviate this  
problem with NVT.


In fact, I would personally do this in vacuum without pbc and use the  
sd integrator. Then you can sample extended conformations. I'd  
consider this a method to generate unfolded conformations that  
probably have no relation to the true temperature denatured state. But  
let's be honest, you're never going to come close to Boltzmann  
sampling of an unfolded state in 300 ns so why bother with the water?


Your link ( http://manual.gromacs.org/online/protunf.html ) is only  
concerned with analysis, so it is irrelevant.


If you are determined to use NPT then I suppose that you could use the  
Berendsen barostat (more stable in simulations but you get the wrong  
ensemble) with a very low compressibility. That's about all I can say,  
perhaps somebody else will comment.


Good luck,

Chris.

-- original message --

Dear Chris,

Thank you for your reply.
My protein is very stable.

I simulated it in 300ns in 300K before but there was almost no change.
That's why I want to do denaturation now.

I want to start a new simulation on the protein which is unfolded.
And I'll read the paper you told me.
Thank you.

But I don't understand what you recommend me to do.
You use 3000K on your protein and the protein unfolded.
I use 600K on my protein but the system explode(box become very big  
and protein locate at corner).


I saw the second way on this web,
http://manual.gromacs.org/online/protunf.html
But it also useless to me.

And in my mdp I use thermostat and barostat, which means my system is NPT.
Does anything wrong in my methods?

Sincerely yours,

Hsin-Lin



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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_hbond

[chris . neale]

From g_hbond -h

"OH and NH groups are regarded as donors, O is an acceptor always, N  
is an acceptor by default, but this can be switched using -nitacc.  
Dummy hydrogen atoms are assumed to be connected to the first  
preceding non-hydrogen atom."


The two groups are to find hbonds between group A and group B. If you  
want to find hbonds within group A, then you specify an index group  
that contains group A two times.


I think that you want:

[ group ]
100 101 200

Then give group 2x as the identical index group.

PS, there were problems with g_hbond a few months ago reported on  
list, I suggest that you check to ensure that these problems were fixed.


Chris.

-- original message --

Hi,
To help clarify my previous g_hbond question, I have these two groups  
in my index file

[ a_100_a_101]  # 100 is N, #101 is H
   100 101
[ a_200 ]# 200 is O
   200

when I ran g_hbond, I got "Found 1 donors and 2 acceptors". The part  
that I don't quite understand is why it says 2 acceptors. Did I set up  
the groups correctly?


Again, hope this clarifies my question.
Thanks for your insight.

Simon



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] gmx4.5.4 genion problem: No line with moleculetype 'SOL' found

[chris . neale]

For the quick fix:

1. run genion on your topology that does work. Look at this to see the  
format of the ion atom and residue names


2. Pick a few waters in the structure containing the ligand and  
replace the OW by the ion and remove the hydrogens, then fix the  
number of atoms on the second line of the file and run editconf on the  
file.


** But I wouldn't be satisfied with that if I were you. there is some  
problem here that you should uncover now before you move farther.


Show us your commands exactly (copy and paste) and perhaps we can help more.

As a side note: your topology that you presented  is going to cause  
you problems if you try to put position restraints on the water, as it  
is incorrectly ordered:


 ; Include water topology
 #include "gromos43a1.ff/spc.itp"

 ; Include ligand topoloty
 #include "drg.itp"
 #ifdef POSRES_WATER
 ; Position restraint for each water oxygen
 [ position_restraints ]
 ;  i funct   fcxfcyfcz
11   1000   1000   1000
 #endif


-- original message --

Thank you for your suggestions.
I have checked the topology file with vim, and it looks pefectly, also, the
only problem happens when I use genion. One thing that might be possible is
that I use the double precision version of gromacs, because when I solve it
in water, the written of topology file looks weired(like I showed).
What do you mean by work-around? You mean just manually change some water
molecules in topology file and coordinate file into ions? Will the program
recognize something like CL, NA? Or can I just add ions to the system and
solve it in water?

Thank you very much

Best
Wishes


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] is lincs used with virtual hydrogens?

[chris . neale]

Dear Users:

If I create the topology of a peptide like this:

pdb2gmx -f protein.gro -vsite hydrogens

And then simulate it in vacuum, is lincs used at all? I believe that  
it is, as if I use a timestep that is too large then I get LINCS  
warnings about angles rotating more than 30 degrees, but that warning  
message could possibly have been written with the assumption that I  
used LINCS and not virtual hydrogens.


Finally, is there a method that needs to be named or cited in relation  
to the fact that the angles are now constrained? Is that also done  
with P-LINCS?


Thank you,
Chris.


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

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

[gmx-users] is lincs used with virtual hydrogens?

[chris . neale]

Thank you Roland.

I did use:

constraints = all-bonds
lincs-iter =  1
lincs-order =  6
constraint_algorithm =  lincs

From looking at the manual, I figured that angle and bond constraints  
would all be done by LINCS if I had done (A):


pdb2gmx -vsite none
constraints = h-angles
(a combination that I have never tried)

But when I use (B):

pdb2gmx -vsite hydrogen
constraints = all-bonds

It seems possible to me that LINCS is not used but instead the  
position of the atom is simply built from a mathematical function.  
Perhaps this all stems from my lack of thorough understanding of  
LINCS, but it seems to me that there need be no iteration to simply  
place an atoms based on virtual_sites3 (which are constructed by  
pdb2gms -hydrogen)


For now, I'll simply add a line to state that I built virtual sites  
for hydrogen atoms to make it clear, but I'd still like to understand  
the difference between options A and B, above, if you have some time.


Thank you again,
Chris.





On Sat, Jun 18, 2011 at 4:13 PM, chris.neale at utoronto.ca <
chris.neale at utoronto.ca> wrote:


Dear Users:

If I create the topology of a peptide like this:

pdb2gmx -f protein.gro -vsite hydrogens

And then simulate it in vacuum, is lincs used at all? I believe that
it is, as if I use a timestep that is too large then I get LINCS
warnings about angles rotating more than 30 degrees, but that warning
message could possibly have been written with the assumption that I
used LINCS and not virtual hydrogens.


Probably. To make sure check the constraint-algorithm selected in your mdp.
BTW: If you want to use large timecheck you should normally use
constraints=all-bonds and lincs-order=6.




Finally, is there a method that needs to be named or cited in relation
to the fact that the angles are now constrained? Is that also done
with P-LINCS?


This is also done with P-LINCS. Not sure whether one should sign something
regarding the construction/usage of v-sites.

Roland




Thank you,
Chris.



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

[gmx-users] is lincs used with virtual hydrogens?

[chris . neale]

Thank you Mark, I really appreciate it.

-- original message --

On 20/06/2011 1:23 AM, chris.neale at utoronto.ca wrote:

Dear Mark:

Now I am confused. Your first post indicated that P-LINCS did the  
angle constraints. But here you indicate that the v-site algorithm  
does it.


No... my second email observed that LINCS warnings can result when
v-sites are being used, and I ascribed that effect to the size of the
time step. I didn't say v-sites were doing anything with angle constraints.

This is probably because my first post was incomplete about the  
method that I used.


Can you please confirm my current understanding?

Therefore (A):

 pdb2gmx -vsite none
 constraints = h-angles
 constraint_algorithm =  lincs

--> P-LINCS constrains all bonds and angles involving hydrogens

But when I use (B):

 pdb2gmx -vsite hydrogen
 constraints = all-bonds
 constraint_algorithm =  lincs

--> Hydrogens are built as V-sites and P-LINCS constrains all bonds  
that do not involve hydrogen.


Of course, water is treated by SETTLE in both of the above.

Please note that I am not concerned about getting LINCS warnings. I  
probably should not have mentioned that. I am simply trying to  
figure out exactly how to write my methods section.


Both the above are correct understandings.

I think you misread my second email. I'd intended to put my
parenthetical observation about LINCS in the first email, but forgot to
(however you seem to have understood my context fine). Sorry for
whatever. :)

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] is lincs used with virtual hydrogens?

[chris . neale]

Dear Mark:

Now I am confused. Your first post indicated that P-LINCS did the  
angle constraints. But here you indicate that the v-site algorithm  
does it. This is probably because my first post was incomplete about  
the method that I used.


Can you please confirm my current understanding?

Therefore (A):

 pdb2gmx -vsite none
 constraints = h-angles
 constraint_algorithm =  lincs

--> P-LINCS constrains all bonds and angles involving hydrogens

But when I use (B):

 pdb2gmx -vsite hydrogen
 constraints = all-bonds
 constraint_algorithm =  lincs

--> Hydrogens are built as V-sites and P-LINCS constrains all bonds  
that do not involve hydrogen.


Of course, water is treated by SETTLE in both of the above.

Please note that I am not concerned about getting LINCS warnings. I  
probably should not have mentioned that. I am simply trying to figure  
out exactly how to write my methods section.


Thank you,
Chris.

-- original message --

Mark

On Sat, Jun 18, 2011 at 4:40 PM, chris.neale at utoronto.ca  
 > wrote:


Thank you Roland.

I did use:

constraints = all-bonds
lincs-iter =  1
lincs-order =  6
constraint_algorithm =  lincs

 From looking at the manual, I figured that angle and bond constraints
would all be done by LINCS if I had done (A):

pdb2gmx -vsite none
constraints = h-angles
(a combination that I have never tried)

But when I use (B):

pdb2gmx -vsite hydrogen
constraints = all-bonds

It seems possible to me that LINCS is not used but instead the
position of the atom is simply built from a mathematical function.
Perhaps this all stems from my lack of thorough understanding of
LINCS, but it seems to me that there need be no iteration to simply
place an atoms based on virtual_sites3 (which are constructed by
pdb2gms -hydrogen)

For now, I'll simply add a line to state that I built virtual sites
for hydrogen atoms to make it clear, but I'd still like to understand
the difference between options A and B, above, if you have some time.

Thank you again,
Chris.





On Sat, Jun 18, 2011 at 4:13 PM, chris.neale at utoronto.ca
 <
chris.neale at utoronto.ca > wrote:

> Dear Users:
>
> If I create the topology of a peptide like this:
>
> pdb2gmx -f protein.gro -vsite hydrogens
>
> And then simulate it in vacuum, is lincs used at all? I believe that
> it is, as if I use a timestep that is too large then I get LINCS
> warnings about angles rotating more than 30 degrees, but that
warning
> message could possibly have been written with the assumption that I
> used LINCS and not virtual hydrogens.
>
Probably. To make sure check the constraint-algorithm selected in
your mdp.
BTW: If you want to use large timecheck you should normally use
constraints=all-bonds and lincs-order=6.


>
> Finally, is there a method that needs to be named or cited in
relation
> to the fact that the angles are now constrained? Is that also done
> with P-LINCS?
>
This is also done with P-LINCS. Not sure whether one should sign
something
regarding the construction/usage of v-sites.

Roland


>
> Thank you,
> Chris.





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

[gmx-users] Programs to add residues

[Chris Neale]

Try "Loopy". You can get it to build termini in addition to loops.

http://wiki.c2b2.columbia.edu/honiglab_public/index.php/Software:Loopy

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
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] restraints in PMF (Justin's tutorial)

[chris . neale]

Actually, this depends on what you are trying to achieve. If you  
simply want to obtain the standard binding free energy, and somehow  
you know that the bound state is represented by your umbrella at  
dist=0 (via a crystal structure, for example), then using additional  
restraints is common, acceptable, and in many cases desirable.


You can do this in analogy to, for example: D. L. Mobley, J. D.  
Chodera, K. A. Dill. "On the use of orientational restraints and  
symmetry number corrections in alchemical free energy calculations",  
Journal of Chemical Physics 125:084902, 2006. Selected for the Virtual  
Journal of Biological Physics Research 12(5), 2006.


Since the free energy is a state function, you are free to select any  
pathway and it is useful to pick one that converges quickly... and  
those involving rotation of large molecules are sure to converge very  
slowly.


If you do this, then you'll need to account for the restraints in your  
free energy term in a way similar to that outlined in the paper that I  
referenced above.


The only reason that you would want to stick to the simple PMF is if  
you are actually interested in the shape of the unrestrained PMF, in  
which case you can't add additional restraints.


Chris.

-- original message --

Rebeca García Fandiño wrote:

Thanks a lot for your quick answer!
I think they are separated enough, however my monomers are cyclic  
(like discs); I start with a parallel conformation between then, but  
along the Umbrella simulation, both of them rotate and approach.

If I do not use restraints, how could I avoid the rotation?



You don't.  Why wouldn't two molecules rotate freely in solution when  
binding or

unbinding?  It seems like completely natural behavior.  Even in simple systems
of protein-ligand association, part of the binding energy is the entropic
restriction of the ligand into a certain binding-competent pose.  Why wouldn't
that happen here?  Sounds like an artificial notion to me.

-Justin


I am using the following md_umbrella.mdp:

title   = Umbrella pulling simulation
define  =
define  =
; Run parameters
integrator  = md
dt  = 0.002
tinit   = 0
nsteps  = 50   ; 1 ns
nstcomm = 10
; Output parameters
nstxout = 5000nstvout = 5000
nstfout = 5000
nstxtcout   = 5000nstenergy   = 5000
; Bond parameters
constraint_algorithm= lincs
constraints = all-bonds
continuation= yes
; Single-range cutoff scheme
nstlist = 5
ns_type = grid
rlist   = 1.4
rcoulomb= 1.4
rvdw= 1.4
; PME electrostatics parameters
coulombtype = PME
fourierspacing  = 0.12
fourier_nx  = 0
fourier_ny  = 0
fourier_nz  = 0
pme_order   = 4
ewald_rtol  = 1e-5
optimize_fft= yes
; Berendsen temperature coupling is on in two groups
Tcoupl  = Nose-Hoover
tc_grps = r_1_r_2  CL3
tau_t   = 0.5   0.5
ref_t   = 300   300
; Pressure coupling is on
Pcoupl  = Parrinello-Rahman
pcoupltype  = isotropic
tau_p   = 1.0
compressibility = 4.5e-5
ref_p   = 1.0
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr= EnerPres
; Pull code
pull= umbrella
pull_geometry   = distance
pull_dim= N N Y
pull_start  = yes
pull_ngroups= 1
pull_group0 = r_1
pull_group1 = r_2
pull_init1  = 0
pull_rate1  = 0.0
pull_k1 = 1000  ; kJ mol^-1 nm^-2
pull_nstxout= 1000  ; every 2 ps
pull_nstfout= 1000  ; every 2 ps

Thanks a lot again for your help.

Best wishes,

Rebeca.


 > Date: Wed, 22 Jun 2011 10:53:16 -0400
 > From: jalemkul at vt.edu
 > To: gmx-users at gromacs.org
 > Subject: Re: [gmx-users] restraints in PMF (Justin's tutorial)
 >
 >
 >
 > Rebeca García Fandiño wrote:
 > > Hello,
 > > I am trying to obtain the PMF from Umbrella Sampling of the process of
 > > separating two monomers of a dimer, following Justin 's tutorial
 > >  
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html

 > >
 > > I have done the Umbrella Sampling simulations without using any
 > > restraints in any of the monomers, however, I can see that they move and
 > > gyrate so that although the c.o.m are separated from each other, there
 > > are parts of both interacting, in such way that they are not separated
 > > as they should be.
 > >
 > > Would it be correct if I apply restraints to both monomers in all the
 > > Umbrella Sampling simulations?. I have seen that in Justin's tutorial
 > > they applied restraints to one of the chains, but in my case I think I
 > > will need to restrain both of the units. Would that be correct for the
 > > PMF calculation?
 > >
 >
 > There are no position restraints applied during the umbrella sampling
 > simulations. They are unnecessary. The umbrella potential is itself a
 > restraint to maintain COM se

[gmx-users] Timing variability

[chris . neale]

Dear Users:

Has anybody else looked at simulation speed (ns/day) over the segments  
of long runs? I always benchmark and optimize my systems carefully,  
but it was only recently that I realized how much variability I am  
obtaining over long runs. Perhaps this is specific to my cluster,  
which is one reason for my post.


Here is the performance (in ns/day) that I obtain with a single run.  
This is representative of what I see with 30 other runs (see list  
below). First, the distribution is bimodal, with peaks at around 80  
and 90 ns/day. Second, the values go as low as 70 ns/day (and I have  
seen as low as 50 ns/day when I look through all 30 run directories  
that differ only by the position of the umbrella restraint).


I am using gromacs 4.5.3 and I run it like this:
mpirun mdrun_mpi -deffnm md1 -dlb yes -npme 16 -cpt 30 -maxh 47 -cpi  
md1.cpt -cpo md1.cpt -px coord.xvg -pf force.xvg -noappend


I have obtained similar behaviour also with another system.

I have attempted to correlate timing with the node on which the job  
runs, the output of the "^Grid:" in the .log file, the DD division of  
PME and real-space and the value of pme mesh/force, all to no avail. I  
have found one correlation whereby the very slow runs also indicate a  
high relative load for PME.


I suspect that it is some value that is being determined at run start  
time that is affecting my performance, but I am not sure what this  
could be. Perhaps the fourier-spacing is leading to different initial  
grids?


Thank you (timing information and my .mdp file follows),
Chris

### Output the ns/day obtained by the run segments
$ grep Perf *log|awk '{print $4}'

78.286
82.345
81.573
83.418
92.423
90.863
85.833
91.131
91.820
71.246
76.844
91.805
92.037
85.702
92.251
89.706
88.590
89.381
90.446
81.142
76.365
76.968
76.037
79.286
79.895
79.047
78.273
79.406
78.018
78.645
78.172
80.255
81.032
81.047
77.414
78.414
80.167
79.278
80.892
82.796
81.300
77.392
71.350
73.134
76.519
75.879
80.684
81.076
87.821
90.064
88.409
80.803
88.435

### My .mdp file

constraints = all-bonds
lincs-iter =  1
lincs-order =  6
constraint_algorithm =  lincs
integrator = sd
dt = 0.004
tinit = 100
init_step=  0
nsteps = 10
nstcomm = 1
nstxout = 10
nstvout = 10
nstfout = 10
nstxtcout = 25000
nstenergy = 25000
nstlist = 5
nstlog=0 ; reduce log file size
ns_type = grid
rlist = 1
rcoulomb = 1
rvdw = 1
coulombtype = PME
ewald-rtol = 1e-5
optimize_fft = yes
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
tc_grps =  System
tau_t   =  1.0
ld_seed =  -1
ref_t = 300
gen_temp = 300
gen_vel = yes
unconstrained_start = no
gen_seed = -1
Pcoupl = berendsen
pcoupltype = semiisotropic
tau_p = 4 4
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0

; COM PULLING
pull = umbrella
pull_geometry= position
pull_dim = N N Y
pull_start   = no
pull_nstxout = 250
pull_nstfout = 250
pull_ngroups = 1
pull_group0  = POPC
pull_pbcatom0= 338
pull_group1  = KSC
pull_pbcatom1= 0
pull_init1   = 0 0 0.0
pull_rate1   = 0
pull_k1  = 3000.0
pull_vec1= 0 0 0
;;;EOF


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] T->A mutation for a dimer protein

[chris . neale]

Dear Lishan:

First, it would be great to see some evidence that you have tried to  
do this yourself before posting. Your "I think" makes it a possible  
waste of time for us to suggest a resolution to a problem that may or  
may not exist.


Second, if indeed it is a problem, perhaps you could use two identical  
topologies with different molecule names and treat them separately.


Chris

-- original message --

Hi All,
I have a protein dimer and I want to calculate a T to A mutation free  
energy change using TI method. Since it is a dimer, it is very  
convenient (and advantageous) to mutate the T in both monomer  
simultaneously. Gromacs will write out dH/dl sum for the two mutations  
together, I think. My question is how can I change Gromacs so than  
dH/dl for T->A mutation in each monomer is written out separately?


Best,
Lishan


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] pulling code

[chris . neale]

Dear Adam:

I have a number of questions listed below. But first, I think it would  
be useful if you explain exactly what you want to accomplish and how  
you tried to do this with non-modified code and what problem you ran  
into that caused you to modify the code. After you outline that,  
please address the following questions:


1. Why not use one terminus as the reference group and the other  
terminus as the pull_group0 ?


2. How did you modify the code to use the shortest vector? I was under  
the impression that the code did already do this (shortest vector  
accounting for PBC) and that this was actually something you wanted to  
avoid.


3. Why can you not just use A-B even if A is positive and B is  
negative? If you are trying to avoid determining a PBC-wraparound  
distance, then it seems like A-B does give you exactly what you want  
as long as the entire protein is inside the unit cell.


4. If you really need to know the box size, t is available in the code  
at runtime. Why can you not just use the appropriate variable?


5. How are you fixing the non-terminal groups that you say need to be  
fixed? I don't understand why you can't fix groups with position  
restraints and the pull with the pull code. Also, why do you need to  
fix certain non-terminal groups?


6. This confused me a lot: "I don't think I want to freeze or restrain  
the groups that need to be fixed, because they must still be allowed  
flexibility to twist and deform--I only want to roughly fix their  
centers of mass." ... so then how did you fix them?


Chris.

-- original message --

Hi all,
I am trying to measure the energy curve of unwinding the SNARE protein
complex.  This requires pulling the C-termini of two helices apart from each
other, while fixing certain other parts of the complex.  To do this I
excluded the reference group, so that all groups were pulled to absolute
coordinates.  However, as the protein was centered in the box, and because
the box center is where the coordinates transition from -boxlength/2 to
+boxlength/2 when u7. Why do you need to fix certain non-sing pbc (ie.  
the dateline), every time a group passes

through the center of the box in any of the 3 dimensions, it is suddenly
yanked violently through the box because it thinks it is now a full
boxlength away from its destination.  I modified the pull code to take the
shortest vector between the current position and destination, which seems to
fix this problem.  I also modified g_wham so that it can take the difference
between two groups' positions as the energy curve coordinate (since several
groups need to be fixed, neither of these two groups can be used as the
reference group).  However, as I use pressure coupling, the box dimensions
change during the simulation, and if one group has a positive and the other
a negative coordinate because they straddle the center of the box, the box
dimensions must be known at each timestep to get their actual separation.
 Is there an easier way to handle this type of simulation?  I don't think I
want to freeze or restrain the groups that need to be fixed, because they
must still be allowed flexibility to twist and deform--I only want to
roughly fix their centers of mass.  Any suggestions would be appreciated,

Thanks,

Adam


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] gromacs 2.1

[chris . neale]

Thank you Rossen and Roland.

I was able to obtain the source but not compile it. It would be  
appreciated if you can offer any suggestions.


I obtained gromacs via:

git clone git://git.gromacs.org/gromacs.git
git checkout release-2-0-01

But I can not compile it. When I cd src and run make or gmake I get  
the errors listed below.


...  ...
gmake[348]: Entering directory  
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd  
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;  
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;

/bin/sh: line 0: cd: include: No such file or directory
gmake[349]: Entering directory  
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd  
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;  
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;

...  ...

I tried the documentation, but the FAQ file at  
file:///home/cneale/work/June2011/gromacs/share/html/gmxfaq.html#install links  
to a non-existant file  
file:///home/cneale/work/June2011/gromacs/share/html/readme2.0.html#install  
.


I moved the install directory into src/ and that solved the error  
messages that I was obtaining, but now make errors with:


/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/include/smalloc.h:36:2: warning: #ident is a deprecated GCC  
extension
dlist.c:197:1: error: pasting "." and "H" does not give a valid  
preprocessing token
dlist.c:197:1: error: pasting "." and "N" does not give a valid  
preprocessing token
dlist.c:197:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token
dlist.c:197:1: error: pasting "." and "C" does not give a valid  
preprocessing token
dlist.c:200:1: error: pasting "." and "N" does not give a valid  
preprocessing token
dlist.c:200:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token
dlist.c:200:1: error: pasting "." and "C" does not give a valid  
preprocessing token
dlist.c:200:1: error: pasting "." and "O" does not give a valid  
preprocessing token
dlist.c:203:1: error: pasting "." and "minO" does not give a valid  
preprocessing token
dlist.c:203:1: error: pasting "." and "minC" does not give a valid  
preprocessing token
dlist.c:203:1: error: pasting "." and "N" does not give a valid  
preprocessing token
dlist.c:203:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid  
preprocessing token

make[1]: *** [dlist.o] Error 1
make[1]: Leaving directory  
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/tools'


Also moving the lib directory into the src directory did not help any further.

The compiler was gcc version 4.4.0

Thank you very much,
Chris.

-- original message --

You can try also

$ git tag
release-1-6
release-1-6-01
release-2-0
release-2-0-01
release-2-0-02
release-3-0
release-3-1
release-3-1-1
release-3-1-2
release-3-1-3
release-3-1-4
release-3-2
v4.0.7
v4.5
v4.5-beta1
v4.5-beta2
v4.5-beta3
v4.5-beta4
v4.5.1
v4.5.2
v4.5.3
v4.5.4

# git checkout release-2-0-01


On 6/12/11 11:49 AM, Roland Schulz wrote:

Hi,

you can get it from git. "git log {file}" shows the history and git  
show {version}:{file} gives you a specific version.

In this case you can use:
git show 74b20ce9:src/tools/g_order.c



Is that the "official" 2.1? We should tag it then.

Rossen


Roland

On Sat, Jun 11, 2011 at 10:11 PM, > wrote:


Dear users:

Does anybody have the source code for gromacs version 2.1? I would
like to check the original source of g_order, but only versions
3.0 and above are available on the gromacs website.

Thank you,
Chris.

-- gmx-users mailing list gmx-users at gromacs.org

http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search 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-request at 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_Li

[gmx-users] gromacs 2.1

[chris . neale]

In the middle of my last post, it should have stated:

"I moved the **include** directory into src/ and that solved the error
messages that I was obtaining, but now make errors with:"

Sorry for the confusion,
Chris.

Quoting chris.ne...@utoronto.ca:


Thank you Rossen and Roland.

I was able to obtain the source but not compile it. It would be
appreciated if you can offer any suggestions.

I obtained gromacs via:

git clone git://git.gromacs.org/gromacs.git
git checkout release-2-0-01

But I can not compile it. When I cd src and run make or gmake I get
the errors listed below.

...  ...
gmake[348]: Entering directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;
/bin/sh: line 0: cd: include: No such file or directory
gmake[349]: Entering directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;
...  ...

I tried the documentation, but the FAQ file at
file:///home/cneale/work/June2011/gromacs/share/html/gmxfaq.html#install   
links

to a non-existant file
file:///home/cneale/work/June2011/gromacs/share/html/readme2.0.html#install
.

I moved the install directory into src/ and that solved the error
messages that I was obtaining, but now make errors with:

/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/include/smalloc.h:36:2: warning: #ident is a deprecated   
GCC

extension
dlist.c:197:1: error: pasting "." and "H" does not give a valid
preprocessing token
dlist.c:197:1: error: pasting "." and "N" does not give a valid
preprocessing token
dlist.c:197:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
dlist.c:197:1: error: pasting "." and "C" does not give a valid
preprocessing token
dlist.c:200:1: error: pasting "." and "N" does not give a valid
preprocessing token
dlist.c:200:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
dlist.c:200:1: error: pasting "." and "C" does not give a valid
preprocessing token
dlist.c:200:1: error: pasting "." and "O" does not give a valid
preprocessing token
dlist.c:203:1: error: pasting "." and "minO" does not give a valid
preprocessing token
dlist.c:203:1: error: pasting "." and "minC" does not give a valid
preprocessing token
dlist.c:203:1: error: pasting "." and "N" does not give a valid
preprocessing token
dlist.c:203:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
dlist.c:212:1: error: pasting "." and "Cn" does not give a valid
preprocessing token
make[1]: *** [dlist.o] Error 1
make[1]: Leaving directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/tools'

Also moving the lib directory into the src directory did not help   
any further.


The compiler was gcc version 4.4.0

Thank you very much,
Chris.

-- original message --

You can try also

$ git tag
release-1-6
release-1-6-01
release-2-0
release-2-0-01
release-2-0-02
release-3-0
release-3-1
release-3-1-1
release-3-1-2
release-3-1-3
release-3-1-4
release-3-2
v4.0.7
v4.5
v4.5-beta1
v4.5-beta2
v4.5-beta3
v4.5-beta4
v4.5.1
v4.5.2
v4.5.3
v4.5.4

# git checkout release-2-0-01


On 6/12/11 11:49 AM, Roland Schulz wrote:

Hi,

you can get it from git. "git log {file}" shows the history and git
show {version}:{file} gives you a specific version.
In this case you can use:
git show 74b20ce9:src/tools/g_order.c



Is that the "official" 2.1? We should tag it then.

Rossen


Roland

On Sat, Jun 11, 2011 at 10:11 PM, mailto:chris.neale at utoronto.ca>> wrote:

Dear users:

Does anybody have the source code for gromacs version 2.1? I would
like to check the original source of g_order, but only versions
3.0 and above are available on the gromacs website.

Thank you,
Chris.

-- gmx-users mailing list gmx-users at gromacs.org

http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search 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-request at 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 

[gmx-users] g_msd bug

[chris . neale]

600 GB of memory? I highly doubt that you have that much memory  
available. Are you sure that this is not a typo? Can you please post  
evidence that you have >=600 GB of memory available? It is common for  
clusters to disallow an individual process from using >10% of the  
total memory on a head-node, which makes 600 MB more likely, in which  
case you can try submitting your job to a compute node.


-- original message --

Hello,
thank you for your reply. I have used following command :

g_msd  -n POPC.ndx  -lateral z -o POPC_msd.xvg -mol POPC_diff.xvg

Trajectory has 1 frames and the system it was ran on is Fedora Red  
Hat 5.4.

Indeed my network administrator was very unhappy about comsumed memory.

Regards,
   Slawomir





Wiadomo¶æ napisana przez Tsjerk Wassenaar w dniu 2011-06-27, o godz. 15:40:

[Hide Quoted Text]
Hi Slawomir,

That's quite a usage of memory! Can you provide more information? Like
the number of frames in the trajectory, the command line you used, and
the system you ran on?

Cheers,

Tsjerk

2011/6/27 S³awomir Stachura :
Hi GMX Users,
I am writting this email, beacause I think the g_msd program in  
Gromacs 4.5.4 bears a problem. I was calculating the MSD od center of  
mass of POPC in membrane (system contains 274 POPC lipid molecules in  
all-atom force field) from 50 ns trajectory and it seems to consume  
great amount of memory. With  time of calculations the memory reserves  
are gradually devoured to the extent, in my case,  of over 600 GB  
(than my administrator of cluster killed the process). It seems that  
it does not release memory and it's pilling results up with steps  in  
memory. Have you heard of such case?

Best wishes,

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Announcement: Large biomolecule benchmark report

[chris . neale]

You should absolutely publish this. it would be of great interest. You  
can mitigate your chances of running into problems with the overview  
by sending a version of the manuscript to the developers of each  
software and asking them to provide a short paragraph, each of which  
you could include in a final section of responses from the developers.


A manuscript such as this (and indeed the information that you have  
already made available) will be very useful for many reasons. One  
reason is that when new PhD students start learning about simulations,  
they tend to use the package that has been adopted by their research  
group and trust the (probably biased and partially uninformed)  
statements of their senior colleagues.


Chris.

-- original message --

Thanks a a lot to you and also to Szilárd for your feedback and
encouragement.  I am very happy to see that this work is indeed useful
especially to developers.

We have no plans to make this into a 'proper' publication.  I am not
sure how much interested the simulation community would be because, to
be honest, I have no overview what has been done in this area (besides
the few benchmarks studies I have cited).

Thanks again,
Hannes.


On Thu, 15 Mar 2012 22:02:21 +0100
David van der Spoel  wrote:

[Hide Quoted Text]
On 2012-03-15 14:37, Hannes Loeffler wrote:
Dear all,

we proudly announce our third benchmarking report on (large)
biomolecular systems carried out on various HPC platforms.  We have
expanded our repertoire to five MD codes (AMBER, CHARMM, GROMACS,
LAMMPS and NAMD) and to five protein and protein-membrane systems
ranging from 20 thousand to 3 million atoms.

Please find the report on
http://www.stfc.ac.uk/CSE/randd/cbg/Benchmark/25241.aspx
where we also offer the raw runtime data.  We also plan to release
the complete joint benchmark suite at a later date (as soon as we
have access to a web server with sufficient storage space).

We are open to any questions or comments related to our reports.
It looks very interesting, and having benchmarks done by independent  
researchers is the best way to avoid bias. The differences are quite

revealing, but it's also good that you point to problems compiling
gromacs. Is this going to be submitted for publication somewhere too?

Thanks for doing this, it must have been quite a job!

Kind regards,
Hannes Loeffler
STFC Daresbury
--
Scanned by iCritical.

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] crazy temperatures

[chris . neale]

I disagree.

What one is generally trying to obtain with elevated temperatures is  
enhanced sampling, not temperature-dependent properties. I believe  
that even TIP4P-EW is not very good at getting the properties of water  
correct at 600 K, temperatures that are commonly used during replica  
exchange simulations (not to mention that nobody has any idea how  
accurate protein forcefields are at temperatures other than the one at  
which they were parameterized).


So I think that doing simulations at massively elevated temperatures  
can possibly be useful.


That said, while doing simulated annealing, I have found previously  
using charmm that once you get to about 3,000 K you will get chiral  
inversions that can not resolve at lower temperature. This is because  
our improper dihedral terms only maintain the given chirality, rather  
than favouring one over the other.


To address your question directly, I believe that chiral inversions  
will be a big problem for you at 20,000 K. Obviously you also have  
simulation stability issues, but one presumes that you could resolve  
those by using a small enough timestep.


Chris.

-- original message --

At that temperature most matter is going to be a plasma, not many  
bonds to be simulated and a lot of free electrons.


Warren Gallin

On 2012-03-28, at 4:43 PM, Mark Abraham wrote:

[Hide Quoted Text]
On 29/03/2012 9:39 AM, Asaf Farhi wrote:
Dear GMCS users

Hi. Does anyone know if MD at 2K is feasible?
Please start new email threads rather than hijacking old ones.

I doubt anybody knows the answer to your question. Force fields are  
parameterized to reproduce data at around 300K. I can't imagine any  
possible use for simulating an MM force field at a temperature hotter  
than the sun.


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] Umbrella sampling along Radius of gyration

[chris . neale]

It is not clear to me how one would do this with MD. The only thing  
that I can think of doing in gromacs is to create a virtual particle  
that is placed at the center of the protein and then apply forces  
along the vector from this virtual atom to each of the Ca atoms. You  
would need to modify the gromacs source code so that the Rg was  
calculated at each step and the magnitude of each force is scaled  
appropriately such that you get a harmonic potential about the desired  
value of Rg. It should be easier to do with MC (although that's a ways  
off for gromacs unless I've missed something).


Some people appear to have done exactly what you want with MD. I  
presume that they used charmm.


http://www.pnas.org/content/94/19/10161.long (and their reference 22).

Chris.

-- original message --

Hi,
  Is there any way in gromacs to use radius of gyration of a polymer  
as reaction coordinate for umbrella sampling ?

Thanks
Sanku



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

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

[gmx-users] Not able to continue with Equilibration

[chris . neale]

I didn't follow this whole thread, but I sometimes need to turn off  
all constraints when doing minimization. In fact, for that reason I  
entirely stopped ever using restraints during energy minimization. In  
extreem cases, I have had success also by forcing the water to be  
flexible with a -DFLEXIBLE (or whatever is appropriate for your water  
model; check the .itp) statement in the .mdp file.


Once EM is done, use rigid water and restraints and everything works out ok.

Chris.

-- original message --

Hallo Felix,
thank you for your answer. I tried the constraints = h-bonds but no
change in the output. If I look at the step.pdb file that is produced
after the running I have some strange outcome. For example some of my
atoms are not recognized as part of my protein any more and my
structure is destroied. Am I using some strange parameters for nvt
without realizing it? I've started from the mdp file provided by the
lysozyme tutorial non the Gromacs webpage.
if anyone has any ideas it is more than welcome.
Francesca

integrator  = md; leap-frog integrator
nsteps  = 1 ; 2 * 1 = 20 ps
dt  = 0.002 ; 2 fs
; Output control
nstxout = 1 ; save coordinates every 0.2 ps
nstvout = 1 ; save velocities every 0.2 ps
nstenergy   = 50; save energies every 0.2 ps
nstlog  = 50; update log file every 0.2 ps
; Bond parameters
unconstrained_start = no
;continuation   = no; first dynamics run
constraint_algorithm = lincs; holonomic constraints
constraints = h-bonds   ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter  = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid  ; search neighboring grid cells
nstlist = 5 ; 10 fs
rlist   = 3.0   ; short-range neighborlist cutoff (in nm)
rcoulomb= 3.0   ; short-range electrostatic cutoff (in nm)
rvdw= 3.0   ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME   ; Particle Mesh Ewald for long-range
electrostatics
pme_order   = 4 ; cubic interpolation
fourierspacing  = 0.16  ; grid spacing for FFT
; Temperature coupling is on
tcoupl  = berendsen ; 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 off
pcoupl  = no; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz   ; 3-D PBC
; Dispersion correction
DispCorr= EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes   ; assign velocities from Maxwell distribution
gen_temp= 300   ; temperature for Maxwell distribution
gen_seed= -1; generate a random seed



2012/3/29 Rausch, Felix :

[Hide Quoted Text]
Hello again,

Well, I once had problems with simulations crashing randomly during  
production runs (sometimes after tens of nanoseconds) with the LINCS  
warnings you described. Switching LINCS from "all-bonds" to only  
"h-bonds" did the trick for me, although I never exactly figured out  
why.

Maybe it's worth a try for you, too?
Cheers,
Felix

-Ursprüngliche Nachricht-
Von: gmx-users-boun...@gromacs.org  
[mailto:gmx-users-boun...@gromacs.org] Im Auftrag von francesca vitalini

Gesendet: Donnerstag, 29. März 2012 15:15
An: Discussion list for GROMACS users
Betreff: Re: [gmx-users] Not able to continue with Equilibration

Hi!
I'm having a similar problem. I have a dimer solvated in a big box of  
water plus ions that I have managed to minimize correctly (see output  
of minimization at the end) but when I try to run NVT equilibration  
(see later) I get LINCS warnings(see below) refearred to atoms which  
are not in a cluster or in a strange position. I have added dihedral  
restraints on them but still the same type of error. I'm using GROMACS  
3.3.1. I have tried to switch to a newer version of GROMACS but still  
the same error.

Does anyone have any suggestions?
Thanks
Francesca


MINIMIZATION OUTPUT

Steepest Descents converged to Fmax < 1000 in 681 steps Potential  
Energy  = -1.9597512e+07

Maximum force =  2.4159846e+02 on atom 13087
Norm of force =  2.1520395e+04


MINIMIZATION.MDP

define = -DEFLEXIBLE ; flexible water
integrator  = steep ; Algorithm (steep = steepest descent
minimization)
emtol   = 1000.0 ; Stop minimization when the maximum
force < 1000.0 kJ/mol/nm
emstep  = 0.001  ; Energy step size
nsteps  = 5 ; Maximum number of (minimization)
steps to perform

; 

[gmx-users] How to compile the template.c when GMX is builded with cmake

[chris . neale]

Dear Users:

I was previously able to compile template.c after I compiled gromacs  
with autoconf. I was unable to compile templae.c, however, I used  
cmake to compile gromacs. This is from gromacs-4.5.4.


I tried "cmake ." in the template directory with no success:

cmake .gpc-f102n084-$ cmake .
-- The C compiler identification is Intel
-- The CXX compiler identification is Intel
-- Check for working C compiler:  
/home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icc
-- Check for working C compiler:  
/home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icc  
-- works

-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working CXX compiler:  
/home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icpc
-- Check for working CXX compiler:  
/home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icpc  
-- works

-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- checking for module 'libmd'
--   package 'libmd' not found
CMake Error at CMakeLists.txt:42 (message):
  libmd not found, source GMXRC.


-- Configuring incomplete, errors occurred!

###

For this test, I used the template.c that came with gromacs.

This issue has previously been posted to the list (see below) but I  
was unable to find any answers to that previous querry.


I also looked at the README that came in the template directory and  
find it to be apparently autoconf specific, with no directions for  
cmake.


Thank you,
Chris.

quoting http://www.mail-archive.com/gmx-users@gromacs.org/msg39347.html

How to compile the template.c when GMX is builded with cmake

Bert
Thu, 31 Mar 2011 10:15:22 -0700

Dear all,

When GMX is builded with cmake, how to compile the template.c? I used
to make the template.c by the command "make -f MakefileXXX", but it
can not work in version 4.5. Thanks for your suggestions.

Best regards,
Bert
--




--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] on the Umbrella Sampling on-line tutorial

[chris . neale]

Dear Acoot:

I'll reply to general topics, not to the tutorial in particular.

If the opening is not large enough to allow the peptide to exit, then  
the surrounding protein will need to change its conformation to permit  
unbinding. This can happen, but you need (a) to have sufficiently long  
sampling times to permit the required conformational change and (b)  
starting structures in which the peptide is near the umbrella center  
and do not crash during simulation.


The US method is always valid, but it is not always the best choice.  
You might try using the free energy code to implement a thermodynamic  
cycle. Nevertheless, one imagines that the protein does actually need  
to open up during peptide binding and so you will still need to sample  
protein opening/closing in order to obtain equilibrium (i.e. correct)  
binding free energies because you need to sample the unbound protein  
at equilibrium. That is to say that a thermodynamic cycle may appear  
converged when in fact it is not, because you have not converged the  
unbound state of the protein. In any event, be careful with your  
convergence analysis.


This would be a good system for which to attempt both US and  
double-decoupling approaches and use the results of each to ensure  
that you are not missing some important conformational states.


That said, I highly doubt that it is possible to converge the free  
energies of induced-fit peptide-protein binding with any available  
atomistic computational method using contemporary computational  
resources. The lower bound of required sampling times is certainly on  
the order of 10 us per umbrella and I bet that the actual value is a  
few orders of magnitude larger. I have less experience with the free  
energy code than I do with US, but I suspect that the required  
sampling times are also very long for a system like this.


Chris.

-- original message --

Dear All,

I planned to use the method introduced in the Umbrella Sampling  
on-line tutorial of Justin Lemkul  
(http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html).


But if a peptide is surrounded by a protein, which means the opening  
of the protein-complex is not large enough to allow the peptide to  
leave the protein without significantly breaking the conformation of  
the protein in the protein-peptide complex, is the Umbrella Sampling  
method still valid for the binding energy calculation?


Will you please also show me in which part of the tutorial the  
direction of pull-apart is defined? We should process it in a  
direction the peptide can leave the protein, not the direction protein  
will bind the peptide much strongly.


I am looking forward to getting a reply from you.

Cheers,

Acoot


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] further discussion on the mdrun -append function

[chris . neale]

Dear Acoot:

You should be able to answer this one yourself. Moreover, you are  
doing yourself a disservice by relying on the mailing list to do your  
work for you because you will eventually need to learn how to find  
answers to these things on your own.


Please remember the following:

1. use a title that matches your question, starting a new thread if  
you have a new question


2. it may sound harsh, but questions like this one in which you have  
obviously not even tried to find the answer yourself for more than a  
couple of minutes tend to annoy the very people that you may want help  
from later on with other issues.


3. whenever you post, please show how you have tried to solve the  
problem yourself. You may find that in the process of writing such a  
question you end up solving the problem yourself.


4. read the manual.

Chris.

-- original message --


Thanks for the detailed explaination.

Will you please explain the function of "-n index.ndx " in "grompp -f  
md.mdp -n index.ndx -p topol.top -c minimized.gro -o md-1ns.tpr" with  
some specific examples?


Cheers,

Acoot


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Different results from identical tpr after MD

[chris . neale]

Dear Acoot:

The idea of convergence is this: start a large number of simulations  
from different conformations, analyze some quantity over time in each  
simulation, and when the deviation of the average value of that  
quantity from each separate simulation is less than the time-variance  
within individual simulations, then you can imply that the simulations  
have converged -- that is, the results of independent simulations  
which started as different are now similar.


There is a huge body of work that uses a single simulations and  
evaluates its so-called convergence using some assumptions and special  
methods. That can also be very useful, but I find it informative to  
think of "convergence" in its standard non-scientific dictionary  
definition as the coming together of previously disparate things.


For simulations, my working definition is this: a set of simulations  
has converged the value of some variable when the simulations were  
initiated from sufficiently distinct conformational basins and then,  
over time, the ensemble distribution of the time-averages of the  
specified variable has a variance that is the same as the mean  
time-averaged variation within independent simulations. The weak point  
here is the part about "sufficiently distinct conformations", but I am  
not sure that this can be stated less vaguely in the general case.


Chris.

-- original message --

Hi Justin,

Can you give me your definition of converged MD and unconverged MD?

Cheers,

Acoot


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] large error bars in PMF

[chris . neale]

I don't see why the box-type makes any difference whatsoever. It is  
possible that if you use a rhombic dodecahedron, you may reduce the  
system size, thus simulate more ns/day, thus converge faster, but that  
should be the only effect. I would be interested to hear more from  
Justin about how the box-shape is expected to affect peptide  
rotation... perhaps I misunderstand this point.


I have a few other comments.

1. If you allow the peptide to rotate freely, then you do indeed need  
to converge all of their different rotational interactions. An  
alternative is to apply orientational restraints during the pulling  
(assuming that you know the bound state) and then to correct to an  
unrestrained state at large separations. You can see, for instance,   
D. L. Mobley, J. D. Chodera, K. A. Dill. "On the use of orientational  
restraints and symmetry number corrections in alchemical free energy  
calculations", ...


2. You are using "pull_dim = N N Y" which, if I haven't entirely  
forgotten how the pull-code works, means that the distance along Z is  
restrained but the distance along X and Y is free to change. What you  
end up with by sampling in this way is pretty strange and will require  
a really weird volumetric correction in the absence of infinite  
sampling time. You must decide to either: (i) pull to a spherical  
distance:


pull_dim = Y Y Y
pull_geometry = distance

or (ii) to pull along a defined vector

pull_dim = Y Y Y
pull_geometry = direction

What you have done:

pull_dim = N N Y
pull_geometry = distance

is only really useful when the system is isotropic along the XY plane  
(at least in the time averaged sense), such as for a lipid bilayer, or  
when the freedom in X and Y is very low, such as in a channel.


3. Finally, just because you sampled *more* doesn't mean that your  
values are converged. Look into block averaging and test to see if  
your binding free energy is drifting over time.


Good luck,
Chris.

-- original message --

Hi again,
I have one doub about the suggestion of using a dodecahedral box for  
my umbrella sampling to remove the problems I am having with the  
peptides rotating. I cannot see why a dodecaheral box is going to  
avoid this. Would it be better a truncated octahedron?

Thanks a lot for your help.
Best wishes,
Rebeca.

[Hide Quoted Text]
Date: Thu, 21 Jul 2011 15:16:52 -0400
From: jalem...@vt.edu
To: gmx-users@gromacs.org
Subject: Re: FW: [gmx-users] large error bars in PMF



Rebeca García Fandiño wrote:


I am trying to achieve the binding energy of the dimer composed by the
two small cyclic peptides, to compare it with experimental. What
advantages would I have using 3D PMF instead only 1D for this calculation?
Intuitively, two molecules diffuse through solution until they find  
one another,
which to me sounds a lot like a 3D path.  Further, using a  
dodecahedral box for

your umbrella sampling removes the problems you're having with the peptides
rotating.  It sounds like you're trying to pull in one direction along a
rectangular box, but the peptides are not playing nice.  I feel like this
discussion has come up at least once or twice before, though...

-Justin
Thanks a lot!
Rebeca.

  > Date: Thu, 21 Jul 2011 14:14:44 -0400
  > From: jalem...@vt.edu
  > To: gmx-users@gromacs.org
  > Subject: Re: [gmx-users] large error bars in PMF
  >
  >
  >
  > Rebeca García Fandiño wrote:
  > > Hi,
  > > thanks a lot for your quick answer.
  > > What I am trying to pull are two small peptides one from another (r_1
  > > and r_2).
  > > I did not understand very well your last suggestion: "...if you want
  > > reasonable error bars you will not lots of well-converged data".
  >
  > Oops, that should have read "you will *need* lots of well-converged
data."
  >
  > > Do you mean I will need also more windows besides extending the
simulations?
  >
  > I doubt you need more windows. Likely you just need more time in each.
  >
  > > I think the problem could be also that the peptides I am using
rotate in
  > > the box and they do not remain flat one respect to the other. They
  > > gyrate freely and some parts of their structure interact along the
  > > pulling...
  >
  > Interactions are part of the dissociation process and are not
problematic per
  > se. But if you're trying to obtain only a one-dimensional PMF then your
  > rotation could be a problem. Is there some reason you need a
one-dimensional
  > PMF and not a three-dimensional PMF? What are you trying to achieve?
  >
  > -Justin
  >
  > > Thanks a lot again for your help.
  > > Best wishes,
  > > Rebeca.
  > >
  > >
  > >
  > >

  > > From: rega...@hotmail.com
  > > To: gmx-users@gromacs.org
  > > Date: Thu, 21 Jul 2011 16:36:59 +
  > > Subject: [gmx-users] large error bars in PMF
  > >
  > >
  > > Hi,
  > > I am trying to calculate the binding energy of two molecules using the
  > > PMF (Umbrella Sampling method) and Groma

[gmx-users] large error bars in PMF

[chris . neale]

I see now what you mean. As it happens, I doubt that this would have  
caused the problem since no force was applied on X and Y dimensions,  
so it would require that there was a PBC-based distance degeneracy  
along Z, although this is of course possible and hopefully Rebecca  
will answer this part.


Also, thanks for the pull_dimension/pull_vec fix.

Chris.

-- original message --


chris.ne...@utoronto.ca wrote:

[Hide Quoted Text]
I don't see why the box-type makes any difference whatsoever. It is
possible that if you use a rhombic dodecahedron, you may reduce the
system size, thus simulate more ns/day, thus converge faster, but that
should be the only effect. I would be interested to hear more from
Justin about how the box-shape is expected to affect peptide rotation...
perhaps I misunderstand this point.
My point was not that the box shape has any effect on protein rotation.  That
will happen regardless of the box shape, of course.  My suggestion for  
this box

type was that since Rebeca has a system that is essentially spherically
symmetric (i.e. two proteins connected by some arbitrary vector, which are at
the same time rotating freely), then she must use a suitable box shape that
reflects this type of symmetry.  I never got a clear answer to whether or not
the weird interactions she cited were due to PBC or not, but if one uses a
rectangular box for a system like this one, there can be artificial  
interactions

very easily.

[Hide Quoted Text]
I have a few other comments.

1. If you allow the peptide to rotate freely, then you do indeed need to
converge all of their different rotational interactions. An alternative
is to apply orientational restraints during the pulling (assuming that
you know the bound state) and then to correct to an unrestrained state
at large separations. You can see, for instance,  D. L. Mobley, J. D.
Chodera, K. A. Dill. "On the use of orientational restraints and
symmetry number corrections in alchemical free energy calculations", ...

2. You are using "pull_dim = N N Y" which, if I haven't entirely
forgotten how the pull-code works, means that the distance along Z is
restrained but the distance along X and Y is free to change. What you
end up with by sampling in this way is pretty strange and will require a
really weird volumetric correction in the absence of infinite sampling
time. You must decide to either: (i) pull to a spherical distance:

pull_dim = Y Y Y
pull_geometry = distance

or (ii) to pull along a defined vector

pull_dim = Y Y Y
pull_geometry = direction
Just a note here - if you set direction geometry, the pull_dim keyword is not
used, but pull_vec is.

-Justin

[Hide Quoted Text]
What you have done:

pull_dim = N N Y
pull_geometry = distance

is only really useful when the system is isotropic along the XY plane
(at least in the time averaged sense), such as for a lipid bilayer, or
when the freedom in X and Y is very low, such as in a channel.

3. Finally, just because you sampled *more* doesn't mean that your
values are converged. Look into block averaging and test to see if your
binding free energy is drifting over time.

Good luck,
Chris.

-- original message --

Hi again,
I have one doub about the suggestion of using a dodecahedral box for my
umbrella sampling to remove the problems I am having with the peptides
rotating. I cannot see why a dodecaheral box is going to avoid this.
Would it be better a truncated octahedron?
Thanks a lot for your help.
Best wishes,
Rebeca.

<... snip...>


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] umbrella sampling

[chris . neale]

format your data for 2D-WHAM with 1D being the distance and the 2nd-D  
being your other coordinate of interest. Specify a value of zero for  
the force constants for your 2nd-D. Run 2D-WHAM. Boltzmann project the  
2D PMF onto your 2nd-D.


I think you can also do essentially the same thing by re-weighting  
using the F-values from 1D-WHAM, but I find the above method to be the  
simplest. It also provides you with a 2D free energy profile, which  
can be informative both biologically and to indicate on sampling  
problems.


Note that you're very likely going to run into convergence problems  
since your 2nd-D will rely on brute-force to converge, and worse: the  
umbrellas in 1D can force the sampling in the 2nd-D to surmount energy  
barriers that might be circumvented in unrestrained sampling.


Chris.

-- original message --

Qian Wang wrote:
Hi,

I used umbrella sampling method to restrain the distance of two
molecules at several distances. Then I can use g_wham to get the free
energy as a function of the distance. Is there any way that I can get
the free energy as a function of another parameter? 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

[gmx-users] umbrella sampling

[chris . neale]

I've never actually reweighted with F-values. These are the values by  
which the component PMFs are shifted during WHAM after debiasing the  
umbrella potentials.


Look at equations 6 and 7 in "The multicanonical weighted histogram  
analysis method for the free-energy landscape along structural  
transition paths" Satoshi Ono, Nobuyuki Nakajima, Junichi Higo and  
Haruki Nakamura


My best guess is that you do it like this:

1. Run 1-D wham and obtain the F-values.

2. create an empty 3D histogram structure for z, rgyrA, and rgyrB.

3. Go through the data from each window and add a value to the  
appropriate location for [z,rgyrA,rgyrB]. This value is not the  
integer count 1. This value that you add to each histogram bin is a  
real number:


  1 * exp(-(1/RT)*-F) * exp(-(1/RT)*0.5K*(z-z0)^2

4. Now convert your histogram of probabilities to free energies dG=-RTln(P)

** To test this, project your 3-D free energy profile onto each  
dimension successively and compare it to (z) the profile from 1-D  
wham, and (rgyrA or rgyrB) the profile from 2-D wham with zero force  
constants as I outlined in my previous post.


PLEASE NOTE: this email outlines an *idea* of how to accomplish what  
you want. You must check the match for yourself.


Chris.

Qian Wang qwang at mail.uh.edu
Sat Jul 30 18:05:33 CEST 2011

* Previous message: [gmx-users] umbrella sampling
* Next message: [gmx-users] OPLSAA parameters
* Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

Thanks very much.
I am sorry I did not describe my question very clear. In my system, I  
am going to obtain three things. (1) the free energy versus distance,  
which I can obtain from g_wham. (2) the free energy versus, eg, the  
radius of gyration of the object A. (3) the 2-D free energy as a  
function of the radius of gyration of A and the radius of gyration of  
B. I do not know how to get (2) and (3).
Your first method can help me obtain (2), but I can not obtain (3)  
because then it needs a 3-D WHAM code, am I correct? For your second  
method, I do not quite understand how to re-weighting using the  
F-values from 1D-WHAM. Could you please explain it? Thanks.


Sincerely,
Qian

- Original Message -
From: chris.neale at utoronto.ca
Date: Friday, July 29, 2011 11:00 pm
Subject: [gmx-users] umbrella sampling
To: gmx-users at gromacs.org

format your data for 2D-WHAM with 1D being the distance and the  
2nd-D being your other coordinate of interest. Specify a value of  
zero for the force constants for your 2nd-D. Run 2D-WHAM. Boltzmann  
project the 2D PMF onto your 2nd-D.


I think you can also do essentially the same thing by re-
weighting using the F-values from 1D-WHAM, but I find the above  
method to be the simplest. It also provides you with a 2D free  
energy profile, which can be informative both biologically and to  
indicate on sampling problems.


Note that you're very likely going to run into convergence problems  
since your 2nd-D will rely on brute-force to converge, and worse:  
the umbrellas in 1D can force the sampling in the 2nd-
D to surmount energy barriers that might be circumvented in  
unrestrained sampling.


Chris.

-- original message --

Qian Wang wrote:



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] spatial distribution function in a binary solvent mixture

[chris . neale]

I think that you have a misconception about what g_spatial does. For a  
system with many type A and many type B, you need to average over all  
of one type as the central solute to compute an rdf, and perhaps that  
is what you want for your sdf. g_spatial, however, does not do any  
fitting. g_sdf did that. You can obtain an old version of the source  
(4.0.7 should have it I think) if you want to use g_sdf. I have never  
used g_sdf myself.


In case that doesn't answer your question, then let me make one more point:

g_spatial requires that a bin exists for every count. Thus either (a)  
find a large memory node somewhere or (b) pre-process your trajectory  
using trjorder to order the solvent molecules and then keep only the N  
closest, then run g_spatial. But again, I think that this will not  
provide what you are looking for but is likely to give you a smear  
over your box.


Chris.

-- original message --

  Dear all,

 I'm working with a binary solvent mixture containing 2000 molecules (1800
 type-A + 200 type-B). Both the types of solvent molecules have similar
structure (they are both diatomic molecules) except the polarity.I'm
trying to calculate sdf of type-B solvent molecules. I followed the
step-by-step instructions from the manual using g_spatial.
 I'm trying to reduce the bin width to 0.05A. The *.cube file generated
with a bin width of 0.09A is already 4.9GB in size. As I'm more
interested in the first solvation shell around the type-B solvent
molecule, I was wondering if I could find a way to control the maximum
radius of the sphere around central molecule (center of coordinates?)?
Also, can anyone please let me know how g_spatial deals with the angular
part of sdf?

(I've searched a lot but could not gather enough meaningful information)

regards,
Debasmita




--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] spatial distribution function in a binary solvent mixture

[chris . neale]

Here's the totally wasteful way for you to get what you want:

1. for each molecule in type A: trjorder everything around that  
molecule and output only the closest N atoms of type A and the closest  
M atoms of type B. Ensure that the coordinate order is: your central  
molecule of type A is first, the remaining type A second, and the type  
B third.


2. concatenate all of these trajectories into one single large trajectory

3. trjconv -fit trans+rot , fitting only to the single central molecule.

4a. g_spatial for the non-central type A
4b. g_spatial for type B

* now repeat switching types A and B in the above steps to get the  
SDFs arouns type B.


A much better idea is for you to modify g_spatial to do this all  
within one loop. But if you don't know C then you are out of luck there.


g_sdf may not work because you are underdetermined (there is a  
symmetry about the long axis that you can not specify). But perhaps  
you could use that to your advantage and select some totally unrelated  
atom as the third atom for the fit. My best guess is that this would  
be no worse than g_spatial.


A final option is for you to use g_mindist to output all of the  
contacts and then run over that with your own script to process the  
data into an SDF and output it in cube format.


Chris.

-- original message --

 Thanks Chris. I presume g_sdf won't be helpful for my system.

[Hide Quoted Text]
I think that you have a misconception about what g_spatial does. For a
system with many type A and many type B, you need to average over all
of one type as the central solute to compute an rdf, and perhaps that
is what you want for your sdf. g_spatial, however, does not do any
fitting. g_sdf did that. You can obtain an old version of the source
(4.0.7 should have it I think) if you want to use g_sdf. I have never
used g_sdf myself.

In case that doesn't answer your question, then let me make one more
point:

g_spatial requires that a bin exists for every count. Thus either (a)
find a large memory node somewhere or (b) pre-process your trajectory
using trjorder to order the solvent molecules and then keep only the N
closest, then run g_spatial. But again, I think that this will not
provide what you are looking for but is likely to give you a smear
over your box.

Chris.

-- original message --

Dear all,

   I'm working with a binary solvent mixture containing 2000 molecules
(1800
   type-A + 200 type-B). Both the types of solvent molecules have similar
structure (they are both diatomic molecules) except the polarity.I'm
trying to calculate sdf of type-B solvent molecules. I followed the
step-by-step instructions from the manual using g_spatial.
   I'm trying to reduce the bin width to 0.05A. The *.cube file generated
with a bin width of 0.09A is already 4.9GB in size. As I'm more
interested in the first solvation shell around the type-B solvent
molecule, I was wondering if I could find a way to control the maximum
radius of the sphere around central molecule (center of coordinates?)?
Also, can anyone please let me know how g_spatial deals with the angular
part of sdf?

(I've searched a lot but could not gather enough meaningful information)

regards,
Debasmita

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Constraints not working in pull code (sometimes, sometimes not)

[chris . neale]

I agree that Justin is probably correct, although constraints should  
technically work just fine with a highly dynamic reference group. Any  
problems should show up as the system blowing up, but not in correctly  
setting the position. I think that the problems can arise, however,  
when you have multiple competing constraints. You might, for example,  
try without constraining the octane bonds -- I think that you should  
get perfect match in that case. Also, it is worth comparing in single  
and double precision.


Still, I am not sure that you actually did a proper comparison for the  
following 2 reasons:


1. you have apparently different masses for the pulled group:

Pull group 1:15 atoms, mass   194.194
Pull group 1:15 atoms, mass   146.146


2. you used very different starting depth offsets: -2.81855 is not  
very close to -1.60929


PS: I would personally prefer if you avoided jumping directly to a bug  
report unless you are sure of it. There is always at least the  
possibility that you are not using the code correctly.


Chris.


Dear all,

we are trying to simulate molecule in specified depth in octanol to
calculate logP.
for this purpose we have used box with slab of octanol, molecule in
specified distance in z-axis and  this mdp:

integrator  = md
dt  = 0.002
tinit   = 0
nsteps  = 5125000   ; 10.250 ns
nstcomm = 1
comm_mode   = Linear
; Output parameters
nstxout = 0 ; every 1 ns
nstvout = 0
nstfout = 0
nstxtcout   = 500   ; every 1 ps
nstenergy   = 500
; Bond parameters
constraints = all-bonds
; Single-range cutoff scheme
nstlist = 5
ns_type = grid
rlist   = 1.4
rcoulomb= 1.4
rvdw= 1.4
; PME electrostatics parameters
coulombtype = PME
fourierspacing  = 0.12
fourier_nx  = 0
fourier_ny  = 0
fourier_nz  = 0
pme_order   = 4
ewald_rtol  = 1e-5
optimize_fft= yes
; Berendsen temperature coupling is on in two groups
Tcoupl  = V-rescale
tc_grps = System
tau_t   = 0.1
ref_t   = 310
; Pressure coupling is on
Pcoupl  = berendsen
pcoupltype  = anisotropic   ; anisotropic pressure coupling
tau_p   = 10
compressibility = 0 0 4.5e-5 0 0 0  ; since octanol tries to separate from
water and shortening xy plane
ref_p   = 1.0 1.0 1.0 0 0 0
; Pull code
pull= constraint;
pull_geometry   = distance  ; Pull along the vector connecting the two
groups. Components can be selected with pull_dim.
pull_dim= N N Y ; put potential only to axis z
pull_start  = yes   ; define initial COM distance > 0
pull_ngroups= 1 ; one group to pull
pull_group0 = OCT   ; reference group (can be empty to set it to
0,0,0)
pull_group1 = DRG
pull_rate1  = 0 ; 0.01 nm per ps = 10 nm per ns = 10 m per s
pull_k1 = 2000  ; kJ mol^-1 nm^-2

however when I list pull-x.xvg I obtain:
@ s0 legend "0 Z"
@ s1 legend "1 dZ"
0.  6.06871 -2.81855
0.0200  6.07793 -2.81861
0.0400  6.08787 -2.81856
0.0600  6.08462 -2.81855
0.0800  6.10028 -2.82007

The situation is even more strange when almost the same mdp with DOPC
membrane works fine as can be seen in pull-x.xvg:
0.  3.72735 -1.60929
0.0200  3.72727 -1.60929
0.0400  3.7272  -1.60929
0.0600  3.72716 -1.60929
0.0800  3.72715 -1.60929
0.1000  3.72717 -1.60929

Differences in the working case in mdp are only this:
compressibility = 4.5e-5 4.5e-5 4.5e-5 0 0 0

ref_p   = 1.0 1.0 1.0 0 0 0
pull_group0 = DOPC

Strangerly enough - the same happens in GMX4.0.7 and in GMX4.5.3

I wonder where the difference is, as both logs states that
Will apply constraint COM pulling in geometry 'distance'
between a reference group and 1 group
Pull group 0:  6912 atoms, mass 100624.859
Pull group 1:15 atoms, mass   194.194

Will apply constraint COM pulling in geometry 'distance'
between a reference group and 1 group
Pull group 0:  9570 atoms, mass 124631.453
Pull group 1:15 atoms, mass   146.146

Where should I look and repair?

Thank you in advance
Karel




--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] Constraints not working in pull code (sometimes, sometimes not)

[chris . neale]

Dear Krapnik:

1. please make the test cases identical, that means doing a simulation  
that you may not really be interested in so that the solutes are the  
same for both lipid and octane. It also means setting the  
displacements to similar values in my opinion (because perhaps the  
problem is based on your box size along z and choice of pull_pbcatom  
for example). There is *something* strange happening, but let's  
eliminate all other possibilities before assuming that there is a  
problem with the source code.


2. Once you have a setup in which the only difference is octane vs.  
lipid (really the only difference please -- even add your octane  
pressure settings to the lipid simulation) and you find strange  
behaviour for octane and not lipid, then please attempt that same  
system with pull=umbrella in place of pull=constraints and see if the  
solute also drifts far from its initial position.


3. Your initial definition of the problem was a little misleading.  
Both Justin and I focused on third decimal place differences, but now  
you are saying that there is nm-scale motion along z. Please report  
those values for your test systems, which are outlined above, after  
you redo your tests with identical settings other than the switch from  
lipid to octane. -- I no longer thing that a double precision test is  
necessary if you are seeing nm-scale movement along z.


4. Perhaps I was a little hasty to complain about a bug-report style  
post. I did like all of the information that you provided. I just got  
the impression from your title and the end of your post that you  
thought it was a problem with gromacs (which it may still be) and I  
think that this is generally a counter-productive attitude because it  
hinders one from making the proper tests.


Chris

-- original message --

Re: Constraints not working in pull code (sometimes, sometimes not)
Krapnik krapnik at gmail.com
Thu Aug 11 09:53:24 CEST 2011

* Previous message: [gmx-users] append files in Gromacs 3.3
* Next message: [gmx-users] Constraints not working in pull code  
(sometimes, sometimes not)

* Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

Dear Justin
Thank you for your responses, but it still is not clear to me why is that
happening


I don't see any real problem here.  The DOPC membrane probably holds its

shape

better than the octanol slab and thus the reference position (which

determines

the constraint) is less mobile.


That is certain. When we were running free simulation, the box started to
narrow. Therefore we have frozen the box in xy plane.


As such, the constraint is more stable in DOPC
than in octanol.  The slight bump in the octanol system amounts to a

change of

0.16%, which I'd say is quite small.


The trouble is, that during next 10 ns the distance was between 2.0 - 4.0
and molecule moved outside of the ocatnol slab and inside, which is then
useless for calculation of PMF.


You're also looking at the first few
frames of the simulation, which may amount to equilibration, but you

haven't

mentioned if you've done anything prior to this run.


Previously, a free simulation with molecule on slab was run for 20 ns, from
which frame with small molecule at specidied depth was taken, so I suppose,
that the start is equilibrated.


As long as there is no systematic change in position, I'd say there's no

problem.

Unfortunately, there is as I have said before.

-Justin

Dear Chris,
thank you for responses as well.


I agree that Justin is probably correct, although constraints should
technically work just fine with a highly dynamic reference group. Any
problems should show up as the system blowing up, but not in correctly
setting the position.


I hope, that the position is set up ok, as it stays steady in DOPC.


I think that the problems can arise, however,
when you have multiple competing constraints. You might, for example,
try without constraining the octane bonds -- I think that you should
get perfect match in that case.


I will try.
The erroneous movement of molecule seems to be enhanced within XY-plane
frozen.
The curious thing is however, that when I did not constrain size of
xy-plane, then the distance to the COM of slab is not changing much (only
last digit changes), which in comparison to frozen XY-plane, where there is
huge movement above the slab and within is rather ok to me. Unfortunately, I
have the issue with narrowing of the box in the unconstrained XYplane
simulation, which also destroys any attempts for PMF as I have mentioned
before.


Also, it is worth comparing in single
and double precision.


I will try this either.


Still, I am not sure that you actually did a proper comparison for the
following 2 reasons:
1. you have apparently different masses for the pulled group:
Pull group 1:15 atoms, mass   194.194
Pull group 1:15 atoms, mass   146.146


The difference is not much in my opinion, while molecules tested are similar
in size and are run in

[gmx-users] vsites with NAC and opls

[chris . neale]

This is fixed, just posting for posterity.

If one wants to run pdb2gmx with -vsites hydrogens on gromacs 4.5.4 or  
4.0.7 while using the oplss ff, there is the error message:


Fatal error:
Can't find dummy mass for type opls_242 bonded to type opls_238 in the  
virtual site database (.vsd files). Add it to the database!


To fix this, one must then add the following line:

   opls_242  opls_238   MCH3A

to the [ CH3 ] section of oplsaa.ff/aminoacids.vsd

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] suggestion that mdrun should ensure npme < the number of processes

[chris . neale]

Currently, gromacs4.5.4 gives a segfault if one runs mpirun -np 8  
mdrun_mpi -npme 120 with no warning of the source of the problem.


Obviously npme>nnodes is a bad setup, but a check would be nice.

Thank you,
Chris.

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

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

[gmx-users] GROMACS 4.5.4 keep crashing all the time.

[chris . neale]

You'll need to provide a much better report than this if you want to  
receive any useful help.


Copy and paste the exact commands of what you did
Copy and paste the exact log file and error messages

Do this for 4.0.7 and 4.5.4, for which I trust that you have been  
using exactly identical test systems. If not, then please try it again  
while conserving the system.


Chris.

-- original message --

Hi Matthew,

Thanks for the reply. First, I don't thinks it is poorly compiled bin, as I
used both GROMACS compiled by me (on my machine) or by the sys-admin on
Linux cluster or blue-gene.

The simulations using 4.5.4 crashed giving LINCS error, which is not the
case with 4.0.7 (using the same mdp and pdb files). Second, I don't thinks
it is poorly compiled bin, as I used both GROMACS compiled with me (on my
machine) or by the sys-admin on Linux cluster/blue-gene.

cheers,
Itamar

On 18 August 2011 01:48, Matthew Zwier  wrote:


Could be a system blowing up, or perhaps a mis-compiled binary.  What
error messages do you get when the crash occurs?

On Tue, Aug 16, 2011 at 9:48 PM, Itamar Kass 
wrote:
> Hi all GROMACS useres and developers,
>
> I am interesting in simulating a small protein (~140 aa) in water, with
and without Ca ions. In order to do so, I had used version 4.5.4. I had
solvate the protein in water, add ions to naturalise the systems,
equilibrated the systems and then tried productive runs. Now, no matter what
I did, it crashed after few ps's of free MD or during the PR runs.
>
> Few of the things I had tried are:
> 1. Running the simulations on different systems (OSX, linux or
blue-gene).
> 2. Using single or double precision versions.
> 3. An equilibration stage, 1ns long with a time-step of 1fs, during which
the restrained forces where gradually reduced from 1000 to 50 kJ mol-1 nm-2.
> 4. Running part of the equilibration stage as NVT and then switched to
NPT.
> 5. Started from different x-ray structures, with resolutions differ from
2.5 to 1.7 Angstrom.
>
> Finally I had moved back to 4.0.7 which worked like charm. I wonder if
someone else had encounter something like this. Attached please find the mdp
files I used.
>
> All the best,
> Itamar.
>
>
>
>



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] GROMACS 4.5.4 keep crashing all the time

[chris . neale]

run an EM with flexible water. I often find that this is the only way  
to get a stable system. 500 steps of steep with define=-DFLEXIBLE (or  
different depending on your water model I think) should be enough.


Chris.

Hi Chris and Justin,

On 18/08/2011, at 9:36 AM, Justin A. Lemkul wrote:




chris.neale at utoronto.ca wrote:
You'll need to provide a much better report than this if you want  
to receive any useful help.

Copy and paste the exact commands of what you did
Copy and paste the exact log file and error messages




The command I had used (for both 4.0.7 and 4.5.4) are:

pdb2gmx -f 1EXR.pdb -o 1EXR.gro -water spc -ignh -p 1EXR.top -i 1EXR.itp
editconf -f 1EXR.gro -o 1EXR_box.gro -c -d 1.4 -bt cubic

genbox -cp 1EXR_box.gro -cs spc216.gro -o 1EXR_solv.gro -p 1EXR.top -seed 5

grompp -f em.mdp -c 1EXR_solv.gro -p 1EXR.top -o ions.tpr
genion -s ions.tpr -o 1EXR_solv_ions.gro -p 1EXR.top -pname NA+ -nname  
CL- -np 17 -seed 66


grompp -f em.mdp -c 1EXR_solv_ions.gro -p 1EXR.top -o em.tpr

mpirun mdrun_d_mpi -v -s em.tpr -x 1EXR_em.xtc -e 1EXR_em.edr -g  
1EXR_em.log -c 1EXR_em.gro


grompp -f pr.mdp -c 1EXR_em.gro -p 1EXR.top -o pr.tpr

mpirun mdrun_d_mpi -v -stepout 1000 -s pr.tpr -x 1EXR_pr.xtc -e  
1EXR_pr.edr -g 1EXR_pr.log -o 1EXR_pr.trr -c 1EXR_pr.gro


grompp -f md_start.mdp -c 1EXR_pr.gro -p 1EXR.top -o md.tpr

mpirun mdrun_d_mpi -v -stepout 1 -s md.tpr -x 1EXR_md.xtc -e  
1EXR_md.edr -g 1EXR_md.log -o 1EXR_md.trr -c 1EXR_md.gro


The em.mdp:
integrator   =  steep
maximum number of steps to integrate
nsteps   =  5
;Energy minimizing stuff
emstep   =  0.001
emtol=  100.0

; OPTIONS FOR BONDS
constraints  =  none


; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout  =  1000
nstvout  =  1000
nstfout  =  1000
; Output frequency and precision for xtc file
nstxtcout=  1000
xtc-precision=  1000
; Energy monitoring
energygrps   =  system

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist  =  5
; ns algorithm (simple or grid)
ns-type  =  Grid
; nblist cut-off
rlist=  0.8

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype  = Reaction-Field
rcoulomb = 1.4
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon-rf   = 62
; Method for doing Van der Waals
vdw-type = Cut-off
; cut-off lengths
rvdw = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = No

; Center of mass control
nstcomm  =  1000
; Periodic boundary conditions
pbc  =  xyz
; Mode for center of mass motion removal
comm-mode=  linear
; Groups for center of mass motion removal
comm-grps=  system

the pr.mdp:
define   =  -DPOSRES

integrator   =  md
dt   =  0.002   ; ps !
nsteps   =  5   ; total 100ps.

; OPTIONS FOR BONDS
; Constrain control
constraints  =  all-bonds
; Do not constrain the start configuration
continuation  = no
; Type of constraint algorithm
constraint-algorithm =  lincs

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout  =  10
nstvout  =  10
nstfout  =  10
; Output frequency and precision for xtc file
nstxtcout=  5000
xtc-precision=  1000
; Energy monitoring
energygrps   =  Protein Non-protein
nstenergy=  5000

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist  =  5
; ns algorithm (simple or grid)
ns-type  =  Grid
; nblist cut-off
rlist=  0.8

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype  = Reaction-Field
rcoulomb = 1.4
epsilon_rf   = 62; As suggested at J Comput Chem. 2004 
Oct;25(13):1656-76.
vdw-type = Cut-off
; cut-off lengths
rvdw = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = no

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl   = berendsen
; Groups to couple separately, time constant (ps) and reference  
temperature (K)

tc-grps  = Protein Non-Protein
tau-t= 0.1 0.1
ref-t= 300 300
; Pressure coupling
Pcoupl   = berendsen
Pcoupltype   = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p=  1.0
compressibility  =  4.5e-5
ref_p=  1.0
; Generate velocites is on at 290K 

[gmx-users] GROMACS 4.5.4 keep crashing all the time

[chris . neale]

OK, here's my last few ideas:

1. Please try to repeat this with gen_vel set to the same value as  
your temperature coupling


2. Can you reproduce this in serial?

3. Can you reproduce this with the sd integrator?

4. Can you reproduce this with a simpler system? protein in vacuum or  
just water or remove the ions, etc?


5. Take the output .gro from 4.0.7 that ran fine for X ns and run it  
under 4.5.4. Do you get the same lincs warnings?


6. Also, note that you are getting warnings and the run does not  
actually crash but just stops after too many warnings. So what are  
atoms 981 and 982? Does their motion look different in an important  
ways between the 4.0.7 and 4.5.4 trajectories?


Chris.

-- original message --

Hi Chris,

thanks for the advice, I have to say I tried this as well without any success.

Itamar

On 18/08/2011, at 11:11 AM, chris.neale at utoronto.ca wrote:

run an EM with flexible water. I often find that this is the only  
way to get a stable system. 500 steps of steep with  
define=-DFLEXIBLE (or different depending on your water model I  
think) should be enough.


Chris.

Hi Chris and Justin,

On 18/08/2011, at 9:36 AM, Justin A. Lemkul wrote:




--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. 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] more than 100% CPU

[chris . neale]

Slightly off-topic, but for all but the smallest systems, I get a  
further 10% efficiency by running a 16-process mpi job on an 8-core  
machine. I suspect that the story is the same with threads. Thus, on a  
16-core node, you could try starting 32 threads. Note that this will  
report a 2x larger "efficiency" as you were discussing before, but by  
checking the ns.day you will see that the benefit is between 5% and 16%


Chris.

It means you are running a lot of parallel processes, but that does  
not translate into a linear increase in speed.  So faster, but not 16X.


Warren Gallin

On 2011-08-18, at 5:36 PM, Park, Jae Hyun nmn wrote:


Thank you, Warren.
Does that mean 16 times faster ?

Jae H. Park

-Original Message-
From: gmx-users-bounces at gromacs.org [mailto:gmx-users-bounces at  
gromacs.org] On Behalf Of Warren Gallin

Sent: Thursday, August 18, 2011 7:26 PM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] more than 100% CPU

I believe that the current version of GROMACS supports threading,  
which does not require mpi.


So mdrun is running threads at 100% of the activity of each of your  
16 nodes, hence 1600%.


Warren Gallin

On 2011-08-18, at 5:19 PM, Park, Jae Hyun nmn wrote:


Hi GMX users,

I installed GMX 4.5.3 recently.
But, when I just execute mdrun (without mpi, I did not installed  
mpi-version of mdrun), the use of CPU appears more than 100% ("top"  
command in LINUX). How is it possible?
For example, I am using 16-node machine. And if I simply run  
"mdrun", then  use of CPU is 1600%!!.

The simulation runs well and the results looks reasonable.
Is there anybody who can teach me what is happening? I would deeply  
appreciate.





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