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

[Mark Abraham]

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.
>
> >> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101111/8e5b80b6/attachment.html>