[gmx-users] mdrun -rerun: bonded interactions of protein

NG HUI WEN HuiWen.Ng at nottingham.edu.my
Mon Nov 22 11:19:41 CET 2010


Dear Mark and Justin,

 

I have a potentially silly question here. I have been trying to work out
the figures from the various options of g_energy. You mentioned before
that the bonded terms are not decomposable/require some hacking with the
.top file. I am a bit confused with that.

 

As mentioned previously, I did 

1)grompp -f new.mdp -n index.ndx -c old.tpr -o rerun1.tpr -p topol.top

Where energy groups are added to the .mdp file

 

2)tpbconv -s rerun1.tpr -n index -o rerun2.tpr -nsteps 0 

To produce individual .tpr files for i) System  ii) Protein only  iii)
Lipid only  iv) solv_ion only  v) protein_lipid  vi) protein_solv_ion
vii) lipid_solv_ion

 

3)trjconv  -f old.trr  -n index.ndx  -s rerun2.tpr  -o rerun.trr

To produce a corresponding .trr files for (i) to (vii) in step 2

 

4) mdrun -s rerun2.tpr  -rerun  rerun.trr

To produce the .edr files for (i) to (vii)

 

With g_energy I derived these different components for (i) to (vii):

a-PE

b-Angle (lipid)/G96 angle (protein)

c-Proper dihedral

d-Improper dihedral

e-Ryckaert Belleman (lipid)

f-LJ 1-4

g-Coul 1-4

h-L J SR

i-Disp. Corr

j-Coul SR

k-Coul. Recip

 

After doing some maths, I got the below:

1)      the PEs for (i) to (vii) are the sum of their corresponding b-k
values

e.g. the PE for the system =  Angle (lipid) + G96 angle (protein) +
Prop. Dihedral + imp dihedral + Ryckaert Belleman (lipid) + LJ14 + Coul
1-4 + LJ SR + Disp Corr + Coul SR + Coul recip. 

2)      For bonded terms such as proper dihedral, I found that the
proper dihedral of the system (i) = (ii) + (iii) i.e sum of proper
dihedral values from protein and lipid. The same goes for the improper
dihedral term. This is where I am confused. Doesn't this show that the
bonded terms can be decomposed?

 

I realize there has been a rather long pause between this email and the
last, sorry about that, it's because I have been trying to work out
myself why you say the bonded terms cannot be decomposed easily.

 

Perhaps as Justin pointed out, it is probably not worth the effort to
derive individual energies of the various components that make up my
system, but I would really appreciate a reply on this before I put the
matter to rest. Thanks!

 

HW

 

 

From: gmx-users-bounces at gromacs.org
[mailto:gmx-users-bounces at gromacs.org] On Behalf Of Mark Abraham
Sent: Friday, November 12, 2010 12:32 AM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] mdrun -rerun: bonded interactions of protein

 

On 11/11/2010 7:46 PM, NG HUI WEN wrote: 

Thanks a lot Justin and Mark for your useful input.

 

Indeed Justin was right, the quest to dissect the total energy of the
system to get that contributed by the protein alone was not trivial at
all! 


Indeed. If I'd observed you were doing PME, I'd have made the same
point. 




I missed out this thread yesterday
http://www.mail-archive.com/gmx-users@gromacs.org/msg34610.html as well
as Mark's trick to decompose bonded terms by hacking the .top file.
(However, I read from
http://www.gromacs.org/Documentation/Gromacs_Utilities/g_energy that the
Coul-recip term is not decomposable regardless of the tricks used...)

 

Mark, following your advice, I added -nsteps 0 (-1 didn't work) to
tpbconv. I managed to get the interactive prompt succesfully this time. 


Yes, nsteps=-1 is a GROMACS 4.5-ism, sorry. 




Reading toplogy and shit from rerun1.tpr
Reading file rerun1.tpr, VERSION 4.0.7 (single precision)
Setting nsteps to 0
Group     0 (      System) has 100581 elements
Group     1 (     Protein) has  4002 elements
Group     2 (   Protein-H) has  3145 elements
Group     3 (     C-alpha) has   412 elements
Group     4 (    Backbone) has  1236 elements
Group     5 (   MainChain) has  1649 elements
Group     6 (MainChain+Cb) has  2027 elements
Group     7 ( MainChain+H) has  2039 elements
Group     8 (   SideChain) has  1963 elements
Group     9 ( SideChain-H) has  1496 elements
Group    10 ( Prot-Masses) has  4002 elements
Group    11 ( Non-Protein) has 96579 elements
Group    12 (        POPE) has 13052 elements
Group    13 (         SOL) has 83523 elements
Group    14 (         CL-) has     4 elements
Group    15 (       Other) has 96579 elements
Group    16 (     SOL_CL-) has 83527 elements
Group    17 (Protein_POPE) has 17054 elements
Select a group: 

 

As suggested, I have also tried to compare the (average) values of
subset.tpr and fullsystem.tpr, these are the total energies of:

System=              -1.28209 e+06

Protein=              -9571.86

Lipid=                    -296121

SOL_CL=              -821353

Other=                 -1.22684e+06

It was found that the: 

Total energies (Protein + lipid + SOL_CL-) not equal to total energy of
the system

Total energies of (Lipid + SOL_CL-) not equal to total energy of other

Total energies of (others + protein) not equal to total energy of the
system

 

I have yet to work out the reasons for the discrepancies, will spend
some time pondering and trying to make sense of these (and other)
values.


Aren't you omitting the inter-group non-bonded terms? Anyway, I was
suggesting this comparison for seeing whether the bonded components can
be decomposed group-wise. It is already known that the non-bonded
decomposition works.

Mark




 

Thanks a bunch again!

 

HW

 

 

From: gmx-users-bounces at gromacs.org
[mailto:gmx-users-bounces at gromacs.org] On Behalf Of Mark Abraham
Sent: Wednesday, November 10, 2010 9:08 PM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] mdrun -rerun: bonded interactions of protein

 

On 10/11/2010 11:29 PM, NG HUI WEN wrote: 

Hi Gmxusers,

 

I have been trying to run mdrun -rerun to get the energy of the protein
in my protein-lipid system. I know similar questions have been raised on
this topic before, I have tried to glean useful information from them to
solve my problem but unfortunately to no avail. Thanks for your
patience!  

 

As I did not have energygrps in the initial .mdp file, I included it
this time  

 

integrator  = md        
nsteps      = 0   
dt          = 0.002           
nstxout     = 50000           
nstvout     = 50000           
nstenergy   = 500       
nstlog      = 500       
Continuation      = yes       
constraint_algorithm = lincs  
constraints = all-bonds 
lincs_iter  = 1         
lincs_order = 4         
ns_type           = grid            
nstlist           = 5         

rlist       = 1.2       
rcoulomb    = 1.2       
rvdw        = 1.2       
coulombtype = PME       
pme_order   = 4         
fourierspacing    = 0.16            
tcoupl            = Nose-Hoover           
tc-grps           = Protein POPE    SOL_CL-     

tau_t       = 0.1 0.1   0.1   
ref_t       = 323       323   323   
pcoupl            = Parrinello-Rahman     
pcoupltype  = semiisotropic         

tau_p       = 5.0             
ref_p       = 1.0 1.0         
compressibility = 4.5e-5      4.5e-5      
pbc         = xyz       
DispCorr    = EnerPres  
gen_vel           = no        
nstcomm         = 1
comm-mode       = Linear
comm-grps       = Protein_POPE SOL_CL-
energygrps  = Protein SOL POPE

 

I then did 

1)grompp -f new.mdp -n index.ndx -c old.tpr -o rerun.tpr -p topol.top

2)trjconv -f old.trr -n index.ndx -s rerun.tpr -o rerun.trr (when
prompted, I selected "0" system)

3)mdrun -s rerun.tpr -rerun rerun.trr

 

I notice that the previous post
http://oldwww.gromacs.org/pipermail/gmx-users/2009-January/038968.html
suggested to use tpbconv (on rerun.tpr) and trjconv (on rerun.trr) to
extract the protein only. While it was possible to do so with trjconv,
it wasn't feasible with tpbconv (I'm using gromacs 4.0.7) - I might have
missed out something as I did not get any prompt/output (see below)

 

tpbconv -s topol.tpr -n index_P.ndx -o rerun2.tpr

 

Reading toplogy and shit from topol.tpr
Reading file topol.tpr, VERSION 4.0.7 (single precision)
0 steps (0 ps) remaining from first run.
You've simulated long enough. Not writing tpr file


Looking at the code, tpbconv checks for whether there are any more steps
to simulate before even considering letting you use it in the "create
subset" mode. You could argue that this is buggy, because the way it
computes whether there are any more steps will always indicate no more
steps in such cases. However, the workaround is to use tpbconv -nsteps
-1 (as well as the other stuff). Let me know how this goes and I'll
update the documentation.





 

Using the g_energy command on the output energy.edr file, I got among
others, these options to choose

49  Coul-SR:Protein-Protein             50  LJ-SR:Protein-Protein


51  Coul-14:Protein-Protein             52  LJ-14:Protein-Protein


 

In order to get the energy of the protein, I reckon I have to add
49,50,51,52 (to account for the nonbonded components) and


I think that you can make a case either way for 1-4 interactions -
they're algorithmically similar to the other non-bonded interactions,
but their parameter values should be tightly coupled to some other of
the bonded parameters, but then in several forcefields those values are
just scaled versions of the normal ones... I'd guess most people call
them non-bonded.





1  Angle            2  G96Angle         3  Proper-Dih.      4
Ryckaert-Bell.  5  Improper-Dih  

for the bonded components. However, I think 1-5 is the bonded terms for
the system and not the protein alone. Can anyone help me with this?


Compare values from reruns on the subset-tpr and the full-tpr. I don't
know whether/how well that works. In extremis, you should be able to
reduce the .tpr to a single interaction.





 

 

Also, on a slightly different note, this post
http://oldwww.gromacs.org/pipermail/gmx-users/2005-July/016307.html
suggested that the force constant of the solvent (DMSO) to be adjusted
to zero. Am I right to think that it does not apply to my case as my
protein-lipid system is solvated with SPC? (I read that SPC is rigid
water. I did not add -DFLEXIBLE in .mdp)


That was in the context of a full-tpr rerun. The point is that atoms
that have been constrained don't contribute to relevant sums. Because
you are using constraints, there's no "Bonds" energy sum for any atom
pair. Because the DMSO was a flexible model, its bonded-interactions
contributions need to be subtracted (i.e. parameters set to zero) to get
a group-wise bonded-interaction value.

Mark





 

Any help would be highly appreciated.Thanks!

 

HW

<< 

This message and any attachment are intended solely for the addressee
and may contain confidential information. If you have received this
message in error, please send it back to me, and immediately delete it.
Please do not use, copy or disclose the information contained in this
message or in any attachment. Any views or opinions expressed by the
author of this email do not necessarily reflect the views of the
University of Nottingham. 

This message has been checked for viruses but the contents of an
attachment may still contain software viruses which could damage your
computer system: you are advised to perform your own checks. Email
communications with the University of Nottingham may be monitored as
permitted by UK & Malaysia legislation.

>> 

 

<< 

This message and any attachment are intended solely for the addressee
and may contain confidential information. If you have received this
message in error, please send it back to me, and immediately delete it.
Please do not use, copy or disclose the information contained in this
message or in any attachment. Any views or opinions expressed by the
author of this email do not necessarily reflect the views of the
University of Nottingham. 

This message has been checked for viruses but the contents of an
attachment may still contain software viruses which could damage your
computer system: you are advised to perform your own checks. Email
communications with the University of Nottingham may be monitored as
permitted by UK & Malaysia legislation.

>> 

 

<< This message and any attachment are intended solely for the addressee and may contain confidential information. If you have received this message in error, please send it back to me, and immediately delete it. Please do not use, copy or disclose the information contained in this message or in any attachment. Any views or opinions expressed by the author of this email do not necessarily reflect the views of the University of Nottingham. 

This message has been checked for viruses but the contents of an attachment may still contain software viruses which could damage your computer system: you are advised to perform your own checks. Email communications with the University of Nottingham may be monitored as permitted by UK & Malaysia legislation. >>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101122/f5dedfd3/attachment.html>


More information about the gromacs.org_gmx-users mailing list