[gmx-users] Reordering of water molecules with Na+ ions coordinates in xtc file

[chris . neale]

Dear Leila:

Perhaps I took the long way around, but I am not aware of any such  
tool. You can make custom modifications to things like trjconv as  
follows.


!!! Please note: I didn't test this except to see that it compiles.  
You should run some analyses before and after switching the order to  
check that the file hasn't been corrupted.


This is done with gromacs 4.0.7.

cd src/tools
cp gmx_trjconv.c ../../exec/share/gromacs/template/my_tool.c

### Add the following text to the bottom of the my_tool.c file:

int
main(int argc,char *argv[])
{
  return gmx_trjconv(argc,argv);
}


### Now go to line 1246 of the file and replace this:

  for(i=0; i### That's not going to modify the forces and velocities, but I  
assume that this will be ok for you. If you need them, you can add  
that part by analogy.


### Now create the makefile:

sed "s/template/my_tool/g" Makefile.x86_64-unknown-linux-gnu >  
Makefile.my_tool


### now compile it:

make -f Makefile.my_tool

-- original message --

Dear gromacs users

My simulation system contains protein, dna, water molecules and Na+
ions respectively.

I want to reorder water molecules with Na+ ions in final xtc file as

at first

1-1867   complex (protein and dna)
1868 - 24085  SOL (water molecules)
24086 - 24099Na+ (ions)

after reordering

1-1867   complex (protein and dna)
1868 -1881 Na+ (ions)
1882 - 24099  SOL (water molecules)

How to do that?

I need to this reordering to obtain amber trajectory file (mdcrd) by VMD.

any help and suggestion will highly appreciated.


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

[gmx-users] Add gromacs forcefield w/ virtual site

[chris . neale]

Sorry Marcelo, that sounds like a job for an author ;) If nobody else  
chimes in, then I suggest that you do some testing and come back to  
the list with specific problems. As an analogy, you might look at how  
the virtual site is handled in tip4p in the absence of settle.


Chris.

-- original message --

Thank you Chris for your answer:

1) The molecule has no net charge because the virtual site in the center
of mass is a point charge twice the charge in the O atom.
2) Until now I've created 5 files but I don't know if I am doing the
right thing:

*   forcefield.itp

#define _FF_OXY

[ defaults ]
; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ
12no
1.01.0

#include "ffnonbonded.itp"
#include "ffbonded.itp"

*ffbonded.itp

[ bondtypes ]
; ij   func
   OO 5  (I've used function 5 because the forcefield has a
fixed length, is it right?) ;

*ffnonbonded.itp

[ atomtypes ]
; name  at.nummass  charge   ptype
sigma  epsilon
  O8   15.99940-0.123   A
3.01300e-014.0780822e-01
COM0   0.00 0.246   V
0.0e+000.000e+00

*   atomtypes.atp

O 15.99940  ; Oxygen Atom
COM   0.0   ; Virtual site (COM charge)

*   residues.rtp

; Oxygen
[ OXY ]
  [ atoms ]
  O O  -0.1231
  COM   COM0.2461
  [ bonds ]
 OO

3) Is it necessary more or less files? And in which file do I
specificate de virtual site?
4) The virtual site in the COM is [ virtual_sites2 ] with a = 0.5 or [
virtual_sitesn ]?



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

[gmx-users] Reordering of water molecules with Na+ ions coordinates in xtc message

[chris . neale]

Dear Leila:

I looked at the source file that you sent me by personal email looks correct.

Did you use the correct version of the template directory? My makefile  
looks like this (see below). If you used I am not sure where you gt  
the makefile since you mentioned many different places in your online  
posts and by email, but be sure that you are using the one produced by  
the 4.0.7 compilation.


Chris.

# Generated automatically from Makefile.in by configure.
#
# This is a Gromacs 3.0 my_tool makefile for your own utility programs.
#
# Copy this file to whatever directory you are using for your own
# software and add more targets like the my_tool one below.
#
# If you are using gmake it is relatively straightforward to add
# an include based on environment variables (like previous Gromacs versions)
# to select compiler flags and stuff automatically, but below it is static:
#

# Variables set by the configuration script:
LIBS = -lmd -lgmx -lnsl -lfftw3f -lm
LDFLAGS  =  
-L/project/pomes/cneale/GPC/exe/intel/gromacs-4.0.7/exec/lib  
-L/project/pomes/cneale/GPC/exe/intel/fftw-3.1.2/exec/lib
CFLAGS   = -O3 -fomit-frame-pointer -finline-functions -Wall  
-Wno-unused -funroll-all-loops   
-I/project/pomes/cneale/GPC/exe/intel/gromacs-4.0.7/exec/include  
-I/project/pomes/cneale/GPC/exe/intel/gromacs-4.0.7/exec/include/gromacs

CC   = cc
LD   = $(CC)

# The real make targets - note that most make programs support
# the shortcut $^ instead of listing all object files a second
# time, but we cannot count on it...

my_tool:my_tool.o
$(LD) $(LDFLAGS) -o $@ my_tool.o $(LIBS)


I found gmx_trjconv.c file, It is in gromacs-4.0.7/src/tools.
since location of template directory is as follows:
gromacs-4.0.7/share/template,
I used cp gmx_trjconv.c ../../share/template/my_tool.c.
is it true?
I did those changes you said in my_tool.c file.
I created makefile without problem using sed "s/template/my_tool/g"
Makefile.x86_64-unknown-linux-gnu > Makefile.my_tool.
but when compile it, I encountered with following errors:

cc -O3 -fomit-frame-pointer -finline-functions -Wall -Wno-unused
-funroll-all-loops   -I/usr/include/libxml2 -I/usr/local/gromacs/include
-I/usr/local/gromacs/include/gromacs -c -o my_tool.o my_tool.c
my_tool.c: In function ?calc_pbc_cluster?:
my_tool.c:88: error: ?bool? undeclared (first use in this function)

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

[gmx-users] Umbrella sampling

[chris . neale]

Nilesh,

In the time it take Justin to answer, you can already have tested some  
of this out for yourself. Pick a single window and run the pull-code  
dynamics and then look at the trajectory with VMD and plot the  
coord.xvg file and check the distance by g_dist. You'll learn a lot  
more like this.


Chris.

-- original message --

 Justin,
I considered single ion- pair.
sorry for typing mistake (I wrote 4 ion-pairs).

How can I define distance range.
ONe more my cation and anion are close in box edges in my equilibrated
strucutre.

By defining pull_rate1  = -0.01 can cation and anion will come close
to each other?

Nilesh


Does will it affect on results?

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

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

[gmx-users] Reordering of water molecules with Na+ ions coordinates in xtc message

[chris . neale]

Dear Leila,

Please address my previous question. *Are you using the correct  
Makefile.template?* Did you do the compilation of gromacs? If not, I  
suggest that you compile gromacs again (and not into /usr/ but into  
./exec or some such local dir) and use your files from there.


If you still can't get the file to work in this way, then just  
download a fresh version of gromacs4.0.7, name it with a special  
directory, make the modification directly to the *original* trjconv,  
and compile it. Now you will have a special version of trjconv.


I can't help you any more because the steps that I suggested worked  
for me. It is unclear exactly what you are doing wrong, but I suspect  
that it is simply that you are not actually following the steps that I  
laid out.


Chris.

-- original message --

I did all of steps again. this tims when compile Makefile.my_tool by both of
make and gmake:

make: *** No rule to make target `my_tool.o', needed by `my_tool'.  Stop.

my Makefile.my_tool is as follows:

# Generated automatically from Makefile.in by configure.
#
# This is a Gromacs 3.0 my_tool makefile for your own utility programs.
#
# Copy this file to whatever directory you are using for your own
# software and add more targets like the my_tool one below.
#
# If you are using gmake it is relatively straightforward to add
# an include based on environment variables (like previous Gromacs versions)
# to select compiler flags and stuff automatically, but below it is static:
#

# Variables set by the configuration script:
LIBS = -lmd -lgmx -lxml2  -L/usr/lib64 -lnsl -lfftw3f -lm   -lSM
-lICE -lX11
LDFLAGS  = -L/usr/local/gromacs/lib
CFLAGS = -O3 -fomit-frame-pointer -finline-functions -Wall
-Wno-unused -funroll-all-loops   -I/usr/include/libxml2
-I/usr/local/gromacs/include -I/usr/local/gromacs/include/gromacs
CC   = cc
LD   = $(CC)

# The real make targets - note that most make programs support
# the shortcut $^ instead of listing all object files a second
# time, but we cannot count on it...

my_tool:my_tool.o
$(LD) $(LDFLAGS) -o $@ my_tool.o $(LIBS)
-- next part --
An HTML attachment was scrubbed...

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

[gmx-users] eordering of water molecules with Na+ ions coordinates in xtc file

[chris . neale]

It seems like Tsjerk has a better solution. It appears that you need  
only to have an index group in which the atoms are ordered in the  
output order that you desire. Therefore, if you wanted to reverse the  
order of the first 3 atoms, you would create an index file like this:


[ my_group ]
3 2 1

Just create something like that for your group. Looks like Tsjerk has  
even showed you how to create the index file that you want.


You would use this with the regular version of trjconv.

I didn't know that this was possible.
Chris.

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

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

[gmx-users] creating mdrun_mpi with cmake

[chris . neale]

Dear users:

I have successfully compiled gromacs 4.5.3 with the description posted  
here:  
http://www.gromacs.org/Developer_Zone/Cmake#An_example_installation


How do I compile the mpi version of mdrun? I tried simply adding  
openmpi libraries to my paths, but this was not enough. Is there a  
special cmake flag? I could not find on one the internet.


Second, Do I even need to compile with mpi in gromacs 4.5.3? It seems  
to me that non-mpi mdrun can do this with threads (at least over cores  
on a single box).


Thank you,
Chris.

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

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

[gmx-users] problem when calculating the solvation free energy of ARG in water using PME.

[chris . neale]

1. How sure are you that you have reproduced the run exactly on the  
different machines?

2. Are there any other differences besides machine? Number of cores?
3. What's your water model?
4. What else is in your system? (counterions?)
5. What are your crash/not-crash statistics? The OS difference might  
just be noise.


I realize that this could be some machine-dependent error (as you seem  
to be suggesting) but I would not be able to help you with that. Let's  
ensure that it is not some other difference and then you could file a  
redmine issue if the problem can not be resolved. You likely have a  
good small test case.


Chris.

-- original message --


I am trying to calculate the solvation free energy of ARG in water, using
GROMACS 4.0.7 (using TI).


I found that when I run the calculations on linux machines (I try few) it
run OK, but when I run the calculations on OSX it crashed in most lambda
points. I run the calculations using PME and RF, and the problem is
consisted only with PME. I also done the same thing for ALA and it works
fine. So I guess it has a connection to the charges.


 Any advices?


Best regards,



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

[gmx-users] Melting Temperatue

[chris . neale]

This one comes to mind:

D. Paschek and A.E. Garcia, ?Reversible temperature and pressure  
denaturation of a protein fragment: Replica exchange molecular  
dynamics simulation study.? Phys. Rev. Let. 93: 238105 (2004)


-- original message --

Hello,

Experimentals often use heating to define strength of contacts.

Has anyone tried to calculate this using MD/GROMACS or can anyone refer me
to a work which was done on this subject?

Thanks,
Ifat

M.Sc. Student
Bar Ilan University


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

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

[gmx-users] creating mdrun_mpi with cmake

[chris . neale]

Thanks for your help Justin. I have updated the instructions at:

http://www.gromacs.org/Developer_Zone/Cmake#An_example_installation

with what has worked for me.


chris.neale at utoronto.ca wrote:

Dear users:

I have successfully compiled gromacs 4.5.3 with the description  
posted here:  
http://www.gromacs.org/Developer_Zone/Cmake#An_example_installation


How do I compile the mpi version of mdrun? I tried simply adding  
openmpi libraries to my paths, but this was not enough. Is there a  
special cmake flag? I could not find on one the internet.




-DGMX_MPI=ON

Other relevant options may or may not include:

-DGMX_DEFAULT_SUFFIX
-DMPI_COMPILER
-DMPI_INCLUDE_PATH

I don't know what you've already found, but I've had to use all of  
these to get

a proper mdrun with a custom suffix (with -DGMX_DEFAULT_SUFFIX=OFF).

Second, Do I even need to compile with mpi in gromacs 4.5.3? It  
seems to me that non-mpi mdrun can do this with threads (at least  
over cores on a single box).




That's my understanding, as well.

-Justin


Thank you,
Chris.



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

[gmx-users] number of cores with a large prime factor

[chris . neale]

Dear users:

I received this error message that gromacs 4.5.3 does not allow me to  
run on 34 cores. I suggest that this gets put in as a warning (so that  
I can avoid it if I choose and actually run on 34 cores).


A run on 32 cores -npme -1 choose to use 24 particle-particle nodes  
and 8 PME nodes. I therefore selected 34 cores and -nmpe 10. I don't  
see how the large prime factor is relevant when -nmpme is not equal to  
zero.


The error message was:

Fatal error:
The number of nodes you selected (34) contains a large prime factor  
17. In most cases this will lead to bad performance. Choose a number  
with smaller prime factors or set the decomposition (option -dd)  
manually.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

Thoughts?


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

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

[gmx-users] Slow Runs

[Chris Neale]

In addition, you're only updating your neighbourlist every 40 ps. If 
you're going to use a 4 fs timestep, I suggest that you use nstlist=5. 
Also, you appear to not be using any constraints while you are using a 4 
fs timestep.


I suggest that you stop worrying about why the run is slower and do some 
tutorials and read the manual and ensure that you understand the 
consequences of your .mdp options.


Chris.

-- original message --

Denny Frost wrote:

/  about 12000 atoms, 8 nodes, CentOS 5.3/Linux, infiniband.  Below is a

/>/  copy of my mdp file.
/>/
/>/  title   =  BMIM+PF6
/>/  cpp =  /lib/cpp
/>/  constraints =  all_bonds
/>/  integrator  =  md
/>/  dt  =  0.004   ; ps !
/>/  nsteps  =  2000   ; total 4ns.
/>/  nstcomm =  1
/>/  nstxout =  5
/>/  nstvout =  5
/>/  nstfout =  0
/>/  nstlog  =  5000
/>/  nstenergy   =  5000
/>/  nstxtcout   =  25000
/>/  nstlist =  10
/>/  ns_type =  grid
/>/  pbc =  xyz
/>/  coulombtype =  PME
/>/  vdwtype =  Shift
/>/  rlist   =  1.0
/>/  rcoulomb=  1.0
/>/  rvdw=  1.0
/>/  fourierspacing  =  0.6
/
This fourierspacing is 5-6 times larger than what is normally accepted as
sufficiently accurate.  A sparse grid will make the PME algorithm faster,
actually, but at the expense of accuracy.

Can you post the domain decomposition statistics from the .log file?  They
appear just above the energies from time 0.  What did grompp tell you about the
relative PME:PP load?

-Justin


/  ;pme_order   =  4

/>/  ewald_rtol  =  1e-5
/>/  ; Berendsen temperature coupling is on in two groups
/>/  Tcoupl  =  berendsen
/>/  tc_grps =  BMI  PF6
/>/  tau_t   =  0.1  0.1
/>/  ref_t   =  300  300
/>/  nsttcouple  =  1
/>/  ; Energy monitoring
/>/  energygrps  =  BMI  PF6
/>/  ; Isotropic pressure coupling is now on
/>/  Pcoupl  =  berendsen
/>/  pcoupltype  =  isotropic
/>/  ;pc-grps =  BMI  PFF
/>/  tau_p   =  1.0
/>/  ref_p   =  1.0
/>/  compressibility =  4.5e-5
/>/
/>/  ; Generate velocites is off at 300 K.
/>/  gen_vel =  yes
/>/  gen_temp=  300.0
/>/  gen_seed=  10
/>/
/>/
/>/  On Thu, Jan 27, 2011 at 4:12 PM, Dallas Warrenhttp://lists.gromacs.org/mailman/listinfo/gmx-users>
/>/  >>  wrote:
/>/
/>/  You will need to provide more details on the system.  How many
/>/  atoms, what sort of computer system is it being run on, how many
/>/  nodes, copy of the mdp file etc.
/>/
/>/
/>/
/>/  Catch ya,
/>/
/>/  Dr. Dallas Warren
/>/
/>/  Medicinal Chemistry and Drug Action
/>/
/>/  Monash Institute of Pharmaceutical Sciences, Monash University
/>/  381 Royal Parade, Parkville VIC 3010
/>/  dallas.warren at monash.edu   
 >
/>/
/>/  +61 3 9903 9304
/>/  -
/>/  When the only tool you own is a hammer, every problem begins to
/>/  resemble a nail.
/>/
/>/
/>/
/>/  *From:*gmx-users-bounces at gromacs.org  

/>/  >
/>/  [mailto:gmx-users-bounces at gromacs.org  

/>/  >] *On Behalf Of *Denny Frost
/>/  *Sent:* Friday, 28 January 2011 9:34 AM
/>/  *To:* Discussion list for GROMACS users
/>/  *Subject:* [gmx-users] Slow Runs
/>/
/>/
/>/
/>/  I am taking over a project for a graduate student who did MD using
/>/  Gromacs 3.3.3.  I now run similar simulations with Gromacs 4.5.1 and
/>/  find that they run only about 1/2 to 1/3 as fast as the previous
/>/  runs done in Gromacs 3.3.3.  The runs have about the same number of
/>/  atoms and both use opls force fields.  The mdp files is virtually
/>/  the same (I copied them).  The only major difference is that my runs
/>/  have difference species and thus have different (although smaller)
/>/  itp files.  The runs are stable and give reasonable thermodynamic
/>/  properties - they're just slow.  Has anyone had any experience with
/>/  something like this?
/>/
/>/
/>/  --
/>/  gmx-users mailing listgmx-users at gromacs.org  

/>/  

[gmx-users] water segment in the z-direction

[chris . neale]

Dear Tusar:

Let's keep this on the mailing list.

You could learn about expr with a simple google search! It should have  
worked though, did you mangle the command somehow when you typed it?  
There should be " > not_last_line.gro" with a redirect at the end but  
you have only a star. Nevertheless, my scripting has improved a bit  
since I wrote that and you could replace 7 by:


head -n -1 initial.gro > not_last_line.gro

And, to answer your second question, you use keep_these_waters.gro in  
part 8, along with other files, to create new_system.gro


Chris.

-- original message --

I am trying to increase the size of the water segment in the z-direction.
For this purpose I am following your recommended steps and the script.

1. run genbox on initial.gro to create solvated.gro
2. cp solvated.gro new_waters.gro
3. use vi to remove everything in new_waters.gro except the new waters (make
sure you remove waters that were in initial.gro)
4. use vi to edit keepbyz.pl
   - upperz and lowerz variables as you please
   - sol to the name of your solvent molecule
5. run keepbyz.pl on new_waters.gro
   ./keepbyz new_waters.gro > keep_these_waters.gro
6. tail -1 initial.gro > last_line.gro
7. head -$(expr $(cat initial.gro | wc -l | awk '{print $1}') - 1 )  
initial.gro

* not_last_line.gro
*8. cat not_last_line.gro new_waters.gro last_line.gro >  
new_system.gro 9. editconf -f new_system.gro -o  
new_system_sequential_numbers.gro


I could successfully go through steps 1 to 6. However, I am unable to
execute step 7. In particular, what is "expr" in step7? Moreover, The file,
"

keep_these_waters.gro" that I have generated in step 5, where do I use
it subsequently?

I shall remain obliged for your kind answer.
Best Regards.


--
Dr. Tusar Bandyopadhyay
Theoretical Chemistry Section,
Chemistry Group
BARC, Trombay

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

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

[chris . neale]

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


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


Chris.

-- original message --


Dear all,

I would like to calculate the potential of mean force between two molecules
in aqueous solution using the pull code. For a start i performed a number
of calulations with a couple of very simple model systems with settings
loosely based on the example in the gmx tutorial section (see mdp file  
included

at the bottom). My box is 8x8x8 nm**3, pmf is calculated from some r_min up to
approx 3.8 nm. Performing simulations with two t-buntanol molecules in  
implicit

solvent, anlyzing the output (x and f*.xvg files) with g_wham, i get rather
unexpected results:

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

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

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

Although what I see suggests that with "pulldim  Y Y Y" the RDF just doesn't
normalized properly, this issue seems to be more involved:
Looking at the average forces and positions (the content of the f*xvg  
and x*xvg

files)
http://www.brunsteiner.net/tbut-f.gif
http://www.brunsteiner.net/tbut-z.gif
suggests that there's already something wrong in the mdrun output
meaning that the problem is further upstream and not connected to anything
done by g_wham.

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

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

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

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

cheers,
Michael



=
mdp file:

integrator  = sd; this is better than md/NHT for systems with very
few DOFs
tau_t   = 2.0 2.0
ref_t   = 290 290   ; a bit lower than 298 since sd with default
parameters
; has problems controlling the temperature.
tc_grps = p1 p2 ; also tried System here, no difference
;
dt  = 0.002
tinit   = 0
nsteps  = 10; played with that, results seem to have converged
nstxout = 1000
nstvout = 0
nstfout = 1000
nstenergy   = 100
;
constraint_algorithm = lincs
constraints  = hbonds
continuation = no
;
comm_mode   = Linear
nstcomm = 1
pbc = xyz   ; also tried "pbc = no" same result
nstlist = 10
ns_type = grid
rlist   = 4.0
rcoulomb= 4.0
rvdw= 4.0
;
coulombtype = Shift
vdwtype = Shift
epsilon_r   = 80
;
pull= umbrella
pull_geometry   = distance
pull_dim= N N Y   ; or "Y Y Y"
pull_start  = yes
pull_ngroups= 1
pull_group0 = p1
pull_group1 = p2
pull_rate1  = 0.0
pull_k1 = 500

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

[chris . neale]

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


Chris.

-- original message --

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

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


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

to enhance sampling...

best regards,
Jianguo





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

[gmx-users] test on g_density

[chris . neale]

Could it be this:

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

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


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


Chris.

-- original message --


Hi all.

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

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

The generated trajectories have 12501 and 12500 steps respectively.

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

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

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

Thank you!

Javier
-- next part --


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

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

[chris . neale]

Michael: your mdp options indicated using US and my stance is that,  
using US, you are incorrect to say that "And you will get the SAME  
result if you do this calculation in one or in three dimensions  
(pulldim NNY or pulldim YYY)". Nevertheless, I'm a little perturbed by  
your call out for help from somebody else "in the know", so I'll leave  
you at this point. Good luck.


Chris.

-- original message --

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

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

chris.neale at utoronto.ca wrote:

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

as per an RDF.


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

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

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

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

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

is, of course, nothing else but the original LJ potential. And you  
will get the

SAME result
if you do this calculation in one or in three dimensions (pulldim NNY  
or pulldim

YYY)
so where does the entropy and the correction term in eqn 6 come in here??

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

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

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

this
distance ... but then this is just wild guess of mine ... anyway, that would
mean

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

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

Yours truly (and maybe some other people) would really appreciate if  
somebody in

the know
would clarify that!

thanks in advance!

Michael





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

[gmx-users] Questions on combination of Berger's united-atom force field for lipid and OPLS_AA force field for protein

[chris . neale]

Your first question is addressed in the following document:

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

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


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


Chris.

-- original message --

Dear All:

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

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

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

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

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

4). make a topology file like this:

My questions are as follows:

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

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

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

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

Regards

Ruo-Xu Gu


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

[gmx-users] using one thermostat for entire system ??

[chris . neale]

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


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


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


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


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


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

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

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


Chris.

-- original message --

Thanks Mark for the advice.

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

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

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

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

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

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

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

Any comments welcome

Thanks

Jenny






Quoting Mark Abraham :


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

Hi,

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


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


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


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

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

[chris . neale]

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


Chris.

-- original message --

Hi,

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

Thank you for your help.

Susana

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



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

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

[gmx-users] g_wham PMF profile

[chris . neale]

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


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


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


Chris.

-- original message --

Thank you!

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

best,

Anna


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

[gmx-users] Pulling

[chris . neale]

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


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


Chris.

-- original message --

Dear All
Afew question about Pulling in Umberella Sampling

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

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

Thanks in advance
-- next part --


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

[gmx-users] Umberella sampling

[chris . neale]

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

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


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


Chris.

-- original message --

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

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

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

Thanks in advance



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

[gmx-users] Pulling

[chris . neale]

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


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


... snip ...


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


sure, but you will have more convergence problems.

Chris.

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


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

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

Chris.

-- original message --

Dear All
Afew question about Pulling in Umberella Sampling

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

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

Thanks in advance
-- next part --



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

[gmx-users] Umberella sampling

[chris . neale]

Mohsen:

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


Chris.

-- original message --

Dear chris
Thanks for your reply

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

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

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


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

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

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

Chris.

-- original message --

Dear All



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

[gmx-users] pick in pull.xvg

[chris . neale]

I did a pulling simulation for protein-drug system.


This could mean almost anything. Please provide details.

Chris.

-- original message --

Dear All

I did a pulling simulation for protein-drug system.

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

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

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

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



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

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

[gmx-users] plot of force via time

[chris . neale]

It depends on what you are trying to do.

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


Chris.

-- original message --

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

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

Dear All

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

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

Thanks ina dvance


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

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

[gmx-users] g_sham

[chris . neale]

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


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


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


Chris.

-- original message --

Dear All

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

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


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

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

[gmx-users] compressing a box of water droplets into a homogeneous solution of liquid water

[chris . neale]

Dear users:

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


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


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


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


Thank you,
Chris.

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

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

[gmx-users] Simulation of slow folding proteins

[chris . neale]

Dear Bharat:

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


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


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


Chris.

-- original message --

Hi,

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

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

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

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

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


--
Bharat
Ph.D. Candidate
Room No. : 7202A, 2nd Floor
Biomolecular Engineering Laboratory
Division of Chemical Engineering and Polymer Science
Pusan National University
Busan -609735
South Korea
Lab phone no. - +82-51-510-3680, +82-51-583-8343
Mobile no. - 010-5818-3680
E-mail : monu46010 at yahoo.com


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

[gmx-users] compressing a box of water droplets into a homogeneous solution of liquid water

[chris . neale]

Thanks Patrick and Andre!

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


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


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


Chris.

-- original message --

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

Patrick

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

[Hide Quoted Text]
Dear users:

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

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

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

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

Thank you,
Chris.

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

[gmx-users] Umbrella sampling

[chris . neale]

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


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


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


Chris.

-- original message --


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

Gavin
XAvier Periole wrote:


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

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

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


Hi Xavier

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

Gavin

XAvier Periole wrote:


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


Hi All

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

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

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

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

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

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


Many thanks

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




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




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

[gmx-users] Umbrella Sampling

[chris . neale]

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


Chris.

-- original message --

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

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

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


--

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

[gmx-users] Umbrella Sampling

[chris . neale]

your comment:

which should be centred around 0.80nm

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


chris.

-- original message --

Cheers Chris

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

Gavin

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

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

Chris.

-- original message --

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


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

[gmx-users] Umbrella Sampling

[chris . neale]

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


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


Chris.


Cheers Chris

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

Gavin
chris.neale at utoronto.ca wrote:

your comment:

which should be centred around 0.80nm

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

chris.

-- original message --

Cheers Chris

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

Gavin

chris.neale at utoronto.ca wrote:

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

Chris.

-- original message --

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



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

[gmx-users] FEP and loss of performance

[chris . neale]

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


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


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


Chris.

-- original message --


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

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

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

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

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

Thanks


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

[gmx-users] minimization and simulation problems

[chris . neale]

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


-- original message --


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

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

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

---

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[gmx-users] FEP and loss of performance

[Chris Neale]

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


There must be some extra communication?

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


Chris.


Luca Bellucci wrote:

/  Hi Chris,

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

-Justin


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

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

[gmx-users] minimization and simulation problems

[Chris Neale]

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


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


Chris.

Quoting politr at fh.huji.ac.il 
:


Dear gromacs users,

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

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

/

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

[gmx-users] FEP and loss of performance

[Chris Neale]

>> Dear Chris and Justin


/  Thank you for your precious suggestions

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

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

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


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

Chris.


-Justin


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

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

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

[gmx-users] FEP and loss of performance

[chris . neale]

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


Good luck,
Chris.

-- original message --


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


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

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

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

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

Chris.


-Justin

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

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


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

[gmx-users] converging value

[chris . neale]

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


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


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


Chris.

-- original message --

Dear All

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

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

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

[gmx-users] rationale behind tcoupling ligand with protein instead of with SOL

[chris . neale]

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


Chris.

-- original message --

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

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

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

[gmx-users] Unexpected results arising from T- and P-coupling methods

[chris . neale]

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

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

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


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

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


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


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


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

Chris.

-- original message --

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


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

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



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

[gmx-users] Different TI free energy values in 4.0.7 and 4.5.3

[chris . neale]

Dear Yan:

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


Chris.

-- original message --

Dear Gromacs users,

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

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

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

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

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

  Does anyone know the reason for the different TI free energy values in
those different version of Gromacs?

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

[gmx-users] Different TI free energy values in 4.0.7 and 4.5.3

[chris . neale]

Based on your use of:

sc-sigma = 0.3

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


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

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


-- original message --

Dear Gromacs users,

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

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


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


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


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


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

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

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

   The values are converged by block averaging.

   Regards,
   Yan

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


Dear Yan:

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

Chris.

-- original message --

Dear Gromacs users,

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

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

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

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

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

 Does anyone know the reason for the different TI free energy values in
those different version of Gromacs?

--



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

[gmx-users] Cylinder Pulling output file pullx.xvg wron

[chris . neale]

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


Chris.

-- original message --

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

Thanks for you help.
Kun

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

Here is the output:

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

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

[gmx-users] Use different fudgeLJ and fudgeQQ values in simulations with the AMBER force field

[chris . neale]

Stephane:

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


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


1. set the fudgeLJ to 1.0 and fudgeQQ to 0.16

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


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


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


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


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


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

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

Chris.

-- original message --

Dear gmx users,

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

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

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

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

Thank you in advance for your response


Stephane





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

[gmx-users] lipid ff decisions

[chris . neale]

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


Thanks,
Chris.

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

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

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

[Chris Neale]

Dear Users:

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


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

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


The same results exist for at least 4.0.5 and 4.5.3.

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

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

[gmx-users] g_select is pretty fantastic

[chris . neale]

Dear g_select developers:

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

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


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


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


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

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


Pretty cool.

Thanks again,
Chris.

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

[gmx-users] free energy

[chris . neale]

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


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


Chris.

-- original message --

Hi all

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

Cheers

Gavin


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

[gmx-users] free energy

[chris . neale]

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


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


Chris.

-- original message --


Hi Chris

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

Gavin


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

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

[gmx-users] free energy

[chris . neale]

Dear Gavin:

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


Chris.

-- original message --

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

Gavin

chris.neale at utoronto.ca wrote:

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

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

Chris.

-- original message --


Hi Chris

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

Gavin



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

[gmx-users] difference of protein dynamics between NPT and NVE ensemble?

[chris . neale]

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


-- original message --

Does anyone have experience of doing NVE simulation? I wish to know, for a
system with box size 6.1*6.1*6.1nm, whether the microcanonical ensemble for
a NVE simulation will converge with the canonical ensemble from a NPT
simulation?



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

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

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

[Chris Neale]

Dear Users:

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


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


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


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


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


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

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

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

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


Parallel run - timing based on wallclock.

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

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


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

#

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


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


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

  S  C  A  M  O  R  G

:-)  VERSION 4.0.5  (-:


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

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

   :-)  gmxcheck  (-:

Option Filename  Type Description

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

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

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

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

Checking file md1.cpt

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

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

---

"Confirmed" (Star Trek)

##

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

[Chris Neale]

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

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

Dear Users:

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


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


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


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


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


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

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

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

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


Parallel run - timing based on wallclock.

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

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


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

#

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


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


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

  S  C  A  M  O  R  G

:-)  VERSION 4.0.5  (-:


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

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

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

   :-)  gmxcheck  (-:

Option Filename  Type Description

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

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

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

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

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

Checking file md1.cpt

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

[gmx-users] constraints

[Chris Neale]

generally, look at mdout.mdp from grompp

-- original message --

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

What is default

constraints=none

constraints=all bonds


NIlesh




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

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

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

[chris . neale]

Dear Users:

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


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


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


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


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


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


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

[ system ]
; name
tables

[ molecules ]
; name  number
MOL_A  1
MOL_B  1

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


[ atomtypes ]
 aaa   AA  8 15.999400.0   A0.0  0.0

[ moleculetype ]
; molname   nrexcl
MOL_A 1

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

[ moleculetype ]
; molname   nrexcl
MOL_B 1

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

 And the .gro file:

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

### And my .mdp file is

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



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


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

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

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


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

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


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


Thank you for your assistance,
Chris.



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

[gmx-users] beware gen-vel=no

[chris . neale]

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

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


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


Chris.

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

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

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

[chris . neale]

Dear users:

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


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


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


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


Thank you,
Chris.

-- original message --

Dear Users:

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

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

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

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

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

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

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

[ system ]
; name
tables

[ molecules ]
; name  number
MOL_A  1
MOL_B  1

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

[ atomtypes ]
  aaa   AA  8 15.999400.0   A0.0  0.0

[ moleculetype ]
; molname   nrexcl
MOL_A 1

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

[ moleculetype ]
; molname   nrexcl
MOL_B 1

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

 And the .gro file:

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

### And my .mdp file is

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



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

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

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

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

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

[chris . neale]

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


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


Chris.

-- original message --

Hi,

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

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

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

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

Kind regards,
Luke



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

[gmx-users] How to get angle distribution between two tyrosine stacking residues

[chris . neale]

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

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

-- original message --

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


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

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

[gmx-users] problem during energy minimization

[chris . neale]

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

   integrator = md

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


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

   ;define= -DEFLEXIBLE

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

Chris.

-- original message --



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

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

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

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





*OUTPUT OF MDRUN*

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

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


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

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

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

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

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

regards,
sree
-- next part --


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

[gmx-users] binaries for mac os x 10.6

[chris . neale]

Dear Users:

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


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


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


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


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


Thank you very much for your time,
Chris.

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

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

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

[chris . neale]

Dear Dimitar:

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


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


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


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


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


Chris.


 Here is an example:

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

 Submitted by:

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

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

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

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

fi
 =


This script is not using mdrun -append.



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



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



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

as usual submitted by:

qsub -N  myjob.q





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



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

  


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


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


 Last output files from restarts:

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

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

[gmx-users] Announcement: Large biomolecule benchmark report

[chris . neale]

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


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


Chris.

-- original message --

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

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

Thanks again,
Hannes.


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

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

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

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

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

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

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

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

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

[gmx-users] crazy temperatures

[chris . neale]

I disagree.

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


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


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


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


Chris.

-- original message --

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


Warren Gallin

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

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

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

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


Mark
--


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

[gmx-users] Umbrella sampling along Radius of gyration

[chris . neale]

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


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


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

Chris.

-- original message --

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

Thanks
Sanku



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

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

[gmx-users] Not able to continue with Equilibration

[chris . neale]

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


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

Chris.

-- original message --

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

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



2012/3/29 Rausch, Felix :

[Hide Quoted Text]
Hello again,

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

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

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

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

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

Does anyone have any suggestions?
Thanks
Francesca


MINIMIZATION OUTPUT

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

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


MINIMIZATION.MDP

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

; 

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

[chris . neale]

Dear Users:

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


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

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

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

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


-- Configuring incomplete, errors occurred!

###

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

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


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


Thank you,
Chris.

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

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

Bert
Thu, 31 Mar 2011 10:15:22 -0700

Dear all,

When GMX is builded with cmake, how to compile the template.c? I used
to make the template.c by the command "make -f MakefileXXX", but it
can not work in version 4.5. Thanks for your suggestions.

Best regards,
Bert
--




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

[gmx-users] on the Umbrella Sampling on-line tutorial

[chris . neale]

Dear Acoot:

I'll reply to general topics, not to the tutorial in particular.

If the opening is not large enough to allow the peptide to exit, then  
the surrounding protein will need to change its conformation to permit  
unbinding. This can happen, but you need (a) to have sufficiently long  
sampling times to permit the required conformational change and (b)  
starting structures in which the peptide is near the umbrella center  
and do not crash during simulation.


The US method is always valid, but it is not always the best choice.  
You might try using the free energy code to implement a thermodynamic  
cycle. Nevertheless, one imagines that the protein does actually need  
to open up during peptide binding and so you will still need to sample  
protein opening/closing in order to obtain equilibrium (i.e. correct)  
binding free energies because you need to sample the unbound protein  
at equilibrium. That is to say that a thermodynamic cycle may appear  
converged when in fact it is not, because you have not converged the  
unbound state of the protein. In any event, be careful with your  
convergence analysis.


This would be a good system for which to attempt both US and  
double-decoupling approaches and use the results of each to ensure  
that you are not missing some important conformational states.


That said, I highly doubt that it is possible to converge the free  
energies of induced-fit peptide-protein binding with any available  
atomistic computational method using contemporary computational  
resources. The lower bound of required sampling times is certainly on  
the order of 10 us per umbrella and I bet that the actual value is a  
few orders of magnitude larger. I have less experience with the free  
energy code than I do with US, but I suspect that the required  
sampling times are also very long for a system like this.


Chris.

-- original message --

Dear All,

I planned to use the method introduced in the Umbrella Sampling  
on-line tutorial of Justin Lemkul  
(http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html).


But if a peptide is surrounded by a protein, which means the opening  
of the protein-complex is not large enough to allow the peptide to  
leave the protein without significantly breaking the conformation of  
the protein in the protein-peptide complex, is the Umbrella Sampling  
method still valid for the binding energy calculation?


Will you please also show me in which part of the tutorial the  
direction of pull-apart is defined? We should process it in a  
direction the peptide can leave the protein, not the direction protein  
will bind the peptide much strongly.


I am looking forward to getting a reply from you.

Cheers,

Acoot


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

[gmx-users] further discussion on the mdrun -append function

[chris . neale]

Dear Acoot:

You should be able to answer this one yourself. Moreover, you are  
doing yourself a disservice by relying on the mailing list to do your  
work for you because you will eventually need to learn how to find  
answers to these things on your own.


Please remember the following:

1. use a title that matches your question, starting a new thread if  
you have a new question


2. it may sound harsh, but questions like this one in which you have  
obviously not even tried to find the answer yourself for more than a  
couple of minutes tend to annoy the very people that you may want help  
from later on with other issues.


3. whenever you post, please show how you have tried to solve the  
problem yourself. You may find that in the process of writing such a  
question you end up solving the problem yourself.


4. read the manual.

Chris.

-- original message --


Thanks for the detailed explaination.

Will you please explain the function of "-n index.ndx " in "grompp -f  
md.mdp -n index.ndx -p topol.top -c minimized.gro -o md-1ns.tpr" with  
some specific examples?


Cheers,

Acoot


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

[gmx-users] Different results from identical tpr after MD

[chris . neale]

Dear Acoot:

The idea of convergence is this: start a large number of simulations  
from different conformations, analyze some quantity over time in each  
simulation, and when the deviation of the average value of that  
quantity from each separate simulation is less than the time-variance  
within individual simulations, then you can imply that the simulations  
have converged -- that is, the results of independent simulations  
which started as different are now similar.


There is a huge body of work that uses a single simulations and  
evaluates its so-called convergence using some assumptions and special  
methods. That can also be very useful, but I find it informative to  
think of "convergence" in its standard non-scientific dictionary  
definition as the coming together of previously disparate things.


For simulations, my working definition is this: a set of simulations  
has converged the value of some variable when the simulations were  
initiated from sufficiently distinct conformational basins and then,  
over time, the ensemble distribution of the time-averages of the  
specified variable has a variance that is the same as the mean  
time-averaged variation within independent simulations. The weak point  
here is the part about "sufficiently distinct conformations", but I am  
not sure that this can be stated less vaguely in the general case.


Chris.

-- original message --

Hi Justin,

Can you give me your definition of converged MD and unconverged MD?

Cheers,

Acoot


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

[gmx-users] mdrun -nosum still complains that 15 % of the run time was spent communicating energies

[Chris Neale]

Hello,

I have been running simulations on a larger number of processors 
recently and am confused about the message regarding -nosum that occurs 
at the end of the .log file. In this case, I have included the -nosum 
option to mdrun and I still get this warning (gromacs 4.0.4).


My command was:
mpirun -np $(wc -l $PBS_NODEFILE | gawk '{print $1}') -machinefile 
$PBS_NODEFILE /scratch/cneale/exe/intel/gromacs-4.0.4/exec/bin/mdrun 
-deffnm test -nosum -npme 128


#

To confirm that I am asking mdrun for -nosum, to stderr I get:
...
Option   Type   Value   Description
--
-[no]h   bool   no  Print help info and quit
-niceint0   Set the nicelevel
-deffnm  string testSet the default filename for all file options
-[no]xvgrbool   yes Add specific codes (legends etc.) in the output
   xvg files for the xmgrace program
-[no]pd  bool   no  Use particle decompostion
-dd  vector 0 0 0   Domain decomposition grid, 0 is optimize
-npmeint128 Number of separate nodes to be used for PME, -1
   is guess
-ddorder enum   interleave  DD node order: interleave, pp_pme or 
cartesian

-[no]ddcheck bool   yes Check for all bonded interactions with DD
-rdd real   0   The maximum distance for bonded interactions 
with

   DD (nm), 0 is determine from initial coordinates
-rconreal   0   Maximum distance for P-LINCS (nm), 0 is estimate
-dlb enum   autoDynamic load balancing (with DD): auto, no 
or yes

-dds real   0.8 Minimum allowed dlb scaling of the DD cell size
-[no]sum bool   no  Sum the energies at every step
-[no]v   bool   no  Be loud and noisy
-[no]compact bool   yes Write a compact log file
-[no]seppot  bool   no  Write separate V and dVdl terms for each
   interaction type and node to the log file(s)
-pforce  real   -1  Print all forces larger than this (kJ/mol nm)
-[no]reprod  bool   no  Try to avoid optimizations that affect binary
   reproducibility
-cpt real   15  Checkpoint interval (minutes)
-[no]append  bool   no  Append to previous output files when continuing
   from checkpoint
-[no]addpart bool   yes Add the simulation part number to all output
   files when continuing from checkpoint
-maxhreal   -1  Terminate after 0.99 times this time (hours)
-multi   int0   Do multiple simulations in parallel
-replex  int0   Attempt replica exchange every # steps
-reseed  int-1  Seed for replica exchange, -1 is generate a seed
-[no]glasbool   no  Do glass simulation with special long range
   corrections
-[no]ionize  bool   no  Do a simulation including the effect of an X-Ray
   bombardment on your system
...



And the message at the end of the .log file is:
...
   D O M A I N   D E C O M P O S I T I O N   S T A T I S T I C S

av. #atoms communicated per step for force:  2 x 3376415.3
av. #atoms communicated per step for LINCS:  2 x 192096.6

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

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

NOTE: 7.9 % performance was lost due to load imbalance
 in the domain decomposition.

NOTE: 10.0 % performance was lost because the PME nodes
 had less work to do than the PP nodes.
 You might want to decrease the number of PME nodes
 or decrease the cut-off and the grid spacing.


R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G

Computing: Nodes Number G-CyclesSeconds %
---
Domain decomp.   256 51  337.551  131.2 0.7
Send X to PME256501   59.454   23.1 0.1
Comm. coord. 256501  289.936  112.7 0.6
Neighbor search  256 51 1250.088  485.9 2.8
Force25650116105.584 6259.935.4
Wait + Comm. F   256501 2441.390  948.9 5.4
PME mesh 128501 5552.336 2158.112.2
Wait + Comm. X/F 128501 9586.486 3726.121.1
Wait + Recv. PME F   256501  459.752  178.7 1.0
Write traj.  256  2  223.993   87.1 0.5
Update   256501  777.618  302.2 1.7
Constraints  256   1002 1223.093  475.4 2.7
Comm. energies   256 51 7011.309 2725.115.4
Rest  

[gmx-users] Re: mdrun -nosum still complains that 15 % of the run time was spent communicating energies

[Chris Neale]

I have now tested with and without -nosum and it appears that the option 
is working (see 51 vs. 501 Number of communications) but that the total 
amount of time communicating energies didn't go down by very much. Seems 
strange to me. Anybody have any ideas if this is normal?


At the very least, I suggest adding an if statement to mdrun so that it 
doesn't output the -nosum usage note if the user did in fact use -nosum 
in that run.



Without using -nosum:

   R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G

Computing: Nodes Number G-CyclesSeconds %
---
...
Write traj.  256  2  233.218   93.7 0.5
Update   256501  777.511  312.5 1.7
Constraints  256   1002 1203.894  483.9 2.7
Comm. energies   256501 7397.995 2973.916.5
Rest 256 128.058   51.5 0.3
---
Total384   44897.46818048.0   100.0
---

NOTE: 16 % of the run time was spent communicating energies,
 you might want to use the -nosum option of mdrun


   Parallel run - timing based on wallclock.

  NODE (s)   Real (s)  (%)
  Time: 47.000 47.000100.0
  (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:  13485.788712.634  1.842 13.029
Finished mdrun on node 0 Mon Jul 20 12:53:41 2009

#

And using -nosum:

   R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G
Computing: Nodes Number G-CyclesSeconds %
---
...
Write traj.  256  2  213.521   83.3 0.5
Update   256501  776.606  303.0 1.8
Constraints  256   1002 1200.285  468.2 2.7
Comm. energies   256 51 6926.667 2702.115.6
Rest 256 127.503   49.7 0.3
---
Total384   44296.67017280.0   100.0
---

NOTE: 16 % of the run time was spent communicating energies,
 you might want to use the -nosum option of mdrun


   Parallel run - timing based on wallclock.

  NODE (s)   Real (s)  (%)
  Time: 45.000 45.000100.0
  (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:  14084.547744.277  1.924 12.475

#

Thanks,
Chris.

Chris Neale wrote:

Hello,

I have been running simulations on a larger number of processors 
recently and am confused about the message regarding -nosum that 
occurs at the end of the .log file. In this case, I have included the 
-nosum option to mdrun and I still get this warning (gromacs 4.0.4).


My command was:
mpirun -np $(wc -l $PBS_NODEFILE | gawk '{print $1}') -machinefile 
$PBS_NODEFILE /scratch/cneale/exe/intel/gromacs-4.0.4/exec/bin/mdrun 
-deffnm test -nosum -npme 128


#

To confirm that I am asking mdrun for -nosum, to stderr I get:
...
Option   Type   Value   Description
--
-[no]h   bool   no  Print help info and quit
-niceint0   Set the nicelevel
-deffnm  string testSet the default filename for all file options
-[no]xvgrbool   yes Add specific codes (legends etc.) in the 
output

   xvg files for the xmgrace program
-[no]pd  bool   no  Use particle decompostion
-dd  vector 0 0 0   Domain decomposition grid, 0 is optimize
-npmeint128 Number of separate nodes to be used for 
PME, -1

   is guess
-ddorder enum   interleave  DD node order: interleave, pp_pme or 
cartesian

-[no]ddcheck bool   yes Check for all bonded interactions with DD
-rdd real   0   The maximum distance for bonded 
interactions with
   DD (nm), 0 is determine from initial 
coordinates
-rconreal   0   Maximum distance for P-LINCS (nm), 0 is 
estimate
-dlb enum   autoDynamic load balancing (with DD): auto, no 
or yes
-dds real   0.8 Minimum allowed dlb scaling of the DD cell 
size

-[no]sum bool   no  Sum the energies at every step
-[no]v   bool   no  Be loud and noisy
-[no]compact bool   yes Write a compact log file
-[no]seppot  bool   no  Write separate V and dVdl terms for each
   interaction type and node to the log file(s)
-pforce  real   -1  Print all forces larger than this (kJ/mol nm)
-[no]reprod  bool   

[gmx-users] Re: mdrun -nosum still complains that 15 % of > the runtime was spent communicating energies

[chris . neale]

Thanks mark, I'll respond inline.


Chris Neale wrote:

I have now tested with and without -nosum and it appears that the option
is working (see 51 vs. 501 Number of communications) but that the total
amount of time communicating energies didn't go down by very much. Seems
strange to me. Anybody have any ideas if this is normal?


Seems strange, but perhaps a 45-second test is not sufficiently long to
demonstrate suitable scaling.


Agreed, although I find it a good jumping off point, especially now  
that we need to optimize -npme. If I use quick tests to narrow down  
the range of nodes/npme that is going to scale the best then I fine  
tune it with longer scaling tests.


There's no discussion in the 4.0.5 release

notes of a relevant change to -nosum, but there has been a change:
http://oldwww.gromacs.org/content/view/181/132/.


Thanks, I did see this but don't think that it is related to this  
issue, which I have now confirmed in both 4.0.4 and 4.0.5.





At the very least, I suggest adding an if statement to mdrun so that it
doesn't output the -nosum usage note if the user did in fact use -nosum
in that run.


Without using -nosum:

   R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G

Computing: Nodes Number G-CyclesSeconds %
---
...
Write traj.  256  2  233.218   93.7 0.5
Update   256501  777.511  312.5 1.7
Constraints  256   1002 1203.894  483.9 2.7
Comm. energies   256501 7397.995 2973.916.5
Rest 256 128.058   51.5 0.3
---
Total384   44897.46818048.0   100.0
---

NOTE: 16 % of the run time was spent communicating energies,
 you might want to use the -nosum option of mdrun


   Parallel run - timing based on wallclock.

  NODE (s)   Real (s)  (%)
  Time: 47.000 47.000100.0
  (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:  13485.788712.634  1.842 13.029
Finished mdrun on node 0 Mon Jul 20 12:53:41 2009

#

And using -nosum:

   R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G
Computing: Nodes Number G-CyclesSeconds %
---
...
Write traj.  256  2  213.521   83.3 0.5
Update   256501  776.606  303.0 1.8
Constraints  256   1002 1200.285  468.2 2.7
Comm. energies   256 51 6926.667 2702.115.6
Rest 256 127.503   49.7 0.3
---
Total384   44296.67017280.0   100.0
---

NOTE: 16 % of the run time was spent communicating energies,
 you might want to use the -nosum option of mdrun


   Parallel run - timing based on wallclock.

  NODE (s)   Real (s)  (%)
  Time: 45.000 45.000100.0
  (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:  14084.547744.277  1.924 12.475

#####

Thanks,
Chris.

Chris Neale wrote:

Hello,

I have been running simulations on a larger number of processors
recently and am confused about the message regarding -nosum that
occurs at the end of the .log file. In this case, I have included the
-nosum option to mdrun and I still get this warning (gromacs 4.0.4).

My command was:
mpirun -np $(wc -l $PBS_NODEFILE | gawk '{print $1}') -machinefile
$PBS_NODEFILE /scratch/cneale/exe/intel/gromacs-4.0.4/exec/bin/mdrun
-deffnm test -nosum -npme 128


Perhaps assigning the result of this to a variable and printing it
before executing it would help confirm that -nosum really was there.


I am not sure what you mean... the while line as a variable? I'm  
pretty sure that -nosum is there.




Your mdrun output from your first email was...


#

To confirm that I am asking mdrun for -nosum, to stderr I get:
...
Option   Type   Value   Description
--
-[no]h   bool   no  Print help info and quit
-niceint0   Set the nicelevel
-deffnm  string testSet the default filename for all file options
-[no]xvgrbool   yes Add specific codes (legends etc.) in the
output
   xvg files for the xmgrace program
-[no]pd  bool   no  Use particle decompostion
-dd  vector 0 0 0   Domain decomposition grid, 0 is optimize
-npmeint128 Number 

[gmx-users] Problem compiling gromacs-4.0.5 against mvapich2-1.4rc1

[chris . neale]

Hi Jerry,

did you have any luck solving this? I have run into the same problem  
compiling gromacs-4.0.5 with MVAPICH2-1.4rc1 using the intel compiler.  
I was able to successfully compile gromacs-4.0.5 using openmpi-1.3.2  
and the intel compiler. In both cases, the MPI libraries were compiled  
using the intel compiler v11.0.081.


make errors out with:

...
mpicc -DHAVE_CONFIG_H -I. -I../../src  -I../../include  
-DGMXLIBDIR=\"/scratch/cneale/exe/intel/gromacs-4.0.5/exec/share/top\"  
-I/scratch/cneale/exe/intel/fftw-3.1.2/exec/include  
-I/scinet/gpc/mpi/mvapich2/1.4rc1-3378_intel-v11.0_ofed/include  
-I/scinet/gpc/mpi/mvapich2/1.4rc1-3378_intel-v11.0_ofed/lib  
-I/usr/local/include  -O3 -tpp7 -axW -ip -w -funroll-all-loops  
-std=gnu99 -MT grompp.o -MD -MP -MF .deps/grompp.Tpo -c -o grompp.o  
grompp.c

grompp.c(983): (col. 5) remark: LOOP WAS VECTORIZED.
mv -f .deps/grompp.Tpo .deps/grompp.Po
/bin/sh ../../libtool --tag=CC   --mode=link mpicc  -O3 -tpp7 -axW -ip  
-w -funroll-all-loops -std=gnu99   
-L/scratch/cneale/exe/intel/fftw-3.1.2/exec/lib   -o grompp grompp.o  
libgmxpreprocess_mpi.la ../mdlib/libmd_mpi.la ../gmxlib/libgmx_mpi.la   
-lnsl -lfftw3f
mpicc -O3 -tpp7 -axW -ip -w -funroll-all-loops -std=gnu99 -o grompp  
grompp.o  -L/scratch/cneale/exe/intel/fftw-3.1.2/exec/lib  
./.libs/libgmxpreprocess_mpi.a ../mdlib/.libs/libmd_mpi.a  
/scratch/cneale/exe/intel/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a  
../gmxlib/.libs/libgmx_mpi.a -lnsl  
/scratch/cneale/exe/intel/fftw-3.1.2/exec/lib/libfftw3f.a -lm
/scinet/gpc/mpi/mvapich2/1.4rc1-3378_intel-v11.0_ofed/lib/libmpich.a(ibv_channel_manager.o):(.bss+0x18): multiple definition of  
`debug'
/scratch/cneale/exe/intel/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a(gmx_fatal.o):(.bss+0x530): first defined  
here

make[3]: *** [grompp] Error 1
make[3]: Leaving directory  
`/scratch/cneale/exe/intel/gromacs-4.0.5/src/kernel'

make[2]: *** [all-recursive] Error 1
make[2]: Leaving directory `/scratch/cneale/exe/intel/gromacs-4.0.5/src'
make[1]: *** [all] Error 2
make[1]: Leaving directory `/scratch/cneale/exe/intel/gromacs-4.0.5/src'
make: *** [all-recursive] Error 1


could this 'debug' be the one int calls like this? I at first looked  
for a debug() function but didn't find any.


$ grep debug src/gmxlib/*.c
...
main.c:  if (debug) {
main.c:fprintf(debug,"This is simulation %d",cr->ms->sim);
main.c:  fprintf(debug,", local number of nodes %d, local nodeid %d",
main.c:fprintf(debug,"\n\n");
...

Thanks,
Chris

#3

-- original message --

Hello all,

I have built MVAPICH2-1.4rc1 with Infiniband support on a CentOS 5.3 box and
am trying to build gromacs-4.0.5 with MPI support.

make is failing with the following error message:

mpicc -O3 -fomit-frame-pointer -finline-functions -Wall -Wno-unused
-funroll-all-loops -std=gnu99 -o grompp grompp.o
./.libs/libgmxpreprocess_mpi.a -L/usr/lib64 ../mdlib/.libs/libmd_mpi.a
/<...snip..>/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a
../gmxlib/.libs/libgmx_mpi.a -lxml2 -lnsl -lfftw3f -lm -lX11

/usr/local/lib/libmpich.a(ibv_channel_manager.o):(.bss+0x10): multiple
definition of `debug'
<..snip..>/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a(gmx_fatal.o):(.bss+0x0):
first defined here
collect2: ld returned 1 exit status
make[3]: *** [grompp] Error 1

It looks like 'debug' is conflicting in both MVAPICH2 and in Gromacs.

Any suggestions?

Thanks,

Jerry.
Research Systems Administrator
SBGrid


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

[gmx-users] Re: Problem compiling gromacs-4.0.5 against mvapich2-1.4rc1

[chris . neale]

Found it.

Jerry got the solution here:  
http://mail.cse.ohio-state.edu/pipermail/mvapich-discuss/2009-July/002403.html


Thanks Jerry,
Chris.

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


Hi Jerry,

did you have any luck solving this? I have run into the same problem
compiling gromacs-4.0.5 with MVAPICH2-1.4rc1 using the intel compiler.
I was able to successfully compile gromacs-4.0.5 using openmpi-1.3.2
and the intel compiler. In both cases, the MPI libraries were compiled
using the intel compiler v11.0.081.

make errors out with:

...
mpicc -DHAVE_CONFIG_H -I. -I../../src  -I../../include
-DGMXLIBDIR=\"/scratch/cneale/exe/intel/gromacs-4.0.5/exec/share/top\"
-I/scratch/cneale/exe/intel/fftw-3.1.2/exec/include
-I/scinet/gpc/mpi/mvapich2/1.4rc1-3378_intel-v11.0_ofed/include
-I/scinet/gpc/mpi/mvapich2/1.4rc1-3378_intel-v11.0_ofed/lib
-I/usr/local/include  -O3 -tpp7 -axW -ip -w -funroll-all-loops
-std=gnu99 -MT grompp.o -MD -MP -MF .deps/grompp.Tpo -c -o grompp.o
grompp.c
grompp.c(983): (col. 5) remark: LOOP WAS VECTORIZED.
mv -f .deps/grompp.Tpo .deps/grompp.Po
/bin/sh ../../libtool --tag=CC   --mode=link mpicc  -O3 -tpp7 -axW -ip
-w -funroll-all-loops -std=gnu99
-L/scratch/cneale/exe/intel/fftw-3.1.2/exec/lib   -o grompp grompp.o
libgmxpreprocess_mpi.la ../mdlib/libmd_mpi.la ../gmxlib/libgmx_mpi.la
-lnsl -lfftw3f
mpicc -O3 -tpp7 -axW -ip -w -funroll-all-loops -std=gnu99 -o grompp
grompp.o  -L/scratch/cneale/exe/intel/fftw-3.1.2/exec/lib
./.libs/libgmxpreprocess_mpi.a ../mdlib/.libs/libmd_mpi.a
/scratch/cneale/exe/intel/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a
../gmxlib/.libs/libgmx_mpi.a -lnsl
/scratch/cneale/exe/intel/fftw-3.1.2/exec/lib/libfftw3f.a -lm
/scinet/gpc/mpi/mvapich2/1.4rc1-3378_intel-v11.0_ofed/lib/libmpich.a(ibv_channel_manager.o):(.bss+0x18): multiple definition   
of

`debug'
/scratch/cneale/exe/intel/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a(gmx_fatal.o):(.bss+0x530): first   
defined

here
make[3]: *** [grompp] Error 1
make[3]: Leaving directory
`/scratch/cneale/exe/intel/gromacs-4.0.5/src/kernel'
make[2]: *** [all-recursive] Error 1
make[2]: Leaving directory `/scratch/cneale/exe/intel/gromacs-4.0.5/src'
make[1]: *** [all] Error 2
make[1]: Leaving directory `/scratch/cneale/exe/intel/gromacs-4.0.5/src'
make: *** [all-recursive] Error 1


could this 'debug' be the one int calls like this? I at first looked
for a debug() function but didn't find any.

$ grep debug src/gmxlib/*.c
...
main.c:  if (debug) {
main.c:fprintf(debug,"This is simulation %d",cr->ms->sim);
main.c:  fprintf(debug,", local number of nodes %d, local nodeid %d",
main.c:fprintf(debug,"\n\n");
...

Thanks,
Chris

#3

-- original message --

Hello all,

I have built MVAPICH2-1.4rc1 with Infiniband support on a CentOS 5.3 box and
am trying to build gromacs-4.0.5 with MPI support.

make is failing with the following error message:

mpicc -O3 -fomit-frame-pointer -finline-functions -Wall -Wno-unused
-funroll-all-loops -std=gnu99 -o grompp grompp.o
./.libs/libgmxpreprocess_mpi.a -L/usr/lib64 ../mdlib/.libs/libmd_mpi.a
/<...snip..>/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a
../gmxlib/.libs/libgmx_mpi.a -lxml2 -lnsl -lfftw3f -lm -lX11

/usr/local/lib/libmpich.a(ibv_channel_manager.o):(.bss+0x10): multiple
definition of `debug'
<..snip..>/gromacs-4.0.5/src/gmxlib/.libs/libgmx_mpi.a(gmx_fatal.o):(.bss+0x0):
first defined here
collect2: ld returned 1 exit status
make[3]: *** [grompp] Error 1

It looks like 'debug' is conflicting in both MVAPICH2 and in Gromacs.

Any suggestions?

Thanks,

Jerry.
Research Systems Administrator
SBGrid







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

[gmx-users] LINCS warnings with timestep of 5 fs (but not 4 fs)

[Chris Neale]

Hello,

I am recently experimenting with a 5 fs timestep using virtual hydrogens 
and LINCS. I am getting sporadic LINCS warnings, although the system 
appears stable and does not blow up during the 100ns in which I get such 
warnings. If I reduce the timestep to 4 fs from 5 fs (and increase 
nstlist to 5 from 4), then I do not get any such LINCS warnings. It 
occurs to me that in a larger timestep there could be a larger rotation 
and that this warning may be entirely benign. I am not sure about that, 
however, and would welcome any advice here, especially on the way that 
the 30 deg threshold for a warning message is derived.


Thank you,
Chris.


Here is a snippit from my stderr output during a run:

...
Step 1387503, time 6937.51 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000156, max 0.001672 (between atoms 17832 and 17831)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 17832  17831   30.90.0946   0.0947  0.0945

Step 1387503, time 6937.51 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000139, max 0.001582 (between atoms 17832 and 17828)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 17832  17831   31.10.0946   0.0946  0.0945

Step 3601434, time 18007.2 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.35, max 0.000308 (between atoms 18168 and 18165)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 18168  18167   30.10.0945   0.0945  0.0945

Step 4391822, time 21959.1 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000121, max 0.002041 (between atoms 14190 and 14189)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 14190  14189   30.50.0945   0.0947  0.0945

Step 4796357, time 23981.8 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000128, max 0.001443 (between atoms 15532 and 15531)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 15532  15531   31.50.0945   0.0946  0.0945

Step 4796357, time 23981.8 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000111, max 0.001185 (between atoms 15532 and 15531)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 15532  15531   31.80.0945   0.0946  0.0945

Step 4796358, time 23981.8 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000221, max 0.002694 (between atoms 15532 and 15531)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
 15532  15531   34.00.0946   0.0948  0.0945
...




With a .mdp file that looks like this:

integrator  =  sd
nsteps  =  2000
tinit   =  0
dt  =  0.005
comm_mode   =  linear
nstcomm =  4
comm_grps   =  System
nstxout =  2000
nstvout =  2000
nstfout =  2000
nstlog  =  0
nstlist =  4
nstenergy   =  2
nstxtcout   =  2
ns_type =  grid
pbc =  xyz
coulombtype =  PME
rcoulomb=  0.9
fourierspacing  =  0.12
pme_order   =  4
vdwtype =  cut-off
rvdw_switch =  0
rvdw=  1.4
rlist   =  0.9
DispCorr=  no
Pcoupl  =  Berendsen
pcoupltype  =  semiisotropic
compressibility =  4.5e-5 4.5e-5
ref_p   =  1. 1.
tau_p   =  4.04.0
tcoupl  =  Berendsen
tc_grps =  System
tau_t   =  1.0
ref_t   =  310.
ld-seed =  -1
annealing   =  no
gen_vel =  yes
unconstrained-start =  no
gen_temp=  310.
gen_seed=  -1
constraints =  all-bonds
constraint_algorithm=  lincs
lincs-iter  =  1
lincs-order =  6
;;;EOF


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

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

[gmx-users] Huge acceleration needed to reproduce results!

[Chris Neale]

Gromacs uses nm as the unit of distance. Did you account for that? If so, please add some .mdp file 
snippits and any other relevant files so that we can see directly what you are doing.


Chris.

--- original message ---

Hello users,

I am trying to reproduce a calculation that I carried out in DL_POLY.  
It is to calculate the transport diffusion coefficient for CH4 in a  
frozen mesoporous silica.


In DL_POLY I used an external force of 0.1 kJ mol-1 A-1. (0.1 KJ per  
mole per angstrom). This equates to 10 dl_poly internal units which I  
add in this way at each timestep;


Fsum = Fsum + Fex

In Gromacs, I want to apply the same force as I used in DL_POLY so I  
calculated the required acceleration using F=ma. Where I took the mass  
to be the mass of one molecule of methane (16 a.m.u).


The final value for acceleration that I came up with (which  
corresponds to a force of 0.1kj mol-1 A-1 on each molecule) was 0.0625  
nm ps-2.


The first hiccup was when I used this value, the MSD was negative  
(though linear in the negative region of the graph). I assumed that  
this had something to do with the orientation of the unit cell and  
tried applying 0.0   0.0  -0.0625. The plot then looked much better.


The problem is when I calculate the Mean displacement of the CH4  
molecules. (I do this using a slightly altered version of the g_msd  
code). The Mean displacement  from gromacs is very different to that  
which I calculate using DL_POLY,


Gromacs gives 95.0, dl_poly 21347.0.

The MSD however (where I don?t add an acceleration are similar) so the  
problem lies with the force I am adding.


To test that it wasn?t some bug in my code to calculate the Mean  
displacement, I also looked at how the acceleration/force altered the  
MSD in DL_POLY and gromacs.


In DL_POLY, adding an external force of 0.1kj mol-1 A-1 would change  
the MSD of methane by 3 orders of magnitude compared to a run with no  
force added.


My equivalent acceleration of -0.0625 in gromacs, in comparison,  
barely changes the MSD from that of a run with no acceleration added.  
In fact it takes an acceleration of -200 in the z direction to cause  
such a difference in the Ds coefficients between runs with and without  
acceleration.


Does anyone have any idea what is going on here? An acceleration of  
200 ps nm-2 surely is not reasonable is it?. It seems very large  
compared to the example in the manual of 0.1. This would then imply  
that my back of the envelope calculation for relating force and  
acceleration is wrong. Am I missing something? I?m quite sure that the  
force I am adding in DL_POLY is equivalent to 0.1 kJ mol-1 A-1 so why  
are my methane molecules moving so much less in gromacs in response to  
the equivalent acceleration?


Also I noticed that although in the .mdp file I specify:

; Non-equilibrium MD stuff
acc-grps = CH4
accelerate   = 0.0 0.0 -200.00

In the md.log file I get the following output

acc:0   0-154.549   0   0 45.4514

Can someone clarify what this means?

Any advice/comments appreciated



--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

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

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

[gmx-users] OPLS-AA: LJ problem? Atoms overlapping...

[Chris Neale]

Sounds like your topology/forcefield is incorrect or incomplete.

Does your log file from mdrun or the stderr/stdout from grompp/mdrun indicate 
that you don't have any LJ parameters for
some atom types? Try a gmxdump on your .tpr file and inspect the LJ parameters 
by hand. Are there any zero's where you should have numbers?
Perhaps you mixed up some atom types (e.g. with capitilization) or you simply 
didn't add the necessary values.

Chris.

-- original message --

Hello,

I'm using ffoplsaa with a new atom type (Ag) that I've included in the .atp
file, LJ parameters I've specified in the nb file.  My system topology I've
written by hand in an itp file.  The problem is that after any kind of run
the molecule penetrates the metal surface!  This shouldn't be happening.  I
had been using g45a3 quite successfully for the same system and everything
was working fine.  Do I need to do something different/special with oplsaa?
What am I missing?

Advice appreciated,
Chris Rowan

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

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

[gmx-users] Re: OPLS-AA: LJ problem? Atoms overlapping...

[Chris Neale]

what's your minimum atomic distance achieved between Ag and your other molecule (g_mindist)? 


If it's < 0.09 nm then you definitely do have some error in your setup. For 
proof of that, generate a new system that is identical to yours
but this time use only standard OPLS atom-types and see the problem likely 
disapear.

If not, then the penetration is just because it goes between atoms and in this case, while you may not have an "error" per se, 
you should check your parameters and possibly your assumptions about what should be happening.


Also, it's very difficult to offer advice when you describe a problem without 
showing some of the important files directly.
Can you paste the changes to your .itp files and the ffopsaa.itp main file that 
you are using?

Hi Chris,

-- original message --

I've done a gmxdump and there aren't any zeros for any of the LJ values.
No capital letter problems either.

Chris Rowan



Mon, 27 Jul 2009 12:30:54 -0700

Sounds like your topology/forcefield is incorrect or incomplete.

Does your log file from mdrun or the stderr/stdout from grompp/mdrun indicate
that you don't have any LJ parameters for
some atom types? Try a gmxdump on your .tpr file and inspect the LJ parameters
by hand. Are there any zero's where you should have numbers?
Perhaps you mixed up some atom types (e.g. with capitilization) or you simply
didn't add the necessary values.

Chris.

-- original message --

Hello,

I'm using ffoplsaa with a new atom type (Ag) that I've included in the .atp
file, LJ parameters I've specified in the nb file.  My system topology I've
written by hand in an itp file.  The problem is that after any kind of run
the molecule penetrates the metal surface!  This shouldn't be happening.  I
had been using g45a3 quite successfully for the same system and everything
was working fine.  Do I need to do something different/special with oplsaa?
What am I missing?

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

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

[gmx-users] Re: OPLS-AA: LJ problem? Atoms overlapping...

[chris . neale]

What hydrogen did you use? Without substantially more cut and paste on  
your end I am in the dark and can't help very much.


On the side, I'll mention that your sigma, and especially your  
epsilon, looks very low here. I'm no gold expert though, and I am just  
pointing that out even though that is not likely to be the source of  
your problem.


 AG AG 47   107.8700 0.000   A2.80455e-01  1.73304e-01

I believe that I also asked for your ffoplsaa.itp file. Is your  
comb-rule set correctly? Are you sure that you are getting the  
ffoplsaanb.itp inclusion that you intend? For example, I use


...
[ defaults ]
; nbfunccomb-rule   gen-pairs   fudgeLJ fudgeQQ
1   3   yes 0.5 0.5

#include "/home/cneale/gromacs/oplspope.top/ffoplsaanb.itp"
#include "/home/cneale/gromacs/oplspope.top/ffoplsaabon.itp"
...

in my modified ffoplsaa.itp where the directories are hard coded to  
ensure that I get what I expect (rather than some unmodified  
ffoplsaanb.itp from $GMXLIB)


I can't think of anything else. You need to send a lot more data to  
this list. I am currently working under the assumption that you have  
made a mistake somewhere, and thus also assume that one of your  
descriptive statements is either misleading or incorrect -- hence the  
need for seeing the actual files (on list please). If, on the other  
hand, you are sure that you have done everything correctly then I  
can't help.


Good luck,
Chris

-- original message --

Here are some more details:

My cg minimized structure has a distance of 0.3 Angstroms between Ag
and another atom of another molecule.
Changing all my Ag to hydrogens and running a cg still gives a
distance of 0.3 Angstroms.

I could also add that each Ag belongs to its own charge group, so they
should be "felt".
I've tried starting with a different geometry as well with similar results.

Files were modified as follows:

ffoplsaa.atp  I added Silver as:

 AG107.87

ffoplsaanb.itp  I've added:

 AG AG 47   107.8700 0.000   A2.80455e-01  1.73304e-01

Chris Rowan


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

[gmx-users] POSITION-RESTRAIN

[chris . neale]

Sounds to me like you have a 52 atom lipid and have a posre.itp file  
with more than 52 atoms (or referencing atom numbers greater than 52).  
What one usually does is has a somelipid.itp file that, upon some .mdp  
define, includes a someposre.itp file that looks like this:


[ position_restraints ]
; atom  type  fx  fy  fz
 8 1  1000  1000  1000

where for a 52 atom lipid I assume that you are using Berger lipids  
and atom 8 is the Phosphorous.


On the other hand, if you have a .itp file that specifies N-lipids,  
then you would need lots of [ position_restraints ]. Caution: this is  
not the best way to go as you will need to modify many files every  
time you set up a new bilayer.


As an analogy, look at the way that the SOL (tip4p in this case) posre  
is applied in a generic .top out of pdb2gmx:


#include "tip4p.itp"
#ifdef POSRES_SOL
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct   fcxfcyfcz
   11   1000   1000   1000
#endif

where it is a per-molecule atomic index definition such that "atom 1"  
is applied to each solvent molecule.


If you are still unclear, I strongly suspect that there is a gmx  
manual description of all of this, and suggest that you read the  
manual in search of it. If you can't find one then probably you should  
complain back to this list about the absence of a section that should  
definitely be there.


Chris.

-- original message --

From: Morteza Khabiri 
Subject: [gmx-users] POSITION-RESTRAIN
To: gmx-users at gromacs.org
Date: Monday, 27 July, 2009, 7:35 PM

Dear gmxusers.

I want to restrain the lipids in my system which contains protein,lipid
and water. I make the restraint itp file by genpr then I added it in
toplogy file.
After doing grompp to make tpr file I get the following message:

Fatal error:
[ file "posre_entirelipid1.itp", line 56 ]:
 Atom index (53) in position_restraints out of bounds (1-52)
I found the similar error in the mailing list and they suggested probably
the place of restraint itp file which was included in  topology file is
wrong.  However, I tried several positions for restraint itp but I think
it is not the
solution. Do you have any other suggestion about this problem??

Thanks


hi,
Mark is right.the position restrain.itp file should be included in  
your topology file in the proper place otherwise the problem would  
continue. place the .itp file after the section

; Include Position restraint file

#ifdef POSRES

#include "posre.itp"

#endif

best wishes
Shamik



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

Re: [gmx-users] Huge acceleration needed to reproduce results!

[chris . neale]

Lincs
; Do not constrain the start configuration
continuation = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR= no
; Relative tolerance of shake
shake-tol= 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order  = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter   = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle  = 30
; Convert harmonic bonds to morse potentials
morse= no

; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl   =


; Non-equilibrium MD stuff
acc-grps = CH4
accelerate   = 0.0 0.0 -200.0
freezegrps   = SI_O
freezedim= Y Y Y
cos-acceleration = 0
deform   =

and a shortened version of my top file:

[ defaults ]
; nbfunccomb-rule   gen-pairs   fudgeLJ fudgeQQ
1   2   no  1.01.0
;
;
[ atomtypes ]
;   typemasschargeptype   c6c12
 SI 28.08  1.28 A0.000 0.00
 O  15.999-0.64 A   0.27081.538176
 OH 15.999-0.53 A0.30  1.538176
 H   1.008 0.21 A0.000 0.000

;
; Include forcefield parameters
#include "CH4.itp"
;
;
[ moleculetype ]
;   Namenrexcl
MCM 3
[ atoms ]
;   nr  typeresnr   residue atomcgnrcharge  mass
1   SI  1   MCM SI  1   1.2804993   28.086  
2   SI  1   MCM SI  2   1.2804993   
28.086  ..
287 SI  1   MCM SI  287 1.2804993   28.086
288 SI  1   MCM SI  288 1.2804993   28.086
289 O   1   MCM O   289 -0.64024965 15.9994
290 O   1   MCM O   290 -0.64024965 15.9994
.
1070OH  1   MCM OH  1070-0.52612471 15.9994
1071H   1   MCM H   10710.20599988  1.00797
[ bonds ]
;   ai  aj  funct   c0  c1  c2  c3
 286 289 6   0.16 260510
 8   10586   0.16 260510
 ..
[ constraints ]
;   ai  aj  funct   c0  c1  c2  c3
 796 797 1   0.0945
 798 799 1   0.0945
  [ angles ]
;   ai  aj  ak  funct   c0  c1
 365 1   414 1   109.04  416.8167
 365 1   631 1   109.04  416.8167
?..
 537 288 728 1   109.04  416.8167
 592 288 728 1   109.04  416.8167
 225 289 286 1   72.364  180
 66  290 117 1   72.364  180
  
 9   794 239 1   72.364  180
 175 795 228 1   72.364  180
 218 796 797 1   108.5   460.954
 2   798 799 1   108.5   460.954
  [ dihedrals ]
;   ai  aj  ak  al  funct   c0  c1
  709 8   105810595   0   0   0   0
  403 8   105810595   0   0   0   0
  603 8   105810595   0   0   0   0
  962 4   104410455   0   0   0   0
  366 4   104410455   0   0   0   0
  524 4   104410455   0   0   0   0
  1012259 295 150 5   0. 0. 0.

[ system ]
; Name
CH4 in MCM

[ molecules ]
; Compound#mols
MCM1
CH4340


And my .itp file

[ atomtypes ]
;   type  masschargeptype   c6c12
 CH416.043 0.00 A0.37321.24650457
;
[ moleculetype ]
; name  nrexcl
CH42

[ atoms ]
;   nr  typeresnr   residu  atomcgnrcharge  mass
1   CH4  1   CH4 CH4 10.00  16.043









Quoting Chris Neale :


Gromacs uses nm as the unit of distance. Did you account for that? If
so, please add some .mdp file snippits and any other relevant files so
that we can see directly what you are doing.

Chris.

--- original message ---

Hello users,

I am trying to reproduce a calculation that I carried out in DL_POLY.
It is to calculate the transport diffusion coefficient for CH4 in a
frozen mesoporous silica.

In DL_POLY I used an external force of 0.1 kJ mol-1 A-1. (0.1 KJ per
mole per 

Re: [gmx-users] Huge acceleration needed to reproduce results!

[chris . neale]

105810595   0   0   0   0
 962 4   104410455   0   0   0   0
 366 4   104410455   0   0   0   0
 524 4   104410455   0   0   0   0
 1012259 295 150 5   0. 0. 0.

[ system ]
; Name
CH4 in MCM

[ molecules ]
; Compound#mols
MCM1
CH4340


And my .itp file

[ atomtypes ]
;   type  masschargeptype   c6c12
CH416.043 0.00 A0.37321.24650457
;
[ moleculetype ]
; name  nrexcl
CH42

[ atoms ]
;   nr  typeresnr   residu  atomcgnrcharge  mass
1   CH4  1   CH4 CH4 10.00  16.043









Quoting Chris Neale :


Gromacs uses nm as the unit of distance. Did you account for that? If
so, please add some .mdp file snippits and any other relevant files so
that we can see directly what you are doing.

Chris.

--- original message ---

Hello users,

I am trying to reproduce a calculation that I carried out in DL_POLY.
It is to calculate the transport diffusion coefficient for CH4 in a
frozen mesoporous silica.

In DL_POLY I used an external force of 0.1 kJ mol-1 A-1. (0.1 KJ per
mole per angstrom). This equates to 10 dl_poly internal units which I
add in this way at each timestep;

Fsum = Fsum + Fex

In Gromacs, I want to apply the same force as I used in DL_POLY so I
calculated the required acceleration using F=ma. Where I took the mass
to be the mass of one molecule of methane (16 a.m.u).

The final value for acceleration that I came up with (which
corresponds to a force of 0.1kj mol-1 A-1 on each molecule) was 0.0625
nm ps-2.

The first hiccup was when I used this value, the MSD was negative
(though linear in the negative region of the graph). I assumed that
this had something to do with the orientation of the unit cell and
tried applying 0.0   0.0  -0.0625. The plot then looked much better.

The problem is when I calculate the Mean displacement of the CH4
molecules. (I do this using a slightly altered version of the g_msd
code). The Mean displacement  from gromacs is very different to that
which I calculate using DL_POLY,

Gromacs gives 95.0, dl_poly 21347.0.

The MSD however (where I don?t add an acceleration are similar) so the
problem lies with the force I am adding.

To test that it wasn?t some bug in my code to calculate the Mean
displacement, I also looked at how the acceleration/force altered the
MSD in DL_POLY and gromacs.

In DL_POLY, adding an external force of 0.1kj mol-1 A-1 would change
the MSD of methane by 3 orders of magnitude compared to a run with no
force added.

My equivalent acceleration of -0.0625 in gromacs, in comparison,
barely changes the MSD from that of a run with no acceleration added.
In fact it takes an acceleration of -200 in the z direction to cause
such a difference in the Ds coefficients between runs with and without
acceleration.

Does anyone have any idea what is going on here? An acceleration of
200 ps nm-2 surely is not reasonable is it?. It seems very large
compared to the example in the manual of 0.1. This would then imply
that my back of the envelope calculation for relating force and
acceleration is wrong. Am I missing something? I?m quite sure that the
force I am adding in DL_POLY is equivalent to 0.1 kJ mol-1 A-1 so why
are my methane molecules moving so much less in gromacs in response to
the equivalent acceleration?

Also I noticed that although in the .mdp file I specify:

; Non-equilibrium MD stuff
acc-grps = CH4
accelerate   = 0.0 0.0 -200.00

In the md.log file I get the following output

acc:0   0-154.549   0   0
45.4514

Can someone clarify what this means?

Any advice/comments appreciated






--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.







Dr. Jennifer Williams
Institute for Materials and Processes
School of Engineering
University of Edinburgh
Sanderson Building
The King's Buildings
Mayfield Road
Edinburgh, EH9 3JL, United Kingdom
Phone: ++44 (0)131 650 4 861


--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.







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

[gmx-users] No default Ryckaert-Bell. types error for N-terminus acetylated peptide, opls-aa force field

[chris . neale]

The inability to search the archives is certainly leading to some  
redundancy. I use google now, it does a decent job of digging up stuff  
from the gmx mailing list archives -- Even the oldwww search link  
takes me to the new site inside of the iframe so that is no use.


I suggest that the search link on www.gromacs simply link out to  
google for now as the current search doesn't appear to work.


In any event, the answer is here:
http://oldwww.gromacs.org/pipermail/gmx-users/2006-September/023872.html
http://oldwww.gromacs.org/pipermail/gmx-users/2006-September/023875.html

Chris.

-- original messsage --

I am trying to solvate a peptide which has an acetylated N-terminus in the
opls-aa force field. Running pdb2gmx (v. 4.0.2) gives me the following
error:

Error 1
No default Ryckaert-Bell types.

The guilty atoms in the topology are 4 atoms at the N-terminus. This issue
has been raised before as well in the mailing lists.

Is there a way to fix this?

thank you

-Maria.

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

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

[gmx-users] No default Ryckaert-Bell. types error for N-terminus acetylated peptide, opls-aa force field

[chris . neale]

the archive *can* still be searched from the "old" site:  
http://oldwww.gromacs.org/swish-e/search/search2.php


Thanks Justin, this is very valuable information. I guess that this is  
where a new search link on the new website should point for now.


Chris.

-- original message --

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

[Hide Quoted Text]
The inability to search the archives is certainly leading to some
redundancy. I use google now, it does a decent job of digging up stuff
from the gmx mailing list archives -- Even the oldwww search link takes
me to the new site inside of the iframe so that is no use.

I suggest that the search link on www.gromacs simply link out to google
for now as the current search doesn't appear to work.
Good idea.  I've posted this information a few times, but in case anyone's
missed it, the archive *can* still be searched from the "old" site:

http://oldwww.gromacs.org/swish-e/search/search2.php

Same old familiar interface :)

-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/search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] HEME-cysteine gromacs simulation

[chris . neale]

You always need a correct topology. The main issue here is that you  
need to have correct parameters. Where did you get your heme and Fe  
parameters? Were you careful about the Fe state oxidation state? I  
suspect that most people use Amber because of their antechamber  
program, which seems like a brilliant idea even if it may overstretch  
it's own parameterization without letting you know. I have no specific  
advice for you here beyond saying that it is worth spending a month  
figuring out what parameters you should really be using (hint: why  
would you be using ffG43a1 if most others use Amber here? do you know  
something that they don't?). This will end up saving you time in the  
long run.


Note that ffG43a1 is the GROMOS forcefield, not the  GROMACS  
forcefield. It is unfortunate that many programs (Amber, Charmm,  
gromacs) have their own similarly named forcefield, but that does not  
mean that the forcefield must be used with the associated program.  
There are forcefields that are not associated with a program (OPLS)  
and programs that never developed their own force field (NAMD,  
Desmond, Tinker, LAMMPS), so it is perfectly ok for you to use the  
amber forcefield with the gromacs program.


Bottom line: read, read, read.

Chris.

-- original message --

I want to run a simulation with  heme group(cytochrome p450). i have  
made a 3-d structure model for a cytochrome p450 .I use the ffG43a1  
force field,someone told me that i need both suitable parameters for  
your desired force field, and then to construct a correct topology. I  
read lots of papers,but most of them used amber force field.Befor  
runing the gromacs molecular simulation,i define the heme group as a  
new residue in the .rtp files,following the required  
formats(atoms,bonds exclusions,angles,impropers,dihedrals)and then  
define a new FE-S bong.


i run the simulation like this:

pdb2gmx -ignh -ff G43a1 -f CYP2W1_CYSHEME.pdb -o new1.pdb -p new1.top  
-water spce

editconf -bt cubic -f new1.pdb -o new2.pdb -d 0.9
genbox -cp new2.pdb -cs spc216.gro -o new3.pdb -p new1.top
grompp -f em.mdp -c new3.pdb -p new1.top -o em.tpr
genion -s em.tpr -o new5.pdb -nname CL- -nn 1 -g ion.log; Edit top  
file, delette em.tpr

grompp -f em.mdp -c new5.pdb -p new1.top -o em.tpr
mdrun -s em.tpr -o em.trr -c new7.pdb -g em.log -e em.edr
grompp -f pr.mdp -c new7.pdb -p new1.top -o pr.tpr
mdrun -s pr.tpr -o pr.trr -c new9.mdp -e pr.edr -g pr.log

 1 i run pdb2gmx and get the new1.pdb ,when checking the pdb with  
accelrys ds visualizer or sybyl,the atom FE was recognized as F  
(Fluorine)and the pdb miss the bonds between FE-S &FE-N,i add or  
delete a blank space befor or after the FE in the pdb file,it was read  
as iron(FE) again but still missing;when checking the .top file,there  
was a FE (both atom name and type)and formed the FE-S & FE-N bonds.


if i ignore the above things,go on the command ,when runing the em.mdp  
there is a tip that Steepest Descents failed converged to Fmax < 1000  
in 5001 steps.when run the pr.tpr,a error emerged as segment fatal  
error¡£


2 someone told me that i need both suitable parameters for your  
desired force field, and then to construct a correct topology. I read  
lots of papers,but most of them used amber force field.Does anyone  
known a paper with a gromacs force field for HEME group?


please help me how handle those problems,how can i go on the simulation?

the HEME-cysteine parameter in the rtp files,define themas a new  
residue.(is there any problems?)


[ HEME ]
 [ atoms ]
N N-0.28000 0
H H 0.28000 0
   CA   CH1 0.0 1
   CB   CH2-0.1 2
   SG S-0.4 2
C C 0.38000 3
O O-0.38000 3
   FEFE 0.84700 4
   NANR-0.37000 4
   NBNR-0.42300 4
   NCNR-0.50400 4

 [ bonds ]
N Hgb_2
NCAgb_20
   CA Cgb_26
C Ogb_4
C+Ngb_9
   CACBgb_26
   CBSGgb_30
   SGFEgb_48
   FENAgb_34
   FENBgb_34
   FENCgb_34
   FENDgb_34

.
 [ angles ]
;  aiajak   gromos type
   -C N H ga_31
H NCA ga_17
   -C NCA ga_30
NCA C ga_12
   CA C+N ga_18
   CA C O ga_29
 ...
 [ impropers ]
;  aiajakal   gromos type
N-CCA H gi_1
CCA+N O gi_1
   CA N CCB gi_2
   FESGNAND gi_3
   FESGNANB gi_3
   FESGNBNC gi_3
   FESGNCND gi_3
  .
 [ dihedrals ]
;  aiajakal   gromos type
  .


ûÓйã¸æµÄÖÕÉíÃâ·ÑÓÊÏä,www.yeah.net




___
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/sear

[gmx-users] SIMULATION OF GOLD NANOPARTICLES

[chris . neale]

No information on what you requested, but I do have some tips on how  
to get help from a mailing list. Please understand that I'm sending  
this only to try to help you.


Don't use all caps. Don't double post, we saw the first one. Show us  
that you have done some work and, if possible, approach the list with  
a question that is more well defined. And just in case, don't send a  
personal mail to somebody on this list asking for advice since the  
list is the place to do that.


Chris.

-- original message --

From:  T? 
Subject: [gmx-users] SIMULATION OF GOLD NANOPARTICLES
To: 
Message-ID: <4c581b47098d4bd489c3a2e734841...@laser7>
Content-Type: text/plain; charset="iso-8859-7"


Hi GROMACS users

Does anyone have any information on the simulation of gold  
nanoparticles with GROMACS?


Thank you in advance

George

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

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

[gmx-users] alternative means of calculating PMF

[chris . neale]

Hi Stefan,

If you know the bound state, you can use the double decoupling method  
to decouple the bound state and then recouple into bulk solvent.


Note that both decoupling and umbrella sampling methods will require  
that you make the appropriate corrections to standard state values and  
that this correction may be what you are missing with your umbrella  
sampling results.


I outline the way to get US standard state free energies, and this is  
a very quick description so please refer to peer reviewed papers for  
something that is more complete.


1. Boltzmann integrate your PMF (E_i values) over the bound state.

G_bound=-1/Bln*sum over i [exp(-BE_i)*delta_i]

2. Define your unbound energy (E_ub) from your PMF and Boltzmann  
integrate it over a distance, dist, such that the volume per monomer  
will be 1.66 nm^3 (volume per molecule at 1 M concentration).


G_unbound=-1/Bln[exp(-BE_ub)*dist]

3. dG_binding = G_bound - G_unbound

Where I have here neglected any corrections for constraints (i.e. I  
assumed that you reaction coordinate is a Center of Mass relation or  
something similar).


Also, note that the form of equation 2 assumes that you have already  
extracted the increasing volume portion of your PMF if you simply use  
the distance between centers of mass. Therefore if you have not done  
this then your equation 2 will be of a form like equation 1, except  
that you integrate over an unbound region for a total of 1.66 nm^3.


Finally, make sure that you take both monomers into account in your  
final calculations correctly so that your final concentration of  
unbound monomer is really 1 M, not 0.5 M or 2 M.


And if your monomers are so big that the bound state occupies a volume  
greater than 2*1.66 nm^3 just set your standard state to 0.1 M or  
something and then make a final correction back to 1 M later.


Good luck,
Chris.


--- original message ---

Hello.
I have finished calculating the PMF of dimerization for my system using
umbrella sampling and WHAM. But came up with some free energy values that
are not compatible to the ones found in other works. I have alread had
several instrutions and tips from Justin in past email here in gms users'
list, and since now I believe there are no more variables to change in my
calculations, I believe another method (implemented in gromacs) would help
me compare results and justify my findings.
I read something in the users'list about instead of using umbrella sampling
I would use constraints in order to separate my system's components. But I
was not able to find anything more specific about this type of calculation
on the manual.
Some help on the matter would be of great help
Thank you

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

[gmx-users] pull code with defined negative relative displacements

[chris . neale]

Hello, I am re-running some of our gromacs 3 simulations using gromacs
4, and as far as I can tell the gromacs 4 pull code, while very nicely
enhanced from gromacs 3, has also lost some functionality.

I am calculating the PMF of a peptide across a bilayer and, to  
simplify the issue, what I can't figure out how to do with gromacs 4  
is to pull group1 to a position 1 nm more negative along the z-axis  
than group 0 (the bilayer in this case) for some replicas and to 1 nm  
more positive along the z-axis than group 0 in others. I'll focus on  
the negative displacement here as it is the one that is giving me  
problems.


This used to be possible with the following gromacs 3 pull.ppa options:

runtype   = umbrella
reftype   = com
pulldim   = N N Y
reference_group = group0
group_1   = group1
K1= 500
Pos1  = 0 0 -1.0

Sure, I could start group1 in the correct negative-z displacement from
group0 and use pull_init1=+1.0, but this will not work when, for
instance, I want to restrain it to -0.1 nm (using pull_init1=+0.1),
where the sampling will infrequently jump back and forth about z=0.

Just to be sure, I tried for gromacs 4 are the following pull code  
.mdp options:


pull = umbrella
pull_geometry= distance
pull_dim = N N Y
pull_start   = no
pull_ngroups = 1
pull_group0  = group0
pull_group1  = group1
pull_init1   = -1.0
pull_rate1   = 0
pull_k1  = 500.0

where mdrun complains:
"Pull reference distance for group 1 is negative (-1.00)"

and it is pretty obvious why this doesn't work since it is asking for
a negative displacement. Nevertheless, I tried it and pull_init1
appears to get set to zero.

I also attempted the following.

pull = umbrella
pull_geometry= distance
pull_dim = N N Y
pull_start   = no
pull_ngroups = 1
pull_group0  = group0
pull_group1  = group1
pull_init1   = 1.0
pull_rate1   = 0
pull_k1  = 500.0
pull_vec1= 0 0 -1

where I would then use
pull_init1 = 1.0
pull_vec1  = 0 0 1
for the positive side of the bilayer.

However, when I look at the forces, I am getting the negative of what  
I should get when pull_vec1 = 0 0 -1.


coord.xvg:
610.4.9343  -1.02019
660.00014.91454 -0.949747

force.xvg:
610.-10.0932
660.000125.1265

Although I am getting exactly what I should get when pull_vec1 = 0 0 1  
(for intended positive displacements).


coord.xvg:
660.00014.9014  1.16304

force.xvg:
660.0001-81.5201

Any ideas are greatly appreciated. I can probably mod the code for my
needs, but a standard gromacs binary is always preferable.

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/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

[gmx-users] mdrun

[chris . neale]

Run some more energy minimization (EM). Make sure that you have  
constraints=none during EM (or I believe that you could also increase  
the LINCS/SHAKE iterations to the appropriate number). Also, while  
some people disagree in principle with me here, I also find in  
practice that a preliminary EM with define=-DFLEXIBLE (or otherwise  
depending on your water model) will help in some cases. If your EM is  
always stopping after 15 steps with max force=inf, then you need to go  
back to some of your previous setup steps and do them again while  
being a little more fastidious about the details (e.g. is your box  
size appropriate for your coordinates).


that error is telling you about what happened *at* 0.005 ps (t is not  
equal to dt).


Next time you post, please provide more information about how you set  
up and EM'd your system. Also, a better title gets better assistance.


Chris.

-- original message --

Hi

I did command (mdrun -s -o -c -e -g) but following error was came up where
as my dt is 0.001 ps in mdp file:

t = 0.005 ps: Water molecule starting at atom 20937 can not be settled.
Check for bad contacts and/or reduce the timestep.
Wrote pdb files with previous and current coordinates.

I reduced time step but same error was came up again.I enlarged box size but
same error was came up again.

please guide me.


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

[gmx-users] About rotational restraint

[chris . neale]

I haven't done it myself, but both of these methods appear as if they  
could be implemented out of the box in gromacs.


On the use of orientational restraints and symmetry corrections in  
alchemical free energy calculations. Mobley DL, Chodera JD, Dill KA.

http://www.ncbi.nlm.nih.gov/pubmed/16965052

Absolute Binding Free Energy Calculations Using Molecular Dynamics  
Simulations with Restraining Potentials. Jiyao Wang, Yuqing Deng, and  
Benoît Roux

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

Chris

-- original message --

I'm applying COM pull on a two-domain protein by fixing the distance between
the centers of mass of each domain with an umbrella potential. However, I
find that the domains can rotate so that the inner surfaces of each domain
are not parallel any more. How can I apply a rotational restraint so that
each domain doesn't rotate during simulation? 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/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

[gmx-users] pull code with defined negative relative displacements

[chris . neale]

Hi Justin, (comment to Berk below)

thanks for pointing out the g_wham problem. I am personally ok here as  
I use a non-gromacs wham program.


Hi Berk,

unless I misunderstood, your suggestion does not yield the gromacs 3  
behaviour that I am trying to reproduce with gromacs 4.


If I simply misunderstood your suggestion, can you please be a little  
more explicit? What I want to do is to harmonically restrain the COM  
of group 1 to be X nm more negative along Z than the COM of group 0.  
It is important that this still works when the value for X is smaller  
than the standard deviation of the sampled values such that the  
distribution does not become bimodal.


Here is the test using pull_geometry=direction, and I note that this  
is exactly the same behaviour that I observed and posted in my last  
email while using pull_geometry=distance in combination with pull_vec.


pull = umbrella
pull_geometry= direction
pull_dim = N N Y
pull_start   = no
pull_nstxout = 10
pull_nstfout = 10
pull_ngroups = 1
pull_group0  = POPC
pull_pbcatom0= 0
pull_group1  = Protein
pull_pbcatom1= 0
pull_init1   = 4.75
pull_rate1   = 0
pull_k1  = 500.0
pull_vec1   = 0 0 -1


coord.xvg:
1669.4000   4.97979 -4.68988
force.xvg:
1669.4000   30.0622
Where the force should be -30.062 (assuming that I understand  
everything correctly).


coord.xvg:
2.9600  5.16151 -5.14321
force.xvg
2.9600  -196.603
Where the force should be +196.603 (assuming that I understand  
everything correctly).


pull = umbrella
pull_geometry= direction
pull_dim = N N Y
pull_start   = no
pull_nstxout = 10
pull_nstfout = 10
pull_ngroups = 1
pull_group0  = POPC
pull_pbcatom0= 0
pull_group1  = Protein
pull_pbcatom1= 0
pull_init1   = 4.75
pull_rate1   = 0
pull_k1  = 500.0
pull_vec1   = 0 0 1

coord.xvg:
0.5200  5.17421 5.00621
force.xvg:
0.5200  -128.107
Where this is what I intend to occur (i.e. working as I desire it to).

coord.xvg:
1649.4000   4.99895 4.69239
force.xvg:
1649.4000   28.807
Where this is what I intend to occur (i.e. working as I desire it to).




Thank you,
Chris.

--- original message --

Hi,

You should use pull_geometry=direction.
distances don't get negative.

Berk


Date: Sat, 14 Nov 2009 09:21:39 -0500
From: chris.neale at utoronto.ca
To: gmx-users at gromacs.org
Subject: [gmx-users] pull code with defined negative relative displacements

Hello, I am re-running some of our gromacs 3 simulations using gromacs
4, and as far as I can tell the gromacs 4 pull code, while very nicely
enhanced from gromacs 3, has also lost some functionality.

I am calculating the PMF of a peptide across a bilayer and, to   
simplify the issue, what I can't figure out how to do with gromacs 4  
 is to pull group1 to a position 1 nm more negative along the z-axis  
 than group 0 (the bilayer in this case) for some replicas and to 1  
nm  more positive along the z-axis than group 0 in others. I'll  
focus on  the negative displacement here as it is the one that is  
giving me  problems.


This used to be possible with the following gromacs 3 pull.ppa options:

runtype   = umbrella
reftype   = com
pulldim   = N N Y
reference_group = group0
group_1   = group1
K1= 500
Pos1  = 0 0 -1.0

Sure, I could start group1 in the correct negative-z displacement from
group0 and use pull_init1=+1.0, but this will not work when, for
instance, I want to restrain it to -0.1 nm (using pull_init1=+0.1),
where the sampling will infrequently jump back and forth about z=0.

Just to be sure, I tried for gromacs 4 are the following pull code   
.mdp options:


pull = umbrella
pull_geometry= distance
pull_dim = N N Y
pull_start   = no
pull_ngroups = 1
pull_group0  = group0
pull_group1  = group1
pull_init1   = -1.0
pull_rate1   = 0
pull_k1  = 500.0

where mdrun complains:
"Pull reference distance for group 1 is negative (-1.00)"

and it is pretty obvious why this doesn't work since it is asking for
a negative displacement. Nevertheless, I tried it and pull_init1
appears to get set to zero.

I also attempted the following.

pull = umbrella
pull_geometry= distance
pull_dim = N N Y
pull_start   = no
pull_ngroups = 1
pull_group0  = group0
pull_group1  = group1
pull_init1   = 1.0
pull_rate1   = 0
pull_k1  = 500.0
pull_vec1= 0 0 -1

where I would then use
pull_init1 = 1.0
pull_vec1  = 0 0 1
for t

[gmx-users] pull code with defined negative relative displacements

[chris . neale]

Hi Berk, I am aware that pull_init=4.75 should pull to 4.75, not 1.0.  
I was mixing results from my simulations as I was discussing, sorry.


Please allow me to restate the issue: If I am indeed pulling to 4.75  
below group 0, then why is the force positive when group 1 has a  
displacement of -4.68988 (see below):


coord.xvg:
1669.4000   4.97979 -4.68988
force.xvg:
1669.4000   30.0622

I think that pull_geometry=distance and pull_geometry=direction are  
not going to work here.


I have moved to a test system that is the a box of water in which I  
pull the z of one water to a relative displacement with respect to one  
other water.


From this, I have found that pull_geometry=position while passing in  
a negative pull_init1 and a pull_vec1 of 0 0 0, appears to be the only  
way to get this done.


Thank you for all of your time,
Chris.


Hi,

pull_init should be 1, not 4.75, if you want group1 to be 1 nm below group0.
Alternatively you can have vec = 0 0 1 and pull_init = -1.

Berk


Date: Sat, 14 Nov 2009 14:54:54 -0500
From: chris.ne...@utoronto.ca
To: gmx-users@gromacs.org
Subject: [gmx-users] pull code with defined negative relative displacements

Hi Justin, (comment to Berk below)

thanks for pointing out the g_wham problem. I am personally ok here as
I use a non-gromacs wham program.

Hi Berk,

unless I misunderstood, your suggestion does not yield the gromacs 3
behaviour that I am trying to reproduce with gromacs 4.

If I simply misunderstood your suggestion, can you please be a little
more explicit? What I want to do is to harmonically restrain the COM
of group 1 to be X nm more negative along Z than the COM of group 0.
It is important that this still works when the value for X is smaller
than the standard deviation of the sampled values such that the
distribution does not become bimodal.

Here is the test using pull_geometry=direction, and I note that this
is exactly the same behaviour that I observed and posted in my last
email while using pull_geometry=distance in combination with pull_vec.

pull = umbrella
pull_geometry= direction
pull_dim = N N Y
pull_start   = no
pull_nstxout = 10
pull_nstfout = 10
pull_ngroups = 1
pull_group0  = POPC
pull_pbcatom0= 0
pull_group1  = Protein
pull_pbcatom1= 0
pull_init1   = 4.75
pull_rate1   = 0
pull_k1  = 500.0
pull_vec1   = 0 0 -1


coord.xvg:
1669.4000   4.97979 -4.68988
force.xvg:
1669.4000   30.0622
Where the force should be -30.062 (assuming that I understand
everything correctly).

coord.xvg:
2.9600  5.16151 -5.14321
force.xvg
2.9600  -196.603
Where the force should be +196.603 (assuming that I understand
everything correctly).

pull = umbrella
pull_geometry= direction
pull_dim = N N Y
pull_start   = no
pull_nstxout = 10
pull_nstfout = 10
pull_ngroups = 1
pull_group0  = POPC
pull_pbcatom0= 0
pull_group1  = Protein
pull_pbcatom1= 0
pull_init1   = 4.75
pull_rate1   = 0
pull_k1  = 500.0
pull_vec1   = 0 0 1

coord.xvg:
0.5200  5.17421 5.00621
force.xvg:
0.5200  -128.107
Where this is what I intend to occur (i.e. working as I desire it to).

coord.xvg:
1649.4000   4.99895 4.69239
force.xvg:
1649.4000   28.807
Where this is what I intend to occur (i.e. working as I desire it to).




Thank you,
Chris.

--- original message --

Hi,

You should use pull_geometry=direction.
distances don't get negative.

Berk

> Date: Sat, 14 Nov 2009 09:21:39 -0500
> From: chris.neale at utoronto.ca
> To: gmx-users at gromacs.org
> Subject: [gmx-users] pull code with defined negative relative   
displacements

>
> Hello, I am re-running some of our gromacs 3 simulations using gromacs
> 4, and as far as I can tell the gromacs 4 pull code, while very nicely
> enhanced from gromacs 3, has also lost some functionality.
>
> I am calculating the PMF of a peptide across a bilayer and, to
> simplify the issue, what I can't figure out how to do with gromacs 4
>  is to pull group1 to a position 1 nm more negative along the z-axis
>  than group 0 (the bilayer in this case) for some replicas and to 1
> nm  more positive along the z-axis than group 0 in others. I'll
> focus on  the negative displacement here as it is the one that is
> giving me  problems.
>
> This used to be possible with the following gromacs 3 pull.ppa options:
>
> runtype   = umbrella
> reftype   = com
> pulldim   = N N Y
> reference_group = group0
> group_1   = group1
> K1= 500
> Pos1  = 0 0 -1.0
>
> Sure, I could start group1 in the correct negative-z displacement from
> group0 and use pull_init1=+1.0, but this will not work when, for
> instance, I want to restrain it to -0.1 nm (using pull_ini

[gmx-users] pull code with defined negative relative displacements

[chris . neale]

I only now noticed Justin mail on g_wham.
You can probably also use pull_geometry=distance and pull_init=1,
if you starting structure has group1 close to 1 nm below group 0.


Agreed, although this will not work when the force constant is not  
strong enough to inhibit any sampling >0 -- wherein the distribution  
about 0 would become bimodal and this is something that I can not allow.


As far as I can tell, pull_geometry=position, pull_init1<0, and  
pull_vec1=0 0 0 is the solution, as per my previous message.


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/search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] pull code with defined negative relative displacements

[chris . neale]

Hi Berk,

I do not mean the box 0, and I am aware about the pull distance  
needing to be less than half of the smallest box length. Please allow  
me describe a bit more rigorously what I need to avoid. I will use an  
example that is totally fictitious, but is designed to emphasize the  
bimodal sampling issue that is relevant to my actual study.


1. Imagine that one wants to calculate the PMF for a water molecule  
along the normal to an asymmetric bilayer in which one leaflet is POPC  
and the other is lipidA.


2. The free energy for that water residing 0.1 nm away from the  
bilayer center in the POPC leaflet (call it -0.1 nm from the bilayer  
center) needs to be calculated separately from the free energy for  
that water residing 0.1 nm away from the bilayer center in the lipidA  
leaflet (call it +0.1 nm from the bilayer center).


3. Let the force constant be 500 kJ/mol/nm^2 and let the user apply  
pull_geometry=distance and pull_init1=0.1 while starting the replica  
at 0.1 nm from the bilayer center in two simulations where one takes  
(a) a starting water position that is 0.1 nm toward the POPC  
headgroups and the other takes (b) a starting position of the water to  
be 0.1 nm toward the lipidA headgroups.


4. I posit that the probability distributions of the water about the  
bilayer normal from simulations (a) and (b) in part 3 above will  
converge to the same distribution and that this distribution will be  
bimodal with one peak at -0.1 nm and one peak at +0.1 nm. This is  
because the simulation that started at -0.1 (0.1 nm closer to the POPC  
headgroups) will, even under the biasing force, infrequently cross the  
bilayer center of mass and then be closer to the lipidA heagroups. At  
this point. the direction of the applied force will the change to draw  
the water toward the other possible location at +0.1 nm.


5a. For displacements >> 0.1 nm, this will not be a problem.
5b. For larger force constants, this will not be a problem at 0.1 nm,  
but will still be a problem for some closer centers of restraint.


6. pull_geometry=distance is entirely incompatible with this approach.

7. pull_geometry=position does appear to work.


Thanks for your patience,
Chris.

-- original message --

Hi,

Sorry, but I have no clue what you mean with sampling >0
and how you would end up with a bimodal distribution.

You don't mean the box 0, do you?
That is irrelevant.
What can be a problem is that you pull distance should not be more  
than half the box length.


Berk


Date: Sun, 15 Nov 2009 12:33:30 -0500
From: chris.neale at utoronto.ca
To: gmx-users at gromacs.org
Subject: [gmx-users] pull code with defined negative relative displacements

> I only now noticed Justin mail on g_wham.
> You can probably also use pull_geometry=distance and pull_init=1,
> if you starting structure has group1 close to 1 nm below group 0.

Agreed, although this will not work when the force constant is not   
strong enough to inhibit any sampling >0 -- wherein the distribution  
 about 0 would become bimodal and this is something that I can not  
allow.


As far as I can tell, pull_geometry=position, pull_init1<0, and   
pull_vec1=0 0 0 is the solution, as per my previous message.


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/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

[gmx-users] PULLING

[chris . neale]

First let me mention that I only scanned the mammoth manuscript  
snippet. Nevertheless, I'll try to address your questions.


1. You are physically able to use spc in place of tip3p. Whether or  
not that is a good idea is up to you to decide. Read the literature  
including the spc, tip3p, and other relevant forcefield papers in  
addition to any comparative studies.


2. No obvious question posed.

3. use the pull code. I don't know exactly what you want to do here...  
perhaps harmonically restrain the distance between the two mentioned  
groups to some specified value? No need to rename the water, just make  
an index group.


4. Who knows. Probably they did some preliminary simulations and  
discovered that this force constant did what they wanted. Your results  
can be entirely de-biased and so any force constant is technically  
rigorous, but some are better than others. You want it to be strong  
enough to hold your sampling near the defined displacements, but also  
weak enough that there is overlap between adjoining replicas. Probably  
time to read about the methodology here.


Chris.


Hi all.

I have to reproduce an experiment from the article "Identification of the
nucleophilic factors and the productive complex for the editing reaction by
leucyl-tRNA synthetase" (by Yohsuke Hagiwara a,b, Osamu Nureki c, Masaru
Tateno a,b,*). This is one of the steps of education I have to complete. All
actions described in this article were made with Amber9. I'm a newer in
modelling and I have only gromacs tools (v4.0.4), so I have downloaded amber
force field port-files (E.J. Sorin, S. Park). I want to know whether my mdp
files are reflexing the data I need. So this is a chapter from article:



"...1. MD simulations were performed using the AMBER 9 program. The parm99
force field was applied to the atoms of LeuRS and tRNALeu, and the TIP3P
model was used for the solvent water molecules. Electrostatic interactions
were calculated by the particle-mesh Ewald (PME) method, using 1.0 as the
dielectric constant. A cut-off of 12 ? was used to calculate the direct
space sum for PME; the electrostatic interactions beyond 12 ? were
calculated in reciprocal space by the fast Fourier transform method.

Thus, all electrostatic interactions between atoms were calculated. The
SHAKE algorithm was used to restrain the bond lengths involving hydrogen
atoms. The time step for integration was set to 1 fs. The temperature and
pres-sure were set using the Berendsen algorithm.

The initial coordinates used for the present MD simulations were from the
modelled structure previously constructed for LeuRS complexed with
valyl-tRNALeu. For the present system, we immersed the complex in a solvent
box comprising 49 587 water molecules, and used the periodic boundary
condition where the size of the unit box was 103.0 _ 138.3 _ 117.1 ?3. The
total number of atoms in the system was 165 739. In the initial (relaxation)
phase of the MD simulation, the water molecules were relaxed at 300 K for 10
ps, while the atoms of the protein and tRNALeu were positionally constrained
by a harmonic function using a force constant of 500 kcal/mol ?2. The force
constant was then reduced to 250, 125, 50, 25, 10, and 5 kcal/mol ?2 in six
MD simulations. The time consumed by each simulation was 2 ps...(then free
MD)...

2. This equilibrated system was used for further structural modeling to
investigate the hydrated structure relevant to the editing reaction, in
which 1 ns MD simulations were performed in the presence of constraints to
effectively explore the states. First, with respect to the atomic distance
between the carbonyl carbon of the substrate and the oxygen atom of the
identified nucleophilic water, a constrained MD simulation was performed to
investigate the

mechanisms of the approach of the water molecule. This simulation was
started by using a harmonic potential as the distance constraint, where the
force constant was set to 5 kcal/mol ?2. When the atomic distance was not
reduced any further, despite the presence of the harmonic potential in the
MD simulation (at about 470 ps), the force constant was increased up to 200
kcal/mol ?2; this was exploited to mimic the first phase of a nucleophilic
attack in the MD simulations, which enabled us to explore larger
conformational spaces than with the first principles MD simulations. A
similar protocol was applied in other constrained MD simulations for
efficient explorations of conformational spaces. To obtain a PMF with
respect to the rotation of the dihedral angle, C40-C30-O30-HO30 , we
employed the umbrella sampling technique, using 14 windows in the range of
50-180_. Umbrella sampling is a theoretical technique to efficiently search
for not only energetically favourable but also energetically-unfavourable
conformations in a phase space of a system, combined with a bias potential
(e.g. a harmonic potential) to overcome energy barriers. The force constants
of the  potential in those windows were set