[gmx-users] mdrun -rerun: bonded interactions of protein
Mark Abraham
Mark.Abraham at anu.edu.au
Thu Nov 11 17:32:20 CET 2010
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:
>
> 1)System= -1.28209 e+06
>
> 2)Protein= -9571.86
>
> 3)Lipid= -296121
>
> 4)SOL_CL= -821353
>
> 5)Other= -1.22684e+06
>
> It was found that the:
>
> 1)Total energies (Protein + lipid + SOL_CL-) not equal to total energy
> of the system
>
> 2)Total energies of (Lipid + SOL_CL-) not equal to total energy of other
>
> 3)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.
>
> >>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101112/ca3b7d2e/attachment.html>
More information about the gromacs.org_gmx-users
mailing list