[gmx-users] Strategy for chaining jobs using different mdp parameters (long)

Stéphane Téletchéa steletch at jouy.inra.fr
Mon May 21 18:22:20 CEST 2007


Dear colleagues,

I've done many chained simulations up to now but i'm encoutering 
problems and i'm skeptical about some choices i made before. In order to 
clear this up, i'd like you to share your experience, so that'll answer 
my needs and be a good template for the upcoming wiki.

I want (and need) to chain jobs with different conditions for 3 reasons:
1 - cluster needs (need to split jobs), tpbconv could be sufficient, but 
i'm changing parameters
2 - ensemble changing (switching from NVT to NPT for instance)
3 - dynamics fine-control (posre on specific atoms/groups)

Please see below and the attached files for run parameters.

a) NVT then NPT.

I've tried using NVT (pr1) for starting my simulation and then NPT 
(pr2). The major differences between pr1 (NVT) and pr2 (NVT) are:

pr1:
< Pcoupl                 = no
< gen_vel                = yes
---
pr2:
 > Pcoupl                 = berendsen
 > gen_vel                = no

The runs are chained like this:
(ini directory contains the top and gro file from previously minimised 
structure, mdp directory contains the reference mdp)
pr1 step:
grompp \
     -f ../mdp/pr1.mdp \
     -p ../ini/molecule_ini.top \
     -pp molecule_pr1.top \
     -c ../ini/molecule_ini.gro \
     -o molecule_pr1.tpr \
     -po molecule_pr1.mdp

mdrun \
    -s molecule_pr1.tpr \
    -o molecule_pr1.trr \
    -c molecule_pr1.gro \
    -g molecule_pr1.log \
    -e molecule_pr1.edr \
    -nice 0

pr2 step:
grompp \
     -f ../mdp/pr2.mdp \
     -p ../ini/molecule_ini.top \
     -pp molecule_pr2.top \
     -c ../ini/molecule_ini.gro \
     -o molecule_pr2.tpr \
     -t ../pr1/molecule_pr1.trr \
     -po molecule_pr2.mdp

(note there is no -e ../pr1/molecule_pr1.edr since we're switching from 
NVT to NPT, as pointed out by Mark Abraham on the list 
http://www.gromacs.org/pipermail/gmx-users/2007-May/027193.html)

mdrun \
    -s molecule_pr2.tpr \
    -o molecule_pr2.trr \
    -c molecule_pr2.gro \
    -g molecule_pr2.log \
    -e molecule_pr2.edr \
    -nice 0

As you may have noticed, i'm doing positionnal restraints (pr) runs at 
theses steps.
Do i need to specify "unconstrained_start    = yes"?
I've not used it for the moment...

Is there any specific precaution to take while doing this?

In particular, i'm seing for the NPT (pr2) start:
#####
Step 1  Warning: pressure scaling more than 1%, mu: -1.59647e+20 
-1.59647e+20 -1.59647e+20
Correcting invalid box:
old box (3x3):
    old box[    0]={-9.82630e+20,  0.00000e+00, -0.00000e+00}
    old box[    1]={ 0.00000e+00, -8.90833e+20, -0.00000e+00}
    old box[    2]={ 0.00000e+00,  0.00000e+00, -1.12456e+21}
new box (3x3):
    new box[    0]={-9.82630e+20,  0.00000e+00, -0.00000e+00}
    new box[    1]={ 0.00000e+00, -8.90833e+20, -0.00000e+00}
    new box[    2]={ 0.00000e+00,  8.90833e+20, -1.12456e+21}
#####

Of course the system exploses (but this is another story, in the ml, 
David van der Spoel stated this could be due to a system setup problem, 
i'll open a new thread on this specific problem).

I'm not completely sure if the problem arises from my specific system or 
if i'm missing something while switching from NVT to NPT (see below for 
reference files).

b) Reference files while continuing the run (not using tbpconv)

Even if i use the same NPT conditions, i still see some parameters 
problem, for instance:

#####
Warning: atom names in ../pr1/test_min_deflated25_sol_ions.top and 
../pr1/OR1G1_min_formd.gro don't match (OW - HW2)
Warning: atom names in ../pr1/test_min_deflated25_sol_ions.top and 
../pr1/OR1G1_min_formd.gro don't match (HW1 - OW)
Warning: atom names in ../pr1/test_min_deflated25_sol_ions.top and 
../pr1/OR1G1_min_formd.gro don't match (HW2 - HW1)
Warning: atom names in ../pr1/test_min_deflated25_sol_ions.top and 
../pr1/OR1G1_min_formd.gro don't match (OW - HW2)
#####

I've understood this comes from the usage of -shuffle *and* -sort option 
(and actually, to get the proper ordering, there is only g_desort, as 
mentionned in 
http://www.gromacs.org/pipermail/gmx-users/2007-January/025584.html), 
but i'm suspicious about the "-renum" and "-renumbs" options in grompp 
(on by default), if it alters the resulting -processed- topology, or not.

To be clear, i'm not completely sure if i need to do:

starting conformation = ../ini/molecule.top and ../ini/molecule.gro

directories :
mkdir {ini,pr1,pr2}

pr1 step:
grompp \
     -f ../mdp/pr1.mdp \
     -p ../ini/molecule_ini.top \
     -pp molecule_pr1.top \
     -c ../ini/molecule_ini.gro \
     -o molecule_pr1.tpr \
     -po molecule_pr1.mdp

mdrun \
    -s molecule_pr1.tpr \
    -o molecule_pr1.trr \
    -c molecule_pr1.gro \
    -g molecule_pr1.log \
    -e molecule_pr1.edr \
    -nice 0

pr2 step:
grompp \
     -f ../mdp/pr2.mdp \
     -p ../ini/molecule_ini.top \
     -pp molecule_pr2.top \
     -c ../ini/molecule_ini.gro \
     -o molecule_pr2.tpr \
     -t ../pr1/molecule_pr1.trr \
     -po molecule_pr2.mdp

mdrun \
    -s molecule_pr2.tpr \
    -o molecule_pr2.trr \
    -c molecule_pr2.gro \
    -g molecule_pr2.log \
    -e molecule_pr2.edr \
    -nice 0

Please note here, i'm always referring to ../ini/molecule_ini.{top,gro} 
for grompp.

In order to chain the jobs, do i need to reference from the processed 
top and last gro from the previous step. In other words, do i need to use :

pr2 step:
grompp \
     -p ../pr1/molecule_pr1.top \
     -c ../pr1/molecule_pr1.gro \

instead of (as above):
pr2 step:
grompp \
     -p ../ini/molecule_ini.top \
     -c ../ini/molecule_ini.gro \

I presume i need to use the previous step top and gro parameters, but 
using g_desort (not done yet).
Again in this case, is g_desort sufficient ?
Are -renum and -rmbdunbs only used for the run (but not permanently 
affecting the top and gro files)?

I'm trying to figure out the various errors i'm seing (which are 
probably not all related to the chaining setup by itself) and since this 
usage seems trivial (at least as written in papers) and i've not found 
specific answers on these points during my various researches, i hope to 
find out the correct answers here, in order to list them extensively on 
the wiki.

Thanks a lot for your patience,

Stéphane

Joining the two mdp parameters in case the chaining setup needs 
additionnal parameter:
pr1 is for the NVT simulation,
pr2 is for the NPT simulation,
i'm not using the same pr for the different pr steps, but this should 
not be a problem for the "generic" way of chaining simulations steps.

-- 
Stéphane Téletchéa, PhD.                  http://www.steletch.org
Unité Mathématique Informatique et Génome http://migale.jouy.inra.fr/mig
INRA, Domaine de Vilvert                  Tél : (33) 134 652 891
78352 Jouy-en-Josas cedex, France         Fax : (33) 134 652 901
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: pr1.mdp
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20070521/d5f72029/attachment.ksh>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: pr2.mdp
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20070521/d5f72029/attachment-0001.ksh>


More information about the gromacs.org_gmx-users mailing list