Re: [gmx-users] floating point exception in .xtc file

[francesco oteri]

Hi Christopher,
you can try to use the program gmx_rescue, by Marc Baaden to recovery your
trajectory.

Below there is the adderess:
http://baaden.free.fr/soft/compchem.html

Francesco

2012/6/9 Mark Abraham 

>  On 9/06/2012 7:27 AM, Christopher Neale wrote:
>
>  Dear Users:
>
>  I have a 500 ns trajectory of 65 GB that gives a floating point
> exception when I run it through gmxcheck or trjcat (generated and analyzed
> with gromacs 4.5.5). Has anybody encountered this? I ran mdrun with -append
> so this is the xtc resulting from months of simulation of a 1,000,000 atom
> system. If I run trjconv -f md.xtc -b 20, where the floating point
> exception occurred around t=18 ps in gmxcheck, then I can extract the
> readable frames and repair around the damaged section. Still, I'd rather
> not lose any data and I had thought that the new default -append option to
> mdrun checked for these types of problems at runtime.
>
>
> I've no idea what might happen when some file-system transient occurs
> mid-simulation, but if mdrun has managed to compute a checksum on an
> incomplete file and stored that in the checkpoint, then the append
> mechanism can be none the wiser. The check upon restart is that the
> checksum matches, not that the checksum is computed on a file whose
> properties would satisfy gmxcheck.
>
> Note also that some file systems that do not support file locking and this
> is known to cause issues (Redmine 924), but I don't know if this is related
> to your observation.
>
> 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
>



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

Re: [gmx-users] floating point exception in .xtc file

[Mark Abraham]

On 9/06/2012 7:27 AM, Christopher Neale wrote:

Dear Users:

I have a 500 ns trajectory of 65 GB that gives a floating point 
exception when I run it through gmxcheck or trjcat (generated and 
analyzed with gromacs 4.5.5). Has anybody encountered this? I ran 
mdrun with -append so this is the xtc resulting from months of 
simulation of a 1,000,000 atom system. If I run trjconv -f md.xtc -b 
20, where the floating point exception occurred around t=18 ps 
in gmxcheck, then I can extract the readable frames and repair around 
the damaged section. Still, I'd rather not lose any data and I had 
thought that the new default -append option to mdrun checked for these 
types of problems at runtime.




I've no idea what might happen when some file-system transient occurs 
mid-simulation, but if mdrun has managed to compute a checksum on an 
incomplete file and stored that in the checkpoint, then the append 
mechanism can be none the wiser. The check upon restart is that the 
checksum matches, not that the checksum is computed on a file whose 
properties would satisfy gmxcheck.


Note also that some file systems that do not support file locking and 
this is known to cause issues (Redmine 924), but I don't know if this is 
related to your observation.


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

Re: [gmx-users] Re: Re: energy conservation: shift vs shifted user potential

[Mark Abraham]

On 9/06/2012 2:05 AM, Anja Kuhnhold wrote:

Ok, I try to explain it again :)
I know the format of the user tables (x, f(x), -f'(x), g(x), -g'(x), h(x), 
-h'(x);
f is coulomb, g dispersion, h repulsion).
For testing the simple cut-off version, my table looks like the
example in the installation.
So I run a simulation with vdwtype=cut-off
and the same with vdwtype=user.
The results are the same.

The problem is, that I want to use shifted forces and potentials within my 
tables.
(Aim is better energy conservation)
In the manual (ch. 4.1.5) is explained how the shift functions
for forces and potentials have to be.
So for -g(x) and h(x) I use the potential from the manual, Phi, with alpha=6 
and 12, respectivly.
And for g'(x) and -h'(x) I use the shifted force from the manual, Fs, with 
alpha=6 and 12, respectivly.
But when I run the simulation with this table on the one hand
and with vdwtype=Shift on the other hand,
the results are not comparable, mostly the energy is not conserved by using the 
table.


The force is the negative derivative of the potential. If you change 
only the force columns then you have made the interaction piecewise 
discontinuous at every segment. So when GROMACS uses both the values and 
derivatives to convert your table into the internal format used (section 
6.7.1), it comes up with junk. There should be a warning printed if the 
force column is not close enough to the numerical derivative of the 
potential column, but maybe this change is close enough that you avoid 
that. Anyway, this makes energy basically impossible to conserve. You 
should generate all your table values from the actual function you wish 
to implement, e.g. with a spreadsheet or scripting language.


Mark



Can you follow me?

Anja



- Ursprüngliche Nachricht -
Von: gmx-users-requ...@gromacs.org
Datum: Freitag, Juni 8, 2012 15:55
Betreff: gmx-users Digest, Vol 98, Issue 52
An: gmx-users@gromacs.org


- Ursprüngliche Nachricht -
Von Mark Abraham
Datum Fri, 08 Jun 2012 23:40:35 +1000
An Discussion list for GROMACS users
Betreff Re: [gmx-users] Re: energy conservation: shift vs shifted user
potential
On 8/06/2012 5:59 PM, Anja Kuhnhold wrote:

Hello all,

can someone give me a hint, please?
Do you need more information?
Has anyone had a similar problem.

I really need to figure that out, because I will simulate other systems
which can be run only with user tables.

Anja



-

- Original Message -
- Original Message -
 From Anja Kuhnhold
Date Wed, 06 Jun 2012 15:07:40 +0200
To gmx-users@gromacs.org
Subject [gmx-users] energy conservation: shift vs shifted user potential
Dear gmx-users,

I have a problem concerning energy conservation when using user
potentials (tables).
I'm using gromacs 4.5.4
I simulate a simple Lennard-Jones(6-12) +Fene polymer melt (1600
chains a 10 beads in a 26x26x26 periodic box).

I tried different vdwtypes (cutoff always 3.24):
The cut-off version does not conserve energy -- okay.
The shift and switch versions conserve energy -- fine.

Now I wanted to do the same with user tables:
Simple Lennard-Jones table gives really the same results as the
cut-off version -- good.

But if I use a table with shifted Lennard-Jones potential it is not
comparable to the shift version
and the energy is not conserved -- ?

I use a shift function as written in the manual (chapter 4.1.5) --
there must be a factor alpha added in the constants A and B --
(r1=0).

The parameters are the same for shift version and shifted user version.

Has someone an idea why the shifted user potential doesn't work in
this way?

We've no real idea what you've done... manual 6.7.2 describes the
required format and there are (unshifted) examples in your
installation
under $GMXLIB/share/top/table*.xvg.

Makr


Here is the mdp:


integrator  = md-vv
dt  = 0.0035
nsteps  = 1000
nstxout = 1
nstvout = 1
nstfout = 1
nstlog  = 1
ns_type = grid
pbc = xyz
rvdw= 3.24
rlist   = 3.6
tcoupl  = no
tc-grps = System
tau_t   = 2.0
ref_t   = 127.2717
vdwtype = user;Shift
rcoulomb= 3.6;2.24;1.12
coulombtype = Cut-off
rvdw-switch = 0.0

energygrps  = bead
energygrp_table = bead bead


Thanks in advance
Anja



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

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

[gmx-users] vsite as interaction site for COM of benzene

[Thomas Schlesier]

Dear Richard,

thanks for the idea, i will try this.

Relating to your side note:
I perform pulling simulations of a calixarene-catenane dimer (both together 600 
atoms)
and i mesitylene as a solvent. To access longer time-scales I wanted to 
coarse-grain the
solvent only. My supervisor said that the mesitylene rotates on such a short 
time-scale,
that it would be ok, so coarse-grain it as one particle. But i also thought 3 
particles
would be an option.
But since it is easier to coarse-grain the mesitylene as one particle, i 
started with that,
at this stage, my first interest to get a somewhat reasonable simulation, which 
runs, since
it is my first project involving coarse-graining and vsites.
If the results are not good enough, i still can make the system more complex :) 
but it is easier
to start with the simple case.
I'm now at the to construct the AA-CG interaction. For this is started with an 
atomistic mesitylene
in coarse-grained mesitylene (strictly speaking it's still an CG-CG 
interaction, but the single mesitylene
has an atomistic representation, as later the calixarene).

For the question here i used benzene instead of mesitylene, since it is smaller 
but the concept is the same.






Dear Thomas,

Whilst in principle a constraint should do what you are asking, you might be
better off using multiple vsites to solve the problem so that there is an
algebraic as apposed to numerical mapping of the forces. My immediate
thought is to use two type 3 (or if necessary type 3fd) virtual sites to
define two points near the centre of the benzene. You could then place your
tabulated potential on a type 2 virtual site between these two. I believe
that should algebraically link all the atomic sites together making the
process more stable and it should be quite easy to transfer it to other
systems.

On a separate note why if you are coarse graining the benzenes to a single
point particle are you then maintaining all the atoms? Could you go a stage
further and map the atoms onto 3 non interacting point masses? That would
maintain angular momentum etc. whilst reducing the computational overhead.
Although it becomes very complex to implement if your benzene's are in the
backbone of a polymer.

Hope that helps,

Richard

On 08/06/2012 18:18, "Thomas Schlesier"http://lists.gromacs.org/mailman/listinfo/gmx-users>>  wrote:


/  Hi all,

/>/
/>/  i have a more conceptional question, for using vsites as
/>/  interaction-centers for coarse-grained particles:
/>/
/>/  First the simple case:
/>/  I want to simulate one benzene molecule (atomistic - aa) in
/>/  coarse-grained (cg-) benzene (each benzene molecule as a single
/>/  particle). For the cg-cg interaction of the benzene i have calculated a
/>/  table potential.
/>/
/>/  But for the aa-cg interaction I'm a little bit clueless:
/>/  I want to use the COM of the aa-benzene as an interaction-site. But
/>/  there is no vsite, which i could construct from 6 atoms.
/>/  Side-note: The vsite would not interact with the with the atoms, only
/>/  with cg-benzenes
/>/
/>/  So there are two ideas:
/>/
/>/  1) Using one '3out'-vsite, which is construted from 3 non-neighbour
/>/  atoms (if we had mesitylene instead of benzene, it would be the three
/>/  C-atoms to which only hydrogens are bound). The the force from the aa-cg
/>/  interaction would be distributed to these three atoms and i would hope
/>/  that the bond-constraints would do 'the rest' (i.e. they would mimick
/>/  that the whole aa-benzene molecule interacts with the cg-particles)
/>/
/>/  2) Using two '3out'-vsites. First viste the same as in (1), the second
/>/  comes from the set of the other three C-atoms. These two vsites, should
/>/  be at the same position (more or less, for a perfect energy-minimized
/>/  benzene molecule it would be so; but i would assume the error is only
/>/  minimal). Since the potential is pair-additive, i would use for this
/>/  table potential only half of the potential, so that the potential of
/>/  both particles would add up to the 'true' potential. The forces would be
/>/  calculated from the 'half-potential'.
/>/
/>/  Big questions is now:
/>/  Are (1) and (2) equivalent?
/>/  (2) Seems like the 'right way' to do.
/>/  But (1) would be easier to implement :) But the big question is, if the
/>/  bond-constraints would correct the approximation, that only 3 out of 6
/>/  atoms have a direct interaction with the cg-particles.
/>/
/>/  Side note:
/>/  The real system is a little more complex. There are fracments of
/>/  molecules where i could not really construct two vsites which would have
/>/  more or less the same position. In this case, method (1) would really help.
/>/
/>/  Any ideas?
/>/  Greetings
/>/  Thomas
/

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

Re: [gmx-users] RE: Energy Minimization - not getting correct lipid area

[Justin A. Lemkul]

On 6/8/12 6:54 PM, Erica Hicks wrote:

How would I know what coordinate file to use? I have been following the
tutorial word for word. However, I did not get the confout.gro file that I
should have gotten after the minimization (mdrun -v -deffnm em) so I have
been using em.gro. Would this be a problem?


No, confout.gro assumes you're using default names.  The em.gro file you 
obtained is correct.  I have no idea how it's got 128 lipids in it though. 
Check its contents.  If it has 128 lipids, you used an incorrect coordinate file 
for grompp to prepare em.tpr.  That's about the only thing I can think that has 
gone wrong.


-Justin



From: Justin A. Lemkul [via GROMACS]
[ml-node+s5086n4998247...@n6.nabble.com] Sent: Friday, June 08, 2012 5:42 PM
To: Hicks, Erica Subject: RE: Energy Minimization - not getting correct lipid
area



On 6/8/12 6:33 PM, Erica Hicks wrote:


Hi,

bash-3.2$ perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5
area.dat Reading. Scaling lipids There are 128 lipids... with 50
atoms per lipid..

Determining upper and lower leaflet... 64 lipids in the upper... 64 lipids
in the lower leaflet

Centering protein Checking for overlap ...this might actually take
a while 100 % done... There are 2 lipids within cut-off range... 1 will
be removed from the upper leaflet... 1 will be removed from the lower
leaflet...

Writing scaled bilayer&   centered protein...


Calculating Area per lipid... Protein X-min/max: 2640 Protein
Y-min/max: 2539 X-range: 14 AY-range: 14 A Building 14 X 14 2D grid
on protein coordinates... Calculating area occupied by protein.. full
TMD.. upper TMD lower TMD Area per protein: 2 nm^2 Area per lipid:
10.4716089904762 nm^2

Area per protein, upper half: 1.75 nm^2 Area per lipid, upper leaflet :
10.475577244 nm^2

Area per protein, lower half: 2 nm^2 Area per lipid, lower leaflet :
10.4716089904762 nm^2

Writing Area per lipid... Done!

bash-3.2$ mdrun -v -deffnm em Back Off! I just backed up em.log to
./#em.log.6# Getting Loaded... Reading file em.tpr, VERSION 4.5.4 (single
precision) Starting 4 threads Loaded with Money

Making 1D domain decomposition 4 x 1 x 1

Back Off! I just backed up em.trr to ./#em.trr.6#

Back Off! I just backed up em.edr to ./#em.edr.6#

Steepest Descents: Tolerance (Fmax)   =  1.0e+03 Number of steps=
5 Step=0, Dmax= 1.0e-02 nm, Epot= -3.06212e+05 Fmax= 5.37269e+03,
atom= 3275 Step=1, Dmax= 1.0e-02 nm, Epot= -3.10498e+05 Fmax=
2.40668e+03, atom= 3024 Step=3, Dmax= 6.0e-03 nm, Epot= -3.10958e+05
Fmax= 2.89119e+03, atom= 4225 Step=5, Dmax= 3.6e-03 nm, Epot=
-3.12589e+05 Fmax= 1.46767e+03, atom= 1365 Step=6, Dmax= 4.3e-03 nm,
Epot= -3.13033e+05 Fmax= 3.83969e+03, atom= 1365 Step=7, Dmax= 5.2e-03
nm, Epot= -3.13599e+05 Fmax= 2.88720e+03, atom= 1365 Step=8, Dmax=
6.2e-03 nm, Epot= -3.13688e+05 Fmax= 4.93154e+03, atom= 1365 Step=9,
Dmax= 7.5e-03 nm, Epot= -3.14101e+05 Fmax= 4.72541e+03, atom= 1365 Step=
11, Dmax= 4.5e-03 nm, Epot= -3.14800e+05 Fmax= 9.90227e+02, atom= 1365

writing lowest energy coordinates.

Back Off! I just backed up em.gro to ./#em.gro.6#

Steepest Descents converged to Fmax<   1000 in 12 steps Potential Energy  =
-3.1480012e+05 Maximum force =  9.9022711e+02 on atom 1365 Norm of
force =  1.2180042e+02


bash-3.2$ perl inflategro.pl em.gro 0.95 DPPC 0 system_shrink1.gro 5
area_shrink1.dat Reading. Scaling lipids There are 128 lipids...


Something is wrong here.  You should only have 126 lipids since 2 were
removed before.  Did you minimize the correct coordinate file?


with 50 atoms per lipid..

Determining upper and lower leaflet... 64 lipids in the upper... 64 lipids
in the lower leaflet

No protein coordinates found...


This is also strange.  This, coupled with the lines above regarding the
number of lipids, suggest you're using the wrong coordinate file.


Writing scaled bilayer&   centered protein...


Calculating Area per lipid... Protein X-min/max: 1- Protein
Y-min/max: 1- X-range: -1 AY-range: -1 A Building
-1 X -1 2D grid on protein coordinates... Calculating area occupied
by protein.. full TMD.. upper TMD.. lower TMD..


The following areas are meaningless, for reasons associated with incorrect
content or processing thereof.


Area per protein: 0 nm^2 Area per lipid: 0.583197761890625 nm^2

Area per protein, upper half: 0 nm^2 Area per lipid, upper leaflet :
0.583197761890625 nm^2

Area per protein, lower half: 0 nm^2 Area per lipid, lower leaflet :
0.583197761890625 nm^2

Writing Area per lipid... Done!

To me, it looks like it is not looking at the 126 lipids that I updated in
the topology file. H


OK, I seem to remember 124 being what my system produced, but that's
irrelevant. Check out the lines above.


Erica

 From: Justin A. Lemkul [via GROMACS]
[[hidden email]] Sent: Friday, Ju

[gmx-users] RE: Energy Minimization - not getting correct lipid area

[Erica Hicks]

How would I know what coordinate file to use? I have been following the 
tutorial word for word. However, I did not get the confout.gro file that I 
should have gotten after the minimization (mdrun -v -deffnm em) so I have been 
using em.gro. Would this be a problem?

From: Justin A. Lemkul [via GROMACS] [ml-node+s5086n4998247...@n6.nabble.com]
Sent: Friday, June 08, 2012 5:42 PM
To: Hicks, Erica
Subject: RE: Energy Minimization - not getting correct lipid area



On 6/8/12 6:33 PM, Erica Hicks wrote:

> Hi,
>
> bash-3.2$ perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 
> area.dat
> Reading.
> Scaling lipids
> There are 128 lipids...
> with 50 atoms per lipid..
>
> Determining upper and lower leaflet...
> 64 lipids in the upper...
> 64 lipids in the lower leaflet
>
> Centering protein
> Checking for overlap
> ...this might actually take a while
> 100 % done...
> There are 2 lipids within cut-off range...
> 1 will be removed from the upper leaflet...
> 1 will be removed from the lower leaflet...
>
> Writing scaled bilayer&  centered protein...
>
>
> Calculating Area per lipid...
> Protein X-min/max: 2640
> Protein Y-min/max: 2539
> X-range: 14 AY-range: 14 A
> Building 14 X 14 2D grid on protein coordinates...
> Calculating area occupied by protein..
> full TMD..
> upper TMD
> lower TMD
> Area per protein: 2 nm^2
> Area per lipid: 10.4716089904762 nm^2
>
> Area per protein, upper half: 1.75 nm^2
> Area per lipid, upper leaflet : 10.475577244 nm^2
>
> Area per protein, lower half: 2 nm^2
> Area per lipid, lower leaflet : 10.4716089904762 nm^2
>
> Writing Area per lipid...
> Done!
>
> bash-3.2$ mdrun -v -deffnm em
> Back Off! I just backed up em.log to ./#em.log.6#
> Getting Loaded...
> Reading file em.tpr, VERSION 4.5.4 (single precision)
> Starting 4 threads
> Loaded with Money
>
> Making 1D domain decomposition 4 x 1 x 1
>
> Back Off! I just backed up em.trr to ./#em.trr.6#
>
> Back Off! I just backed up em.edr to ./#em.edr.6#
>
> Steepest Descents:
> Tolerance (Fmax)   =  1.0e+03
> Number of steps=5
> Step=0, Dmax= 1.0e-02 nm, Epot= -3.06212e+05 Fmax= 5.37269e+03, atom= 3275
> Step=1, Dmax= 1.0e-02 nm, Epot= -3.10498e+05 Fmax= 2.40668e+03, atom= 3024
> Step=3, Dmax= 6.0e-03 nm, Epot= -3.10958e+05 Fmax= 2.89119e+03, atom= 4225
> Step=5, Dmax= 3.6e-03 nm, Epot= -3.12589e+05 Fmax= 1.46767e+03, atom= 1365
> Step=6, Dmax= 4.3e-03 nm, Epot= -3.13033e+05 Fmax= 3.83969e+03, atom= 1365
> Step=7, Dmax= 5.2e-03 nm, Epot= -3.13599e+05 Fmax= 2.88720e+03, atom= 1365
> Step=8, Dmax= 6.2e-03 nm, Epot= -3.13688e+05 Fmax= 4.93154e+03, atom= 1365
> Step=9, Dmax= 7.5e-03 nm, Epot= -3.14101e+05 Fmax= 4.72541e+03, atom= 1365
> Step=   11, Dmax= 4.5e-03 nm, Epot= -3.14800e+05 Fmax= 9.90227e+02, atom= 1365
>
> writing lowest energy coordinates.
>
> Back Off! I just backed up em.gro to ./#em.gro.6#
>
> Steepest Descents converged to Fmax<  1000 in 12 steps
> Potential Energy  = -3.1480012e+05
> Maximum force =  9.9022711e+02 on atom 1365
> Norm of force =  1.2180042e+02
>
>
> bash-3.2$ perl inflategro.pl em.gro 0.95 DPPC 0 system_shrink1.gro 5 
> area_shrink1.dat
> Reading.
> Scaling lipids
> There are 128 lipids...

Something is wrong here.  You should only have 126 lipids since 2 were removed
before.  Did you minimize the correct coordinate file?

> with 50 atoms per lipid..
>
> Determining upper and lower leaflet...
> 64 lipids in the upper...
> 64 lipids in the lower leaflet
>
> No protein coordinates found...

This is also strange.  This, coupled with the lines above regarding the number
of lipids, suggest you're using the wrong coordinate file.

> Writing scaled bilayer&  centered protein...
>
>
> Calculating Area per lipid...
> Protein X-min/max: 1-
> Protein Y-min/max: 1-
> X-range: -1 AY-range: -1 A
> Building -1 X -1 2D grid on protein coordinates...
> Calculating area occupied by protein..
> full TMD..
> upper TMD..
> lower TMD..

The following areas are meaningless, for reasons associated with incorrect
content or processing thereof.

> Area per protein: 0 nm^2
> Area per lipid: 0.583197761890625 nm^2
>
> Area per protein, upper half: 0 nm^2
> Area per lipid, upper leaflet : 0.583197761890625 nm^2
>
> Area per protein, lower half: 0 nm^2
> Area per lipid, lower leaflet : 0.583197761890625 nm^2
>
> Writing Area per lipid...
> Done!
>
> To me, it looks like it is not looking at the 126 lipids that I updated in 
> the topology file. H

OK, I seem to remember 124 being what my system produced, but that's irrelevant.
  Check out the lines above.

> Erica
>
> 
> From: Justin A. Lemkul [via GROMACS] [[hidden email]]
> Sent: Friday, June 08, 2012 5:21 PM
> To: Hicks, Erica
> Subject: RE: Energy Minimization - not getting correct lipid area
>
>
>
> On 6/8/12 6:18 PM, Erica Hicks wr

Re: [gmx-users] RE: Energy Minimization - not getting correct lipid area

[Justin A. Lemkul]

On 6/8/12 6:33 PM, Erica Hicks wrote:

Hi,

bash-3.2$ perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat
Reading.
Scaling lipids
There are 128 lipids...
with 50 atoms per lipid..

Determining upper and lower leaflet...
64 lipids in the upper...
64 lipids in the lower leaflet

Centering protein
Checking for overlap
...this might actually take a while
100 % done...
There are 2 lipids within cut-off range...
1 will be removed from the upper leaflet...
1 will be removed from the lower leaflet...

Writing scaled bilayer&  centered protein...


Calculating Area per lipid...
Protein X-min/max: 2640
Protein Y-min/max: 2539
X-range: 14 AY-range: 14 A
Building 14 X 14 2D grid on protein coordinates...
Calculating area occupied by protein..
full TMD..
upper TMD
lower TMD
Area per protein: 2 nm^2
Area per lipid: 10.4716089904762 nm^2

Area per protein, upper half: 1.75 nm^2
Area per lipid, upper leaflet : 10.475577244 nm^2

Area per protein, lower half: 2 nm^2
Area per lipid, lower leaflet : 10.4716089904762 nm^2

Writing Area per lipid...
Done!

bash-3.2$ mdrun -v -deffnm em
Back Off! I just backed up em.log to ./#em.log.6#
Getting Loaded...
Reading file em.tpr, VERSION 4.5.4 (single precision)
Starting 4 threads
Loaded with Money

Making 1D domain decomposition 4 x 1 x 1

Back Off! I just backed up em.trr to ./#em.trr.6#

Back Off! I just backed up em.edr to ./#em.edr.6#

Steepest Descents:
Tolerance (Fmax)   =  1.0e+03
Number of steps=5
Step=0, Dmax= 1.0e-02 nm, Epot= -3.06212e+05 Fmax= 5.37269e+03, atom= 3275
Step=1, Dmax= 1.0e-02 nm, Epot= -3.10498e+05 Fmax= 2.40668e+03, atom= 3024
Step=3, Dmax= 6.0e-03 nm, Epot= -3.10958e+05 Fmax= 2.89119e+03, atom= 4225
Step=5, Dmax= 3.6e-03 nm, Epot= -3.12589e+05 Fmax= 1.46767e+03, atom= 1365
Step=6, Dmax= 4.3e-03 nm, Epot= -3.13033e+05 Fmax= 3.83969e+03, atom= 1365
Step=7, Dmax= 5.2e-03 nm, Epot= -3.13599e+05 Fmax= 2.88720e+03, atom= 1365
Step=8, Dmax= 6.2e-03 nm, Epot= -3.13688e+05 Fmax= 4.93154e+03, atom= 1365
Step=9, Dmax= 7.5e-03 nm, Epot= -3.14101e+05 Fmax= 4.72541e+03, atom= 1365
Step=   11, Dmax= 4.5e-03 nm, Epot= -3.14800e+05 Fmax= 9.90227e+02, atom= 1365

writing lowest energy coordinates.

Back Off! I just backed up em.gro to ./#em.gro.6#

Steepest Descents converged to Fmax<  1000 in 12 steps
Potential Energy  = -3.1480012e+05
Maximum force =  9.9022711e+02 on atom 1365
Norm of force =  1.2180042e+02


bash-3.2$ perl inflategro.pl em.gro 0.95 DPPC 0 system_shrink1.gro 5 
area_shrink1.dat
Reading.
Scaling lipids
There are 128 lipids...


Something is wrong here.  You should only have 126 lipids since 2 were removed 
before.  Did you minimize the correct coordinate file?



with 50 atoms per lipid..

Determining upper and lower leaflet...
64 lipids in the upper...
64 lipids in the lower leaflet

No protein coordinates found...


This is also strange.  This, coupled with the lines above regarding the number 
of lipids, suggest you're using the wrong coordinate file.



Writing scaled bilayer&  centered protein...


Calculating Area per lipid...
Protein X-min/max: 1-
Protein Y-min/max: 1-
X-range: -1 AY-range: -1 A
Building -1 X -1 2D grid on protein coordinates...
Calculating area occupied by protein..
full TMD..
upper TMD..
lower TMD..


The following areas are meaningless, for reasons associated with incorrect 
content or processing thereof.



Area per protein: 0 nm^2
Area per lipid: 0.583197761890625 nm^2

Area per protein, upper half: 0 nm^2
Area per lipid, upper leaflet : 0.583197761890625 nm^2

Area per protein, lower half: 0 nm^2
Area per lipid, lower leaflet : 0.583197761890625 nm^2

Writing Area per lipid...
Done!

To me, it looks like it is not looking at the 126 lipids that I updated in the 
topology file. H


OK, I seem to remember 124 being what my system produced, but that's irrelevant. 
 Check out the lines above.



Erica


From: Justin A. Lemkul [via GROMACS] [ml-node+s5086n4998244...@n6.nabble.com]
Sent: Friday, June 08, 2012 5:21 PM
To: Hicks, Erica
Subject: RE: Energy Minimization - not getting correct lipid area



On 6/8/12 6:18 PM, Erica Hicks wrote:

Hi,

Everything should be updated correctly. I added in the DPPC 128 underneath
the protein in [molecules] but I still get the same results. What do you mean
by " Further manual modifications are necessary based on how many lipids
InflateGRO removes"? After I generated the new position restrain file, I
updated the minim.mdp withe the correct info. I then scaled by a factor of 4,
updated [molecules], ran energy minimization (mdrun -v -deffnm em), and


This is all my comment meant.  If InflateGRO removes 4 lipids, then you have 124
left, which it sounds like you have accounted for.


scaled by a factor 0.95. I have no idea what could possibly going wrong.
Maybe something with

Re: [gmx-users] vsite as interaction site for COM of benzene

[Broadbent, Richard]

Dear Thomas,

Whilst in principle a constraint should do what you are asking, you might be
better off using multiple vsites to solve the problem so that there is an
algebraic as apposed to numerical mapping of the forces. My immediate
thought is to use two type 3 (or if necessary type 3fd) virtual sites to
define two points near the centre of the benzene. You could then place your
tabulated potential on a type 2 virtual site between these two. I believe
that should algebraically link all the atomic sites together making the
process more stable and it should be quite easy to transfer it to other
systems.

On a separate note why if you are coarse graining the benzenes to a single
point particle are you then maintaining all the atoms? Could you go a stage
further and map the atoms onto 3 non interacting point masses? That would
maintain angular momentum etc. whilst reducing the computational overhead.
Although it becomes very complex to implement if your benzene's are in the
backbone of a polymer.

Hope that helps,

Richard

On 08/06/2012 18:18, "Thomas Schlesier"  wrote:

> Hi all,
> 
> i have a more conceptional question, for using vsites as
> interaction-centers for coarse-grained particles:
> 
> First the simple case:
> I want to simulate one benzene molecule (atomistic - aa) in
> coarse-grained (cg-) benzene (each benzene molecule as a single
> particle). For the cg-cg interaction of the benzene i have calculated a
> table potential.
> 
> But for the aa-cg interaction I'm a little bit clueless:
> I want to use the COM of the aa-benzene as an interaction-site. But
> there is no vsite, which i could construct from 6 atoms.
> Side-note: The vsite would not interact with the with the atoms, only
> with cg-benzenes
> 
> So there are two ideas:
> 
> 1) Using one '3out'-vsite, which is construted from 3 non-neighbour
> atoms (if we had mesitylene instead of benzene, it would be the three
> C-atoms to which only hydrogens are bound). The the force from the aa-cg
> interaction would be distributed to these three atoms and i would hope
> that the bond-constraints would do 'the rest' (i.e. they would mimick
> that the whole aa-benzene molecule interacts with the cg-particles)
> 
> 2) Using two '3out'-vsites. First viste the same as in (1), the second
> comes from the set of the other three C-atoms. These two vsites, should
> be at the same position (more or less, for a perfect energy-minimized
> benzene molecule it would be so; but i would assume the error is only
> minimal). Since the potential is pair-additive, i would use for this
> table potential only half of the potential, so that the potential of
> both particles would add up to the 'true' potential. The forces would be
> calculated from the 'half-potential'.
> 
> Big questions is now:
> Are (1) and (2) equivalent?
> (2) Seems like the 'right way' to do.
> But (1) would be easier to implement :) But the big question is, if the
> bond-constraints would correct the approximation, that only 3 out of 6
> atoms have a direct interaction with the cg-particles.
> 
> Side note:
> The real system is a little more complex. There are fracments of
> molecules where i could not really construct two vsites which would have
> more or less the same position. In this case, method (1) would really help.
> 
> Any ideas?
> Greetings
> Thomas

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

[gmx-users] RE: Energy Minimization - not getting correct lipid area

[Erica Hicks]

Hi,

bash-3.2$ perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat
Reading.
Scaling lipids
There are 128 lipids...
with 50 atoms per lipid..

Determining upper and lower leaflet...
64 lipids in the upper...
64 lipids in the lower leaflet

Centering protein
Checking for overlap
...this might actually take a while
100 % done...
There are 2 lipids within cut-off range...
1 will be removed from the upper leaflet...
1 will be removed from the lower leaflet...

Writing scaled bilayer & centered protein...


Calculating Area per lipid...
Protein X-min/max: 2640
Protein Y-min/max: 2539
X-range: 14 AY-range: 14 A
Building 14 X 14 2D grid on protein coordinates...
Calculating area occupied by protein..
full TMD..
upper TMD
lower TMD
Area per protein: 2 nm^2
Area per lipid: 10.4716089904762 nm^2

Area per protein, upper half: 1.75 nm^2
Area per lipid, upper leaflet : 10.475577244 nm^2

Area per protein, lower half: 2 nm^2
Area per lipid, lower leaflet : 10.4716089904762 nm^2

Writing Area per lipid...
Done!

bash-3.2$ mdrun -v -deffnm em
Back Off! I just backed up em.log to ./#em.log.6#
Getting Loaded...
Reading file em.tpr, VERSION 4.5.4 (single precision)
Starting 4 threads
Loaded with Money

Making 1D domain decomposition 4 x 1 x 1

Back Off! I just backed up em.trr to ./#em.trr.6#

Back Off! I just backed up em.edr to ./#em.edr.6#

Steepest Descents:
   Tolerance (Fmax)   =  1.0e+03
   Number of steps=5
Step=0, Dmax= 1.0e-02 nm, Epot= -3.06212e+05 Fmax= 5.37269e+03, atom= 3275
Step=1, Dmax= 1.0e-02 nm, Epot= -3.10498e+05 Fmax= 2.40668e+03, atom= 3024
Step=3, Dmax= 6.0e-03 nm, Epot= -3.10958e+05 Fmax= 2.89119e+03, atom= 4225
Step=5, Dmax= 3.6e-03 nm, Epot= -3.12589e+05 Fmax= 1.46767e+03, atom= 1365
Step=6, Dmax= 4.3e-03 nm, Epot= -3.13033e+05 Fmax= 3.83969e+03, atom= 1365
Step=7, Dmax= 5.2e-03 nm, Epot= -3.13599e+05 Fmax= 2.88720e+03, atom= 1365
Step=8, Dmax= 6.2e-03 nm, Epot= -3.13688e+05 Fmax= 4.93154e+03, atom= 1365
Step=9, Dmax= 7.5e-03 nm, Epot= -3.14101e+05 Fmax= 4.72541e+03, atom= 1365
Step=   11, Dmax= 4.5e-03 nm, Epot= -3.14800e+05 Fmax= 9.90227e+02, atom= 1365

writing lowest energy coordinates.

Back Off! I just backed up em.gro to ./#em.gro.6#

Steepest Descents converged to Fmax < 1000 in 12 steps
Potential Energy  = -3.1480012e+05
Maximum force =  9.9022711e+02 on atom 1365
Norm of force =  1.2180042e+02


bash-3.2$ perl inflategro.pl em.gro 0.95 DPPC 0 system_shrink1.gro 5 
area_shrink1.dat
Reading.
Scaling lipids
There are 128 lipids...
with 50 atoms per lipid..

Determining upper and lower leaflet...
64 lipids in the upper...
64 lipids in the lower leaflet

No protein coordinates found...
Writing scaled bilayer & centered protein...


Calculating Area per lipid...
Protein X-min/max: 1-
Protein Y-min/max: 1-
X-range: -1 AY-range: -1 A
Building -1 X -1 2D grid on protein coordinates...
Calculating area occupied by protein..
full TMD..
upper TMD..
lower TMD..
Area per protein: 0 nm^2
Area per lipid: 0.583197761890625 nm^2

Area per protein, upper half: 0 nm^2
Area per lipid, upper leaflet : 0.583197761890625 nm^2

Area per protein, lower half: 0 nm^2
Area per lipid, lower leaflet : 0.583197761890625 nm^2

Writing Area per lipid...
Done!

To me, it looks like it is not looking at the 126 lipids that I updated in the 
topology file. H
Erica


From: Justin A. Lemkul [via GROMACS] [ml-node+s5086n4998244...@n6.nabble.com]
Sent: Friday, June 08, 2012 5:21 PM
To: Hicks, Erica
Subject: RE: Energy Minimization - not getting correct lipid area



On 6/8/12 6:18 PM, Erica Hicks wrote:
> Hi,
>
> Everything should be updated correctly. I added in the DPPC 128 underneath
> the protein in [molecules] but I still get the same results. What do you mean
> by " Further manual modifications are necessary based on how many lipids
> InflateGRO removes"? After I generated the new position restrain file, I
> updated the minim.mdp withe the correct info. I then scaled by a factor of 4,
> updated [molecules], ran energy minimization (mdrun -v -deffnm em), and

This is all my comment meant.  If InflateGRO removes 4 lipids, then you have 124
left, which it sounds like you have accounted for.

> scaled by a factor 0.95. I have no idea what could possibly going wrong.
> Maybe something with the manipulation of InflateGRO

Please provide the entire screen output of the first shrinking step that
InflateGRO produces.

-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing list[hidden email]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archi

Re: [gmx-users] RE: Energy Minimization - not getting correct lipid area

[Justin A. Lemkul]

On 6/8/12 6:18 PM, Erica Hicks wrote:

Hi,

Everything should be updated correctly. I added in the DPPC 128 underneath
the protein in [molecules] but I still get the same results. What do you mean
by " Further manual modifications are necessary based on how many lipids
InflateGRO removes"? After I generated the new position restrain file, I
updated the minim.mdp withe the correct info. I then scaled by a factor of 4,
updated [molecules], ran energy minimization (mdrun -v -deffnm em), and


This is all my comment meant.  If InflateGRO removes 4 lipids, then you have 124 
left, which it sounds like you have accounted for.



scaled by a factor 0.95. I have no idea what could possibly going wrong.
Maybe something with the manipulation of InflateGRO


Please provide the entire screen output of the first shrinking step that 
InflateGRO produces.


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

[gmx-users] RE: Energy Minimization - not getting correct lipid area

[Erica Hicks]

Hi,

Everything should be updated correctly. I added in the DPPC 128 underneath the 
protein in [molecules] but I still get the same results. What do you mean by " 
Further manual modifications are necessary based on how many lipids InflateGRO 
removes"? After I generated the new position restrain file, I updated the 
minim.mdp withe the correct info. I then scaled by a factor of 4, updated 
[molecules], ran energy minimization (mdrun -v -deffnm em), and scaled by a 
factor 0.95. I have no idea what could possibly going wrong. Maybe something 
with the manipulation of InflateGRO

Erica

From: Justin A. Lemkul [via GROMACS] [ml-node+s5086n4998242...@n6.nabble.com]
Sent: Friday, June 08, 2012 4:41 PM
To: Hicks, Erica
Subject: RE: Energy Minimization - not getting correct lipid area



On 6/8/12 5:32 PM, Erica Hicks wrote:
> Hi,
>
> I went back and found that a possible error could be that the coordinate file 
> was not updated correctly. After concatenated the Kalp_newbox.gro and 
> dppc128_whole.gro, it said to update the coordinate file with the correct 
> number of atoms. This should be the topology.top, right? But when I look 
> through this, all I find is:
> [ molecules ]
> ; Compound#mols
> Protein 1
> Where should the number of atoms be and is topology.top the correct file?

Make sure you added the correct topology information (step 2 in the tutorial,
note the green text):

http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/membrane_protein/02_topology.html

Step 3.2 tells you:

"Remove unnecessary lines (the box vectors from the KALP structure, the header
information from the DPPC structure) and update the second line of the
coordinate file (total number of atoms) accordingly."

The tutorial assumes you're familiar with such operations, so in reality what
it's implying is that, since you've now concatenated a protein and a 128-lipid
membrane, you need to add:

DPPC 128

to the [molecules] directive (after the Protein line to keep the correct order).
  Further manual modifications are necessary based on how many lipids InflateGRO
removes.

-Justin

>
> Erica
> 
> From: Justin A. Lemkul [via GROMACS] [[hidden email]]
> Sent: Friday, June 08, 2012 4:02 PM
> To: Hicks, Erica
> Subject: Re: Energy Minimization - not getting correct lipid area
>
>
>
> On 6/8/12 4:57 PM, Erica Hicks wrote:
>> It was after just one iteration. The number of lipids removed after scaling
>> the lipids by a factor of 4 was two. I am not sure how to update this as the
>> [molecules] directive in the topology.top shows only the number of moles,
>> i.e. 1.
>>
>
> The number shown in [molecules] is the number of molecules of a given species.
> If you started with 128 lipids, a removal of 4 means you need to list 124 
> lipids.
>
> I don't know how you can get such a small area with only one iteration of
> shrinking.  Your box should be enormous at that point.
>
> -Justin
>
> --
> 
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
> --
> gmx-users mailing list[hidden email]
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to [hidden email].
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>
> 
> If you reply to this email, your message will be added to the discussion 
> below:
>
> NAML
>
>
> --
> View this message in context: 
> http://gromacs.5086.n6.nabble.com/Energy-Minimization-not-getting-correct-lipid-area-tp4998235p4998241.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

Re: [gmx-users] RE: Energy Minimization - not getting correct lipid area

[Justin A. Lemkul]

On 6/8/12 5:32 PM, Erica Hicks wrote:

Hi,

I went back and found that a possible error could be that the coordinate file 
was not updated correctly. After concatenated the Kalp_newbox.gro and 
dppc128_whole.gro, it said to update the coordinate file with the correct 
number of atoms. This should be the topology.top, right? But when I look 
through this, all I find is:
[ molecules ]
; Compound#mols
Protein 1
Where should the number of atoms be and is topology.top the correct file?


Make sure you added the correct topology information (step 2 in the tutorial, 
note the green text):


http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/membrane_protein/02_topology.html

Step 3.2 tells you:

"Remove unnecessary lines (the box vectors from the KALP structure, the header 
information from the DPPC structure) and update the second line of the 
coordinate file (total number of atoms) accordingly."


The tutorial assumes you're familiar with such operations, so in reality what 
it's implying is that, since you've now concatenated a protein and a 128-lipid 
membrane, you need to add:


DPPC 128

to the [molecules] directive (after the Protein line to keep the correct order). 
 Further manual modifications are necessary based on how many lipids InflateGRO 
removes.


-Justin



Erica

From: Justin A. Lemkul [via GROMACS] [ml-node+s5086n4998239...@n6.nabble.com]
Sent: Friday, June 08, 2012 4:02 PM
To: Hicks, Erica
Subject: Re: Energy Minimization - not getting correct lipid area



On 6/8/12 4:57 PM, Erica Hicks wrote:

It was after just one iteration. The number of lipids removed after scaling
the lipids by a factor of 4 was two. I am not sure how to update this as the
[molecules] directive in the topology.top shows only the number of moles,
i.e. 1.



The number shown in [molecules] is the number of molecules of a given species.
If you started with 128 lipids, a removal of 4 means you need to list 124 
lipids.

I don't know how you can get such a small area with only one iteration of
shrinking.  Your box should be enormous at that point.

-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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



If you reply to this email, your message will be added to the discussion below:
http://gromacs.5086.n6.nabble.com/Energy-Minimization-not-getting-correct-lipid-area-tp4998235p4998239.html
To unsubscribe from Energy Minimization - not getting correct lipid area, click 
here.
NAML


--
View this message in context: 
http://gromacs.5086.n6.nabble.com/Energy-Minimization-not-getting-correct-lipid-area-tp4998235p4998241.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.


--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

[gmx-users] RE: Energy Minimization - not getting correct lipid area

[Erica Hicks]

Hi,

I went back and found that a possible error could be that the coordinate file 
was not updated correctly. After concatenated the Kalp_newbox.gro and 
dppc128_whole.gro, it said to update the coordinate file with the correct 
number of atoms. This should be the topology.top, right? But when I look 
through this, all I find is:
[ molecules ]
; Compound#mols
Protein 1
Where should the number of atoms be and is topology.top the correct file?

Erica

From: Justin A. Lemkul [via GROMACS] [ml-node+s5086n4998239...@n6.nabble.com]
Sent: Friday, June 08, 2012 4:02 PM
To: Hicks, Erica
Subject: Re: Energy Minimization - not getting correct lipid area



On 6/8/12 4:57 PM, Erica Hicks wrote:
> It was after just one iteration. The number of lipids removed after scaling
> the lipids by a factor of 4 was two. I am not sure how to update this as the
> [molecules] directive in the topology.top shows only the number of moles,
> i.e. 1.
>

The number shown in [molecules] is the number of molecules of a given species.
If you started with 128 lipids, a removal of 4 means you need to list 124 
lipids.

I don't know how you can get such a small area with only one iteration of
shrinking.  Your box should be enormous at that point.

-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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



If you reply to this email, your message will be added to the discussion below:
http://gromacs.5086.n6.nabble.com/Energy-Minimization-not-getting-correct-lipid-area-tp4998235p4998239.html
To unsubscribe from Energy Minimization - not getting correct lipid area, click 
here.
NAML


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

[gmx-users] floating point exception in .xtc file

[Christopher Neale]

Dear Users:

I have a 500 ns trajectory of 65 GB that gives a floating point exception when 
I run it through gmxcheck or trjcat (generated and analyzed with gromacs 
4.5.5). Has anybody encountered this? I ran mdrun with -append so this is the 
xtc resulting from months of simulation of a 1,000,000 atom system. If I run 
trjconv -f md.xtc -b 20, where the floating point exception occurred around 
t=18 ps in gmxcheck, then I can extract the readable frames and repair 
around the damaged section. Still, I'd rather not lose any data and I had 
thought that the new default -append option to mdrun checked for these types of 
problems at runtime.

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

Re: [gmx-users] Re: Energy Minimization - not getting correct lipid area

[Justin A. Lemkul]

On 6/8/12 4:57 PM, Erica Hicks wrote:

It was after just one iteration. The number of lipids removed after scaling
the lipids by a factor of 4 was two. I am not sure how to update this as the
[molecules] directive in the topology.top shows only the number of moles,
i.e. 1.



The number shown in [molecules] is the number of molecules of a given species. 
If you started with 128 lipids, a removal of 4 means you need to list 124 lipids.


I don't know how you can get such a small area with only one iteration of 
shrinking.  Your box should be enormous at that point.


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

[gmx-users] Re: Energy Minimization - not getting correct lipid area

[Erica Hicks]

It was after just one iteration. The number of lipids removed after scaling
the lipids by a factor of 4 was two. I am not sure how to update this as the
[molecules] directive in the topology.top shows only the number of moles,
i.e. 1.

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

[gmx-users] CECAM ScalaLife Workshop, Oct 03 - Oct 05, 2012

[Rossen Apostolov]

Dear GROMACS Users,

GROMACS is one of the applications participating in the ScalaLife EU 
project (www.scalalife.eu).


ScalaLife is organising a CECAM workshop on "High Performance Computing 
in Computational Chemistry and Molecular Biology: Challenges and 
Solutions provided by ScalaLife project", between October 3 and 5, in 
CECAM-HQ-EPFL, Lausanne, Switzerland.


ScalaLife has invited renowned life scientists from Europe and USA to 
discuss recent advances in large scale simulations and showcase their 
latest results. The workshop will also present the newest developments 
in three life science software packages: GROMACS, DALTON and DISCRETE.


The workshop features multiple hands-on sessions, giving the 
participants the opportunity to practice the modelling capabilities of 
the three software packages. These sessions will also include 
user-clinics, a unique opportunity for users to bring in their own 
problems and discuss the challenges they face in modelling with the 
developers of these packages.


Registration:
=
Life-scientist and researchers from academia and industry are cordially 
invited to participate in the workshop.


To apply for participation please go to 
http://www.cecam.org/workshop-6-672.html or send an email to 
scalal...@pdc.kth.se


* The number of participants is limited. There is no attendance fee. *
* Deadline for registration is 1st of September *

Coffee breaks, lunches and the workshop dinner listed in the program are 
covered by the workshop.


Agenda:
===
The detailed program can be downloaded at 
http://www.scalalife.eu/content/cecam.
The short version can be found online at 
http://www.cecam.org/workshop-2-672.html


Location:
=
A map of EPFL: 
http://www.cecam.org/upload/misc/EPFL-plan-campus-lausanne.jpg.


For more information on how to get to CECAM Headquarters: 
http://www.cecam.org/lausanne2011.html.


Information about nearby hotels will be sent with the confirmation email 
for your registration.


More information about the workshop: 
http://www.cecam.org/workshop-0-672.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

Re: [gmx-users] Energy Minimization - not getting correct lipid area

[Justin A. Lemkul]

On 6/8/12 3:15 PM, Hicks, Erica wrote:

Hi,

I am working through the KALP-15 tutorial and having difficulties getting the
correct area per lipid (~71 A^2). I first scaled the lipid by a factor of 4
(perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat) then
did energy minimization (mdrun -v -deffnm em) and then scaled down by 0.95 (perl
inflategro.pl em.gro 0.95 DPPC 0 system_shrink.gro 5 area_shrink1.dat). The
tutorial then said to do energy minimizations and scale by 0.95 until the
correct area per lipid is reached. When I do, my numbers stay the same, i.e. all
three are 0.583197761890625 nm^2. Any suggestions?



After how many iterations?  If you've hit 0.58 nm^2, you've gone too far. 
You'll likely hit a point where you can't squeeze the lipids any more, but 25-26 
iterations consistently gives the value stated in the tutorial.


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

[gmx-users] Energy Minimization - not getting correct lipid area

[Hicks, Erica]

Hi,

I am working through the KALP-15 tutorial and having difficulties getting the 
correct area per lipid (~71 A^2). I first scaled the lipid by a factor of 4 
(perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat) then 
did energy minimization (mdrun -v -deffnm em) and then scaled down by 0.95 
(perl inflategro.pl em.gro 0.95 DPPC 0 system_shrink.gro 5 area_shrink1.dat). 
The tutorial then said to do energy minimizations and scale by 0.95 until the 
correct area per lipid is reached. When I do, my numbers stay the same, i.e. 
all three are 0.583197761890625 nm^2. Any suggestions?


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

Re: [gmx-users] parameters for bond types for GROMOS force field.

[James Starlight]

Justin,



thanks alot. I'll try to use gb_10 ( for C=N) as well as gb_16 ( for C=C)
for parametrisation in my task.


James

2012/6/8 Justin A. Lemkul 

>
>
> On 6/8/12 2:01 PM, James Starlight wrote:
>
>> Justin,
>>
>> Does the gb_10 suitable for both bonds that I want to parametrise ?
>>
>> As I've told previously I need to parametrise two different bonds
>>
>> 1) is the >C=C< bond wich are not in the ring. This is just double bound
>> in
>> linnear sequence -C=C-.
>>
>> 2) is the >C=N- double bond where both atoms in the 5-m ring system
>>
>>
> This is a better description.  The last message suggested everything was
> in a 5-membered ring.  Had that been the case, my suggestion would have
> been right (see, for instance, a HIS side chain - gb_10 is used for both
> C=C and C=N in that case).
>
> There are no linear C=C bonds in Gromos force fields by default, as such
> groups are not in proteins.  The closest you can likely come is to take a
> C=C bond from the PHE ring.  Bond lengths are simple, and are usually
> derived from spectroscopic information.  If using constraints on all bonds,
> the force constant is irrelevant since the bond is rigid.  If not using
> constraints, then you'll have to find or derive suitable parameters
> yourself.  Welcome to the inconvenient world of parameterization :)
>
> -Justin
>
>  Thanks again,
>>
>> James
>>
>>
>>
>> 2012/6/8 Justin A. Lemkul mailto:jalem...@vt.edu>>
>>
>>
>>
>>
>>On 6/8/12 1:22 PM, James Starlight wrote:
>>
>>I've found that information in ffbonded.itp but I'm not sure about
>> exactly
>>meaning of some types.
>>
>>
>>E.g I'm looking for bond type for simple double bond between C=C
>> as well
>>as C=N
>>( here both atoms in 5-m ring) atoms pairs. Should I use for the
>> first one
>>#define gb_15   0.1390  8.6600e+06
>>; CH2  -  C, CR1 (6-ring)   800
>>
>>   and for the second one ;
>>
>>#define gb_11   0.1340  1.0500e+07
>>; C  -  N, NZ, NE   900
>>?
>>
>>
>>Neither.  gb_10 sounds exactly like what you need:
>>
>>#define gb_10   0.1330  1.1800e+07
>>; C, CR1  -  N, NR, CR1, C (peptide, 5-ring)   1000
>>
>>-Justin
>>
>>--
>>==**__==
>>
>>
>>Justin A. Lemkul, Ph.D.
>>Research Scientist
>>Department of Biochemistry
>>Virginia Tech
>>Blacksburg, VA
>>jalemkul[at]vt.edu  | (540) 231-9080
>>
>> http://www.bevanlab.biochem.__**vt.edu/Pages/Personal/justin
>>
>> 
>> >
>>
>>==**__==
>>
>>--
>>gmx-users mailing list gmx-users@gromacs.org > gmx-users@gromacs.org>
>>
>> http://lists.gromacs.org/__**mailman/listinfo/gmx-users
>>
>>
>> 
>> >
>>Please search the archive at
>>
>> http://www.gromacs.org/__**Support/Mailing_Lists/Search
>>
>>
>> >
>> before posting!
>>Please don't post (un)subscribe requests to the list. Use the www
>> interface
>>or send it to gmx-users-requ...@gromacs.org
>>> >.
>>Can't post? Read 
>> http://www.gromacs.org/__**Support/Mailing_Lists
>>
>> 
>> >
>>
>>
>>
> --
> ==**==
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin
>
> ==**==
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Searchbefore
>  posting!
> Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read 
> http://www.gromacs.org/**Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post

Re: [gmx-users] parameters for bond types for GROMOS force field.

[Justin A. Lemkul]

On 6/8/12 2:01 PM, James Starlight wrote:

Justin,

Does the gb_10 suitable for both bonds that I want to parametrise ?

As I've told previously I need to parametrise two different bonds

1) is the >C=C< bond wich are not in the ring. This is just double bound in
linnear sequence -C=C-.

2) is the >C=N- double bond where both atoms in the 5-m ring system



This is a better description.  The last message suggested everything was in a 
5-membered ring.  Had that been the case, my suggestion would have been right 
(see, for instance, a HIS side chain - gb_10 is used for both C=C and C=N in 
that case).


There are no linear C=C bonds in Gromos force fields by default, as such groups 
are not in proteins.  The closest you can likely come is to take a C=C bond from 
the PHE ring.  Bond lengths are simple, and are usually derived from 
spectroscopic information.  If using constraints on all bonds, the force 
constant is irrelevant since the bond is rigid.  If not using constraints, then 
you'll have to find or derive suitable parameters yourself.  Welcome to the 
inconvenient world of parameterization :)


-Justin


Thanks again,

James



2012/6/8 Justin A. Lemkul mailto:jalem...@vt.edu>>



On 6/8/12 1:22 PM, James Starlight wrote:

I've found that information in ffbonded.itp but I'm not sure about 
exactly
meaning of some types.


E.g I'm looking for bond type for simple double bond between C=C as well
as C=N
( here both atoms in 5-m ring) atoms pairs. Should I use for the first 
one
#define gb_15   0.1390  8.6600e+06
; CH2  -  C, CR1 (6-ring)   800

   and for the second one ;

#define gb_11   0.1340  1.0500e+07
; C  -  N, NZ, NE   900
?


Neither.  gb_10 sounds exactly like what you need:

#define gb_10   0.1330  1.1800e+07
; C, CR1  -  N, NR, CR1, C (peptide, 5-ring)   1000

-Justin

--
==__==

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu  | (540) 231-9080
http://www.bevanlab.biochem.__vt.edu/Pages/Personal/justin


==__==

--
gmx-users mailing list gmx-users@gromacs.org 
http://lists.gromacs.org/__mailman/listinfo/gmx-users

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





--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

Re: [gmx-users] parameters for bond types for GROMOS force field.

[James Starlight]

Justin,

Does the gb_10 suitable for both bonds that I want to parametrise ?

As I've told previously I need to parametrise two different bonds

1) is the >C=C< bond wich are not in the ring. This is just double bound in
linnear sequence -C=C-.

2) is the >C=N- double bond where both atoms in the 5-m ring system

Thanks again,

James



2012/6/8 Justin A. Lemkul 

>
>
> On 6/8/12 1:22 PM, James Starlight wrote:
>
>> I've found that information in ffbonded.itp but I'm not sure about exactly
>> meaning of some types.
>>
>>
>> E.g I'm looking for bond type for simple double bond between C=C as well
>> as C=N
>> ( here both atoms in 5-m ring) atoms pairs. Should I use for the first one
>> #define gb_15   0.1390  8.6600e+06
>> ; CH2  -  C, CR1 (6-ring)   800
>>
>>   and for the second one ;
>>
>> #define gb_11   0.1340  1.0500e+07
>> ; C  -  N, NZ, NE   900
>> ?
>>
>>
> Neither.  gb_10 sounds exactly like what you need:
>
> #define gb_10   0.1330  1.1800e+07
> ; C, CR1  -  N, NR, CR1, C (peptide, 5-ring)   1000
>
> -Justin
>
> --
> ==**==
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin
>
> ==**==
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Searchbefore
>  posting!
> Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read 
> http://www.gromacs.org/**Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] parameters for bond types for GROMOS force field.

[Justin A. Lemkul]

On 6/8/12 1:22 PM, James Starlight wrote:

I've found that information in ffbonded.itp but I'm not sure about exactly
meaning of some types.


E.g I'm looking for bond type for simple double bond between C=C as well as C=N
( here both atoms in 5-m ring) atoms pairs. Should I use for the first one
#define gb_15   0.1390  8.6600e+06
; CH2  -  C, CR1 (6-ring)   800

   and for the second one ;

#define gb_11   0.1340  1.0500e+07
; C  -  N, NZ, NE   900
?



Neither.  gb_10 sounds exactly like what you need:

#define gb_10   0.1330  1.1800e+07
; C, CR1  -  N, NR, CR1, C (peptide, 5-ring)   1000

-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

Re: [gmx-users] parameters for bond types for GROMOS force field.

[James Starlight]

I've found that information in ffbonded.itp but I'm not sure about exactly
meaning of some types.


E.g I'm looking for bond type for simple double bond between C=C as well as
C=N ( here both atoms in 5-m ring) atoms pairs. Should I use for the first
one
#define gb_15   0.1390  8.6600e+06
; CH2  -  C, CR1 (6-ring)   800

  and for the second one ;

#define gb_11   0.1340  1.0500e+07
; C  -  N, NZ, NE   900
?

Thanks for help

James

2012/6/8 Mark Abraham 

> On 9/06/2012 12:00 AM, James Starlight wrote:
>
>> Dear Gromacs Users!
>>
>>
>> I'm looking for description of the parameters of bonds terms ( termed as
>> the gb_# in the topology.top file) . Could you tell me where I could find
>> such descriptions for all possible bond types ?
>>
>
> See manual for literature references for the force field components.
>
> 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/Searchbefore
>  posting!
> Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read 
> http://www.gromacs.org/**Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] vsite as interaction site for COM of benzene

[Thomas Schlesier]

Hi all,

i have a more conceptional question, for using vsites as 
interaction-centers for coarse-grained particles:


First the simple case:
I want to simulate one benzene molecule (atomistic - aa) in 
coarse-grained (cg-) benzene (each benzene molecule as a single 
particle). For the cg-cg interaction of the benzene i have calculated a 
table potential.


But for the aa-cg interaction I'm a little bit clueless:
I want to use the COM of the aa-benzene as an interaction-site. But 
there is no vsite, which i could construct from 6 atoms.
Side-note: The vsite would not interact with the with the atoms, only 
with cg-benzenes


So there are two ideas:

1) Using one '3out'-vsite, which is construted from 3 non-neighbour 
atoms (if we had mesitylene instead of benzene, it would be the three 
C-atoms to which only hydrogens are bound). The the force from the aa-cg 
interaction would be distributed to these three atoms and i would hope 
that the bond-constraints would do 'the rest' (i.e. they would mimick 
that the whole aa-benzene molecule interacts with the cg-particles)


2) Using two '3out'-vsites. First viste the same as in (1), the second 
comes from the set of the other three C-atoms. These two vsites, should 
be at the same position (more or less, for a perfect energy-minimized 
benzene molecule it would be so; but i would assume the error is only 
minimal). Since the potential is pair-additive, i would use for this 
table potential only half of the potential, so that the potential of 
both particles would add up to the 'true' potential. The forces would be 
calculated from the 'half-potential'.


Big questions is now:
Are (1) and (2) equivalent?
(2) Seems like the 'right way' to do.
But (1) would be easier to implement :) But the big question is, if the 
bond-constraints would correct the approximation, that only 3 out of 6 
atoms have a direct interaction with the cg-particles.


Side note:
The real system is a little more complex. There are fracments of 
molecules where i could not really construct two vsites which would have 
more or less the same position. In this case, method (1) would really help.


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

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

[gmx-users] Re: Re: energy conservation: shift vs shifted user potential

[Anja Kuhnhold]

Ok, I try to explain it again :)
I know the format of the user tables (x, f(x), -f'(x), g(x), -g'(x), h(x), 
-h'(x);
f is coulomb, g dispersion, h repulsion).
For testing the simple cut-off version, my table looks like the
example in the installation.
So I run a simulation with vdwtype=cut-off
and the same with vdwtype=user.
The results are the same.

The problem is, that I want to use shifted forces and potentials within my 
tables.
(Aim is better energy conservation)
In the manual (ch. 4.1.5) is explained how the shift functions
for forces and potentials have to be.
So for -g(x) and h(x) I use the potential from the manual, Phi, with alpha=6 
and 12, respectivly.
And for g'(x) and -h'(x) I use the shifted force from the manual, Fs, with 
alpha=6 and 12, respectivly.
But when I run the simulation with this table on the one hand
and with vdwtype=Shift on the other hand,
the results are not comparable, mostly the energy is not conserved by using the 
table.

Can you follow me?

Anja



- Ursprüngliche Nachricht -
Von: gmx-users-requ...@gromacs.org
Datum: Freitag, Juni 8, 2012 15:55
Betreff: gmx-users Digest, Vol 98, Issue 52
An: gmx-users@gromacs.org

> - Ursprüngliche Nachricht -
> Von Mark Abraham 
> Datum Fri, 08 Jun 2012 23:40:35 +1000
> An Discussion list for GROMACS users 
> Betreff Re: [gmx-users] Re: energy conservation: shift vs shifted user 
> potential
> On 8/06/2012 5:59 PM, Anja Kuhnhold wrote:
> > Hello all,
> >
> > can someone give me a hint, please?
> > Do you need more information?
> > Has anyone had a similar problem.
> >
> > I really need to figure that out, because I will simulate other systems
> > which can be run only with user tables.
> >
> > Anja
> >
> >
> >
> > -
> >> - Original Message -
> >> - Original Message -
> >> From Anja Kuhnhold
> >> Date Wed, 06 Jun 2012 15:07:40 +0200
> >> To gmx-users@gromacs.org
> >> Subject [gmx-users] energy conservation: shift vs shifted user potential
> >> Dear gmx-users,
> >>
> >> I have a problem concerning energy conservation when using user
> >> potentials (tables).
> >> I'm using gromacs 4.5.4
> >> I simulate a simple Lennard-Jones(6-12) +Fene polymer melt (1600
> >> chains a 10 beads in a 26x26x26 periodic box).
> >>
> >> I tried different vdwtypes (cutoff always 3.24):
> >> The cut-off version does not conserve energy -- okay.
> >> The shift and switch versions conserve energy -- fine.
> >>
> >> Now I wanted to do the same with user tables:
> >> Simple Lennard-Jones table gives really the same results as the
> >> cut-off version -- good.
> >>
> >> But if I use a table with shifted Lennard-Jones potential it is not
> >> comparable to the shift version
> >> and the energy is not conserved -- ?
> >>
> >> I use a shift function as written in the manual (chapter 4.1.5) --
> >> there must be a factor alpha added in the constants A and B --
> >> (r1=0).
> >>
> >> The parameters are the same for shift version and shifted user version.
> >>
> >> Has someone an idea why the shifted user potential doesn't work in
> >> this way?
> 
> We've no real idea what you've done... manual 6.7.2 describes the 
> required format and there are (unshifted) examples in your 
> installation 
> under $GMXLIB/share/top/table*.xvg.
> 
> Makr
> 
> >>
> >> Here is the mdp:
> >>
> >>
> >> integrator = md-vv
> >> dt = 0.0035
> >> nsteps = 1000
> >> nstxout= 1
> >> nstvout= 1
> >> nstfout= 1
> >> nstlog = 1
> >> ns_type= grid
> >> pbc= xyz
> >> rvdw   = 3.24
> >> rlist  = 3.6
> >> tcoupl = no
> >> tc-grps= System
> >> tau_t  = 2.0
> >> ref_t  = 127.2717
> >> vdwtype= user;Shift
> >> rcoulomb   = 3.6;2.24;1.12
> >> coulombtype= Cut-off
> >> rvdw-switch= 0.0
> >>
> >> energygrps = bead
> >> energygrp_table= bead bead
> >>
> >>
> >> Thanks in advance
> >> Anja
> >>
> 

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

Re: [gmx-users] Electrostatic interaction energy

[Justin A. Lemkul]

On 6/8/12 11:04 AM, Turgay Cakmak wrote:

Hi all,

I have 2 questions related to the electrostatic energy of the system..

1) Is *total electrostatic energy* equal to *Coulomb-(SR)* + *Coul.-recip*. OR
*Coulomb-14* + *Coulomb-(SR)* + *Coul.-recip*. I am so confused.



The total electrostatic energy is the sum of all Coulombic components.  Whether 
that is relevant is up to you.



I specified in my .mdp file;
coulombtype = PME
rvdw = 1.0
rlist = 1.0
rcoulomb = 1.0
fourierspacing = 0.16
pme_order = 4
ewald_rtol = 1e-5

2) I want to calculate electrostatic interaction energy between 2 residues in my
system which has 10 peptides in total. But I think, I should have specified the
energy groups composed of that 2  residues in my mdp file before starting the
simulation. So in this case, is there any way to get this interaction energy OR
can I just calculate the electrostatic energy of the whole system?



Yes, to get an energy decomposition between two residues you need to specify 
appropriate energygrps.  You can create a new .tpr file and back-calculate these 
quantities using the mdrun -rerun feature.  Be aware though, that the PME 
reciprocal term cannot be (easily) decomposed and the best you can do is 
calculate short-range nonbonded energies between your chosen groups.  Whether or 
not that is useful is something you must demonstrate.


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


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

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

[gmx-users] Electrostatic interaction energy

[Turgay Cakmak]

Hi all,

I have 2 questions related to the electrostatic energy of the system..

1) Is *total electrostatic energy* equal to *Coulomb-(SR)* + *Coul.-recip*.
OR *Coulomb-14* + *Coulomb-(SR)* + *Coul.-recip*. I am so confused.

I specified in my .mdp file;
coulombtype = PME
rvdw = 1.0
rlist = 1.0
rcoulomb = 1.0
fourierspacing = 0.16
pme_order = 4
ewald_rtol = 1e-5

2) I want to calculate electrostatic interaction energy between 2 residues
in my system which has 10 peptides in total. But I think, I should have
specified the energy groups composed of that 2  residues in my mdp file
before starting the simulation. So in this case, is there any way to get
this interaction energy OR can I just calculate the electrostatic energy of
the whole system?

Thanks in advance,

Turgay



2012/6/8 Mark Abraham 

>  On 8/06/2012 11:51 PM, Turgay Cakmak wrote:
>
> Thanks Mark.
> I want to calculate Coulomb and Lennard-Jones energies using
> total_dt100.edr .
> But, I am not sure I can get the correct energies?
>
>
> Neither am I. Only you know the history of these .edr files.
>
> Mark
>
>
>
> Turgay
>
> 2012/6/8 Mark Abraham 
>
>>  On 8/06/2012 10:40 PM, Turgay Cakmak wrote:
>>
>>
>> Hi all,
>>
>> I have encountered a problem using eneconv with -dt flag.
>> I have several energy files. Firstly, I concatenated them using the
>> following.
>>
>> eneconv  -f  first.edr   second.edr    -settime  -o  total.edr
>>
>> But, total.edr is so large. To reduce the size of file, I used the below
>> command:
>>
>> eneconv -f  total.edr  -dt  100 -o  total_dt100.edr
>>
>> But, I get the following warning:
>>* .*
>>* .*
>>
>> WARNING: missing energy sums at time 69984.00
>>
>> WARNING: missing energy sums at time 69988.00
>>
>> WARNING: missing energy sums at time 69992.00
>>
>> WARNING: missing energy sums at time 69996.00
>> Reading energy frame  35000 time 7.000
>> WARNING: missing energy sums at time 7.00
>> Last energy frame read 35000 time 7.000
>> Last step written from energy.edr: t 7, step 3500
>>
>> Last frame written was at step at step 3500, time 7.00
>> Wrote 701 frames
>>
>>
>>  At some point something did something to break the usual format of the
>> files... but we can have no idea what.
>>
>>
>>
>> Eventually, I get the "total_dt100.edr". But, I am not sure I could
>> continue using this file.
>> If someone could help me, I would be glad.
>>
>>
>>  Depends what you want to use it for.
>>
>
>
>>
>> Mark
>>
>> --
>> gmx-users mailing listgmx-users@gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-requ...@gromacs.org.
>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>
>
>
>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] qmmm and "real-life" meaning of null lennard-jones parameters

[Javier Cerezo]

The system in you topology must be complete even for QMMM, so you'll 
need FF parameters also for the QM part.


Partition models for QM/MM work with the following energy decomposition:
E(total) = E(QM_system) + E(MM_system) + E(QM-MM interaction)

E(QM_system) and E(MM_system) are easy to define (computed at QM and MM 
levels respectively), but E(QM-MM interaction) is a bit more involved. 
One way to account for this interaction energy is through the ONIOM (and 
the like) partition scheme:
E(total) = E(QM_system@QM level) + E(whole_system@MM level) - 
E(QM_system@MM level)


So, you need to evaluate the whole system at MM level and that's why you 
need a force field for this part as well. Anyway, I think that for the 
QM_system tretaed at MM level only nonbonded interactions mater, since 
bonded interactions (that are not involved with the boundaries) are 
cancelled out when performing this subtraction: E(whole_system@MM level) 
- E(QM_system@MM level). For that reason, I think that bonded parameters 
that are not present in the boundaries, are not relevant and could be 
omitted in the topology.


Javier



El 08/06/12 01:47, Justin A. Lemkul escribió:



On 6/7/12 7:28 PM, Edward Deira wrote:

Dear all,

I'm currently starting to dwell deeper in MD, and I'm taking some 
time to

understand what's going on inside the gromacs "black-box".
In one of those dwellings, I came across an older post
[http://www.mail-archive.com/gmx-users@gromacs.org/msg42568.html] 
which reads:


Question:
4. In ffnonbonded.itp, why are both sigma and epsilon set to zero for HW
(opls_117)? This seems to imply that, as far as Lennard-Jones 
interactions are
concerned, the hydrogens on the waters don't exist. Or, in other 
words, in the
absence of charges, the hydrogens don't "feel" the hydrogens, the 
hydrogens
don't "feel" the oxygens, and the oxygens don't "feel" the hydrogens. 
In other
words, the hydrogens interact with the world only via electrostatic 
(Coulombic)
interactions. Is this a correct interpretation?Correct. Many force 
fields do this.

Answer:

So, my question, if a question at all:

Suppose I have a regular protein and put inside some metal atom that 
will
coordinate with some O and N atoms from the side chains. If the sigma 
and
epsilon for that metal are null, than the metal - sidechains 
interactions are
exclusively electrostatic. Does this make sense ? What are the 
implications of
this for the "coordination chemistry" of that "metal - sidechain 
complex" ?


On the side: suppose I want some non parameterized metal atom, say W, 
for which
I will compute all the other parameters in the same/similar way 
described in the
force field papers, but for which no experimental data are available 
for me to
compare computable meaningful sigma and epsilon values. Can I just 
sigma and

epsilon to zero ? Or should I do qmmm to have W in the qm part ?



The fact that the LJ parameters for H are zero derives from its size.  
The environment is more strongly influenced by the heavy atom to which 
H is bonded.  In the case of a larger metal ion, I would seriously 
doubt that setting LJ parameters to zero is valid.  It's quite 
convenient, but in most force fields, all metal ions have some LJ 
parameters.  Perhaps investigating how those parameters were derived 
would be useful.  For what it's worth, I believe the origin of the 
zero-LJ H parameters comes from this work:


http://pubs.acs.org/doi/abs/10.1021/ja00824a004

Also, from the few tutorials and from the manual, I have the 
impression that
even for qmmm with gromacs and mopac i still need force field 
parameters for the
qm part, is this true ? Or i just need to include a qmmm section in 
all mdp
files, including the first ion adding and energy minimization steps ? 
Sorry for

the naivety in this, but i've only made "regular protein" MD so far.



I've never done any QM/MM, but my assumption would be that you have to 
have some valid topology to start with.  Perhaps someone else can 
comment on this methodological issue.


-Justin



--
Javier CEREZO BASTIDA
PhD Student
Physical Chemistry
Universidad de Murcia
Murcia (Spain)
Tel: (+34)868887434
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] Total energy of a group of atoms

[Markus Kaukonen]

Dear All,

I'm planning to build a qm/mm interface using ase package
(https://wiki.fysik.dtu.dk/ase/ ).

Question:
Is it possible to divide the system in two (or more) groups say
group1: qm
group2 : mm
and get the total energy assoaciated to each of the groups?
(both bonded and non-bonded interactions).

terveisin, Markus


-- 
--www=http://www.iki.fi/markus.kaukonen
--markus.kauko...@iki.fi
--office: 1st year student at Helsinki University
--home: Viinirinne 3 F 12, 02630 Espoo, FIN
--tel: h 045-1242068
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] parameters for bond types for GROMOS force field.

[Mark Abraham]

On 9/06/2012 12:00 AM, James Starlight wrote:

Dear Gromacs Users!


I'm looking for description of the parameters of bonds terms ( termed 
as the gb_# in the topology.top file) . Could you tell me where I 
could find such descriptions for all possible bond types ?


See manual for literature references for the force field components.

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

Re: [gmx-users] Furthest distance between 2 atoms

[Mark Abraham]

On 8/06/2012 11:49 PM, Catarina Santos wrote:

Hi everybody,

I am trying to find an easy way to calculate the furthest distance 
between the atoms of two groups in a system. I know that g_dist can 
give me all the distances I want, but I would have to run g_dist for 
each possible atom pair (which is not a fast or a straightforward 
process).


I have already checked the manual and tried some tools, I did not 
find anything that returns all the distances between all the possible 
pairs of atoms though. Am I missing something?




g_mindist -max with both groups the same?

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] parameters for bond types for GROMOS force field.

[James Starlight]

Dear Gromacs Users!


I'm looking for description of the parameters of bonds terms ( termed as
the gb_# in the topology.top file) . Could you tell me where I could find
such descriptions for all possible bond types ?



thanks,

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

Re: [gmx-users] eneconv with -dt flag

[Mark Abraham]

On 8/06/2012 11:51 PM, Turgay Cakmak wrote:

Thanks Mark.
I want to calculate Coulomb and Lennard-Jones energies using 
total_dt100.edr .

But, I am not sure I can get the correct energies?


Neither am I. Only you know the history of these .edr files.

Mark



Turgay

2012/6/8 Mark Abraham >


On 8/06/2012 10:40 PM, Turgay Cakmak wrote:


Hi all,

I have encountered a problem using eneconv with -dt flag.
I have several energy files. Firstly, I concatenated them using
the following.

eneconv-f first.edr   second.edr    -settime -o total.edr

But, total.edr is so large. To reduce the size of file, I used
the below command:

eneconv -f  total.edr  -dt  100 -o  total_dt100.edr

But, I get the following warning:
*.*
*.*

WARNING: missing energy sums at time 69984.00

WARNING: missing energy sums at time 69988.00

WARNING: missing energy sums at time 69992.00

WARNING: missing energy sums at time 69996.00
Reading energy frame  35000 time 7.000
WARNING: missing energy sums at time 7.00
Last energy frame read 35000 time 7.000
Last step written from energy.edr: t 7, step 3500

Last frame written was at step at step 3500, time 7.00
Wrote 701 frames


At some point something did something to break the usual format of
the files... but we can have no idea what.




Eventually, I get the "total_dt100.edr". But, I am not sure I
could continue using this file.
If someone could help me, I would be glad.


Depends what you want to use it for.


Mark

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

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






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

Re: [gmx-users] eneconv with -dt flag

[Turgay Cakmak]

Thanks Mark.
I want to calculate Coulomb and Lennard-Jones energies using
total_dt100.edr .
But, I am not sure I can get the correct energies?

Turgay

2012/6/8 Mark Abraham 

>  On 8/06/2012 10:40 PM, Turgay Cakmak wrote:
>
>
> Hi all,
>
> I have encountered a problem using eneconv with -dt flag.
> I have several energy files. Firstly, I concatenated them using the
> following.
>
> eneconv  -f  first.edr   second.edr    -settime  -o  total.edr
>
> But, total.edr is so large. To reduce the size of file, I used the below
> command:
>
> eneconv -f  total.edr  -dt  100 -o  total_dt100.edr
>
> But, I get the following warning:
>* .*
>* .*
>
> WARNING: missing energy sums at time 69984.00
>
> WARNING: missing energy sums at time 69988.00
>
> WARNING: missing energy sums at time 69992.00
>
> WARNING: missing energy sums at time 69996.00
> Reading energy frame  35000 time 7.000
> WARNING: missing energy sums at time 7.00
> Last energy frame read 35000 time 7.000
> Last step written from energy.edr: t 7, step 3500
>
> Last frame written was at step at step 3500, time 7.00
> Wrote 701 frames
>
>
> At some point something did something to break the usual format of the
> files... but we can have no idea what.
>
>
>
> Eventually, I get the "total_dt100.edr". But, I am not sure I could
> continue using this file.
> If someone could help me, I would be glad.
>
>
> Depends what you want to use it for.
>


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

[gmx-users] Furthest distance between 2 atoms

[Catarina Santos]

Hi everybody,

I am trying to find an easy way to calculate the furthest distance between
the atoms of two groups in a system. I know that g_dist can give me all the
distances I want, but I would have to run g_dist for each possible atom
pair (which is not a fast or a straightforward process).

I have already checked the manual and tried some tools, I did not
find anything that returns all the distances between all the possible pairs
of atoms though. Am I missing something?

Thanks in advance.

Regards,

Catarina Santos
Research Student
Molecular Simulation Group
Instituto de Tecnologia Química e Biológica
Universidade Nova de Lisboa
Av.da República, Apartado 127
2781-901 Oeiras, Portugal
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] qmmm and "real-life" meaning of null

[Gerrit Groenhof]

Hi,

INdeed, you would still need a LJ term for the QM atoms, in order to interact 
with the MM atoms. Only if the radius around your special atom were larger than 
the LJ cut-off could you safely ignore the LJ paramter.

gerrit

> 
> Message: 3
> Date: Thu, 07 Jun 2012 19:47:02 -0400
> From: "Justin A. Lemkul" 
> Subject: Re: [gmx-users] qmmm and "real-life" meaning of null
>   lennard-jones   parameters
> To: Discussion list for GROMACS users 
> Message-ID: <4fd13d76.9040...@vt.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
> 
> 
> On 6/7/12 7:28 PM, Edward Deira wrote:
>> Dear all,
>> 
>> I'm currently starting to dwell deeper in MD, and I'm taking some time to
>> understand what's going on inside the gromacs "black-box".
>> In one of those dwellings, I came across an older post
>> [http://www.mail-archive.com/gmx-users@gromacs.org/msg42568.html] which 
>> reads:
>> 
>> Question:
>> 4. In ffnonbonded.itp, why are both sigma and epsilon set to zero for HW
>> (opls_117)? This seems to imply that, as far as Lennard-Jones interactions 
>> are
>> concerned, the hydrogens on the waters don't exist. Or, in other words, in 
>> the
>> absence of charges, the hydrogens don't "feel" the hydrogens, the hydrogens
>> don't "feel" the oxygens, and the oxygens don't "feel" the hydrogens. In 
>> other
>> words, the hydrogens interact with the world only via electrostatic 
>> (Coulombic)
>> interactions. Is this a correct interpretation?Correct. Many force fields do 
>> this.
>> Answer:
>> 
>> So, my question, if a question at all:
>> 
>> Suppose I have a regular protein and put inside some metal atom that will
>> coordinate with some O and N atoms from the side chains. If the sigma and
>> epsilon for that metal are null, than the metal - sidechains interactions are
>> exclusively electrostatic. Does this make sense ? What are the implications 
>> of
>> this for the "coordination chemistry" of that "metal - sidechain complex" ?
>> 
>> On the side: suppose I want some non parameterized metal atom, say W, for 
>> which
>> I will compute all the other parameters in the same/similar way described in 
>> the
>> force field papers, but for which no experimental data are available for me 
>> to
>> compare computable meaningful sigma and epsilon values. Can I just sigma and
>> epsilon to zero ? Or should I do qmmm to have W in the qm part ?
>> 
> 
> The fact that the LJ parameters for H are zero derives from its size.  The 
> environment is more strongly influenced by the heavy atom to which H is 
> bonded. 
>  In the case of a larger metal ion, I would seriously doubt that setting LJ 
> parameters to zero is valid.  It's quite convenient, but in most force 
> fields, 
> all metal ions have some LJ parameters.  Perhaps investigating how those 
> parameters were derived would be useful.  For what it's worth, I believe the 
> origin of the zero-LJ H parameters comes from this work:
> 
> http://pubs.acs.org/doi/abs/10.1021/ja00824a004
> 
>> Also, from the few tutorials and from the manual, I have the impression that
>> even for qmmm with gromacs and mopac i still need force field parameters for 
>> the
>> qm part, is this true ? Or i just need to include a qmmm section in all mdp
>> files, including the first ion adding and energy minimization steps ? Sorry 
>> for
>> the naivety in this, but i've only made "regular protein" MD so far.
>> 
> 
> I've never done any QM/MM, but my assumption would be that you have to have 
> some 
> valid topology to start with.  Perhaps someone else can comment on this 
> methodological issue.
> 
> -Justin
> 
> -- 
> 
> 
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> 

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

[gmx-users] 4. QM/MM Calculation with Orca (Minos Matsoukas)

[Gerrit Groenhof]

YOucan only use one thread in mdrun, but more than one in ORCA. Try to estimate 
the ratio of computation time spent between the QM and MM calculation to get an 
idea of why we never bothered to parallellize the MM part.

Hope this helps,

Gerrit


> 
>   4. QM/MM Calculation with Orca (Minos Matsoukas)
>   
> 
> Dear GROMACS Users,
> 
> I am using Gromacs with ORCA for a qm calculation, but I can only use
> 1 processor for mdrun.
> 
> If I place the line !PAL4 in the topol.ORCAINFO and run mdrun -nt 1,
> it works, and ORCA parallelizes correctly for 4 processors, although
> mdrun only runs one thread.
> 
> If I place the line !PAL4 in the topol.ORCAINFO and run mdrun -nt 4,
> it doesn‚t work. mdrun crashes with a segmentation fault error :
> 
> Back Off! I just backed up md.log to ./#md.log.35#
> Reading file topol.tpr, VERSION 4.5.5 (single precision)
> Starting 4 threads
> QM/MM calculation requested.
> QM/MM calculation requested.
> QM/MM calculation requested.
> QM/MM calculation requested.
> there we go!
> there we go!
> there we go!
> there we go!
> Layer 0
> nr of QM atoms 73
> QMlevel: B3LYP/STO-3G
> 
> /opt/soft/orca_2_9_1_linux_x86-64/opt/soft/orca_2_9_1_linux_x86-64...
> orca initialised...
> Segmentation fault (core dumped)
> 
> 
> Is this because mdrun can only run with 1 thread when used with ORCA ?
> 
> Here‚s some information :
> OS: CentOS 6.2 x86_64 kernel 2.6.32-220
> GROMACS: 4.5.5 compiled with options : --with-qmmm-orca 
> --without-qmmm-gaussian
> ORCA: 2.9.1 for x86_64
> OPENMPI: 1.4.3
> 
> My minimization mdp file is:
> 
> ;
> ; Input file
> ;
> title   =  Yo ; a string
> cpp =  /lib/cpp   ; c-preprocessor
> integrator  =  cg ;
> ;integrator   = l-bfgs
> 
> rlist   =  1.0; cut-off for ns
> rvdw=  1.0; cut-off for vdw
> rcoulomb=  1.5; cut-off for coulomb
> ;   Energy minimizing stuff
> ;
> nsteps  =  200
> emtol   =  10
> emstep  =  0.1
> ;define  = -DPOSRES
> constraints = none
> QMMM = yes
> QMMM-grps= MOL
> QMMMscheme   = normal
> QMbasis  = STO-3g
> QMmethod = b3lyp
> QMcharge = 1
> QMmult   = 1
> -
> 
> 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

Re: [gmx-users] Re: energy conservation: shift vs shifted user potential

[Mark Abraham]

On 8/06/2012 5:59 PM, Anja Kuhnhold wrote:

Hello all,

can someone give me a hint, please?
Do you need more information?
Has anyone had a similar problem.

I really need to figure that out, because I will simulate other systems
which can be run only with user tables.

Anja



-

- Original Message -
- Original Message -
From Anja Kuhnhold
Date Wed, 06 Jun 2012 15:07:40 +0200
To gmx-users@gromacs.org
Subject [gmx-users] energy conservation: shift vs shifted user potential
Dear gmx-users,

I have a problem concerning energy conservation when using user
potentials (tables).
I'm using gromacs 4.5.4
I simulate a simple Lennard-Jones(6-12) +Fene polymer melt (1600
chains a 10 beads in a 26x26x26 periodic box).

I tried different vdwtypes (cutoff always 3.24):
The cut-off version does not conserve energy -- okay.
The shift and switch versions conserve energy -- fine.

Now I wanted to do the same with user tables:
Simple Lennard-Jones table gives really the same results as the
cut-off version -- good.

But if I use a table with shifted Lennard-Jones potential it is not
comparable to the shift version
and the energy is not conserved -- ?

I use a shift function as written in the manual (chapter 4.1.5) --
there must be a factor alpha added in the constants A and B --
(r1=0).

The parameters are the same for shift version and shifted user version.

Has someone an idea why the shifted user potential doesn't work in
this way?


We've no real idea what you've done... manual 6.7.2 describes the 
required format and there are (unshifted) examples in your installation 
under $GMXLIB/share/top/table*.xvg.


Makr



Here is the mdp:


integrator  = md-vv
dt  = 0.0035
nsteps  = 1000
nstxout = 1
nstvout = 1
nstfout = 1
nstlog  = 1
ns_type = grid
pbc = xyz
rvdw= 3.24
rlist   = 3.6
tcoupl  = no
tc-grps = System
tau_t   = 2.0
ref_t   = 127.2717
vdwtype = user;Shift
rcoulomb= 3.6;2.24;1.12
coulombtype = Cut-off
rvdw-switch = 0.0

energygrps  = bead
energygrp_table = bead bead


Thanks in advance
Anja



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

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

Re: [gmx-users] eneconv with -dt flag

[Mark Abraham]

On 8/06/2012 10:40 PM, Turgay Cakmak wrote:


Hi all,

I have encountered a problem using eneconv with -dt flag.
I have several energy files. Firstly, I concatenated them using the 
following.


eneconv-f first.edr   second.edr    -settime -o total.edr

But, total.edr is so large. To reduce the size of file, I used the 
below command:


eneconv -f  total.edr  -dt  100 -o  total_dt100.edr

But, I get the following warning:
*.*
*.*

WARNING: missing energy sums at time 69984.00

WARNING: missing energy sums at time 69988.00

WARNING: missing energy sums at time 69992.00

WARNING: missing energy sums at time 69996.00
Reading energy frame  35000 time 7.000
WARNING: missing energy sums at time 7.00
Last energy frame read 35000 time 7.000
Last step written from energy.edr: t 7, step 3500

Last frame written was at step at step 3500, time 7.00
Wrote 701 frames


At some point something did something to break the usual format of the 
files... but we can have no idea what.




Eventually, I get the "total_dt100.edr". But, I am not sure I could 
continue using this file.

If someone could help me, I would be glad.


Depends what you want to use it for.

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

Re: [gmx-users] md.log i/o error in futil.c +459

[Mark Abraham]

On 8/06/2012 11:05 PM, Inon Sharony wrote:


Good afternoon!

I'm getting the following error message from mdrun, and despite 
trying, I can't seem to find the file futil.c in which the error is 
supposed to be taking place. Can anyone help me out, please?




It's the source file of GROMACS that produce the code from which the 
error originated.



Thanks, and have a good weekend!

Inon.


 :-)  G  R  O  M  A  C  S  (-:


  GROningen Mixture of Alchemy and Childrens' Stories

:-)  VERSION 4.5.4  (-:

Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
  Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra,
Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff,
   Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz,
Michael Shirts, Alfons Sijbers, Peter Tieleman,

   Berk Hess, David van der Spoel, and Erik Lindahl.

   Copyright (c) 1991-2000, University of Groningen, The Netherlands.
Copyright (c) 2001-2010, The GROMACS development team at
Uppsala University & The Royal Institute of Technology, Sweden.
check out http://www.gromacs.org for more information.

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

  :-)  mdrun_d (double precision)  (-:


---
Program mdrun_d, VERSION 4.5.4
Source code file: futil.c, line: 459

File input/output error:
md.log
For more information and tips for troubleshooting, please check the 
GROMACS

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



mdrun can't do whatever it should be doing with this file. Probably you 
lack permissions or file space.


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] md.log i/o error in futil.c +459

[Inon Sharony]

  
  
Good afternoon!
I'm getting the following error message from mdrun, and despite
  trying, I can't seem to find the file futil.c in which the error
  is supposed to be taking place. Can anyone help me out, please?

Thanks, and have a good weekend!
Inon.


 :-)  G  R  O  M  A  C  S  (-:

    GROningen Mixture of Alchemy and Childrens' Stories
  
      :-)  VERSION 4.5.4  (-:
  
      Written by Emile Apol, Rossen Apostolov, Herman J.C.
  Berendsen,
    Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton
  Feenstra, 
      Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter
  Meulenhoff, 
     Teemu Murtola, Szilard Pall, Sander Pronk, Roland
  Schulz, 
      Michael Shirts, Alfons Sijbers, Peter Tieleman,
  
     Berk Hess, David van der Spoel, and Erik Lindahl.
  
     Copyright (c) 1991-2000, University of Groningen, The
  Netherlands.
      Copyright (c) 2001-2010, The GROMACS development team
  at
      Uppsala University & The Royal Institute of
  Technology, Sweden.
      check out http://www.gromacs.org for more information.
  
   This program is free software; you can redistribute it
  and/or
    modify it under the terms of the GNU General Public
  License
   as published by the Free Software Foundation; either
  version 2
   of the License, or (at your option) any later
  version.
  
    :-)  mdrun_d (double precision)  (-:
  
  
  ---
  Program mdrun_d, VERSION 4.5.4
  Source code file: futil.c, line: 459
  
  File input/output error:
  md.log
  For more information and tips for troubleshooting, please check
  the GROMACS
  website at http://www.gromacs.org/Documentation/Errors
  ---
  
  "Everything Must Go" (Red Hot Chili Peppers)
  

-- 
Inon Sharony
ינון   שרוני
972-3-6407634
Please consider your environmental responsibility before printing this e-mail.
  

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

[gmx-users] QM/MM Calculation with Orca

[Minos Matsoukas]

Dear GROMACS Users,

I am using Gromacs with ORCA for a qm calculation, but I can only use
1 processor for mdrun.

If I place the line !PAL4 in the topol.ORCAINFO and run mdrun -nt 1,
it works, and ORCA parallelizes correctly for 4 processors, although
mdrun only runs one thread.

If I place the line !PAL4 in the topol.ORCAINFO and run mdrun -nt 4,
it doesn’t work. mdrun crashes with a segmentation fault error :

Back Off! I just backed up md.log to ./#md.log.35#
Reading file topol.tpr, VERSION 4.5.5 (single precision)
Starting 4 threads
QM/MM calculation requested.
QM/MM calculation requested.
QM/MM calculation requested.
QM/MM calculation requested.
there we go!
there we go!
there we go!
there we go!
Layer 0
nr of QM atoms 73
QMlevel: B3LYP/STO-3G

/opt/soft/orca_2_9_1_linux_x86-64/opt/soft/orca_2_9_1_linux_x86-64...
orca initialised...
Segmentation fault (core dumped)


Is this because mdrun can only run with 1 thread when used with ORCA ?

Here’s some information :
OS: CentOS 6.2 x86_64 kernel 2.6.32-220
GROMACS: 4.5.5 compiled with options : --with-qmmm-orca --without-qmmm-gaussian
ORCA: 2.9.1 for x86_64
OPENMPI: 1.4.3

My minimization mdp file is:

;
;   Input file
;
title   =  Yo   ; a string
cpp =  /lib/cpp ; c-preprocessor
integrator  =  cg   ;
;integrator   = l-bfgs

rlist   =  1.0  ; cut-off for ns
rvdw=  1.0  ; cut-off for vdw
rcoulomb=  1.5  ; cut-off for coulomb
;   Energy minimizing stuff
;
nsteps  =  200
emtol   =  10
emstep  =  0.1
;define  = -DPOSRES
constraints = none
QMMM = yes
QMMM-grps= MOL
QMMMscheme   = normal
QMbasis  = STO-3g
QMmethod = b3lyp
QMcharge = 1
QMmult   = 1
-

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] eneconv with -dt flag

[Turgay Cakmak]

Hi all,

I have encountered a problem using eneconv with -dt flag.
I have several energy files. Firstly, I concatenated them using the
following.

eneconv  -f  first.edr   second.edr    -settime  -o  total.edr

But, total.edr is so large. To reduce the size of file, I used the below
command:

eneconv -f  total.edr  -dt  100 -o  total_dt100.edr

But, I get the following warning:
   * .*
   * .*

WARNING: missing energy sums at time 69984.00

WARNING: missing energy sums at time 69988.00

WARNING: missing energy sums at time 69992.00

WARNING: missing energy sums at time 69996.00
Reading energy frame  35000 time 7.000
WARNING: missing energy sums at time 7.00
Last energy frame read 35000 time 7.000
Last step written from energy.edr: t 7, step 3500

Last frame written was at step at step 3500, time 7.00
Wrote 701 frames

Eventually, I get the "total_dt100.edr". But, I am not sure I could
continue using this file.
If someone could help me, I would be glad.

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

Re: [gmx-users] Pulling - ligand outside the box

[Mark Abraham]

On 8/06/2012 8:35 PM, Steven Neumann wrote:

Dear Gmx Users,

I pulled my ligand away from the protein and I found out that after 
getting rid of PBC (trjconv -pbc whole)


There is no magic way to "get rid of PBC", and "-pbc whole" is certainly 
not it. Read trjconv -h and see 
http://www.gromacs.org/Documentation/Terminology/Periodic_Boundary_Conditions 
for suggested workflows.



my ligand goes outside the box?


There is no outside. Since your ligand and your protein were probably 
not part of the same [moleculetype], "-pbc whole" can't magically 
understand you think they should be represented relative to the same 
periodic cell. Correct usage of "-pbc cluster" might help.



Is it fine or should I increase my box as I want to run umbrella sampling?


We can't possibly know.

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] Pulling - ligand outside the box

[Steven Neumann]

Dear Gmx Users,

I pulled my ligand away from the protein and I found out that after getting
rid of PBC (trjconv -pbc whole) my ligand goes outside the box? Is it fine
or should I increase my box as I want to run umbrella sampling?

Thanks,

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

[gmx-users] Re: energy conservation: shift vs shifted user potential

[Anja Kuhnhold]

Hello all,

can someone give me a hint, please?
Do you need more information?
Has anyone had a similar problem.

I really need to figure that out, because I will simulate other systems
which can be run only with user tables.

Anja



-
> - Original Message -
> - Original Message -
> From Anja Kuhnhold 
> Date Wed, 06 Jun 2012 15:07:40 +0200
> To gmx-users@gromacs.org
> Subject [gmx-users] energy conservation: shift vs shifted user potential
> Dear gmx-users,
> 
> I have a problem concerning energy conservation when using user 
> potentials (tables).
> I'm using gromacs 4.5.4
> I simulate a simple Lennard-Jones(6-12) +Fene polymer melt (1600 
> chains a 10 beads in a 26x26x26 periodic box).
> 
> I tried different vdwtypes (cutoff always 3.24):
> The cut-off version does not conserve energy -- okay.
> The shift and switch versions conserve energy -- fine.
> 
> Now I wanted to do the same with user tables:
> Simple Lennard-Jones table gives really the same results as the 
> cut-off version -- good.
> 
> But if I use a table with shifted Lennard-Jones potential it is not 
> comparable to the shift version
> and the energy is not conserved -- ?
> 
> I use a shift function as written in the manual (chapter 4.1.5) -- 
> there must be a factor alpha added in the constants A and B --
> (r1=0).
> 
> The parameters are the same for shift version and shifted user version.
> 
> Has someone an idea why the shifted user potential doesn't work in 
> this way?
> 
> Here is the mdp:
> 
> 
> integrator= md-vv
> dt= 0.0035
> nsteps= 1000
> nstxout   = 1
> nstvout   = 1
> nstfout   = 1
> nstlog= 1
> ns_type   = grid
> pbc   = xyz
> rvdw  = 3.24
> rlist = 3.6
> tcoupl= no
> tc-grps   = System
> tau_t = 2.0
> ref_t = 127.2717
> vdwtype   = user;Shift
> rcoulomb  = 3.6;2.24;1.12
> coulombtype   = Cut-off
> rvdw-switch   = 0.0
> 
> energygrps= bead
> energygrp_table   = bead bead
> 
> 
> Thanks in advance
> Anja
> 

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