[gmx-users] mdrun -rerun: bonded interactions of protein
Justin A. Lemkul
jalemkul at vt.edu
Mon Nov 22 20:29:06 CET 2010
Quoting NG HUI WEN <HuiWen.Ng at nottingham.edu.my>:
> 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.
>
>
The bonded terms can be decomposed by the procedure you've gone through. I (and
others) always caution against thinking that the simple use of energygrps will
do the trick. An iterative approach is indeed required.
>
> 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!
>
I still do not see the value of trying to claim "the energy of the protein is
(some number)" by decomposing various potential energy terms. The PME term
cannot be decomposed, for one. Even though you've gotten a value of Coul-recip
for your various calculations, I don't know how valid the results are.
-Justin
>
>
> 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. >>
>
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list