[gmx-users] Half double pair list method in GROMACS [update]
ABEL Stephane 175950
Stephane.ABEL at cea.fr
Sat Dec 3 14:54:25 CET 2011
Hi Sai
I am also agree with Christ, your approach seems also correct to me. But don't to forget to test your change to be sure, by doing a single point energy test as it is said in the Christ's message.
Good luck
Stephane
----------------------------------------------------------------------
Message: 1
Date: Fri, 02 Dec 2011 16:07:19 -0500
From: Chris Neale <chris.neale at utoronto.ca>
Subject: [gmx-users] Half double pair list method in GROMACS [update]
To: gmx-users at gromacs.org
Message-ID: <4ED93E07.5070401 at utoronto.ca>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
What you have done seems alright. I didn't look closely enough to be
sure though. One of the good things about this method is that you can
easily test it yourself. To do this, create two different .gro files,
one containing the atoms from one ff and the other containing the other
atoms. For each, do separate zero-step mdrun evaluations of their
energies under (a) their original ff, and (b) the combined HEDP ff that
you have constructed. The energy evaluations should be the same once you
account for rounding errors. Then all that is left is for you to be
certain that you didn't cause any problems for nonbonded interactions.
Note that you'll need to return to the original atomtypes for the first
half of this test.
Note that the only thing that the HEDP method is intended to perturb is
the 1-4 interactions (and by perturb I mean that they will now be correct).
Chris.
-- original message --
Hi Stephane & Chris,
I followed all the threads posted by you two.
I have a protein using ff99SB and GLYCAM for sugars. I have a disaccharide
bound to protein. In xleap of AMBERTOOLS, I use the GLYCAM and ff99SB to
generate the topology and coordinate files. I did the tests as Chris'
original posts and by Stephane.
The 1-4 interaction terms match for sugar alone and protein alone with HEDP
method.
When I generate topology and coordinate files with xleap of AMBER, there
are three atom types that are common to protein and sugar. The atom types
for my case are H1, OH and HO.
Since for protein pairtypes using ff99SB, the epsilon has to be divided by
10 and pairs section be replicated five times and for sugar the epsilon in
the pairtypes be divided by six and replicated six times, I am a little
concerned about the three atomtypes that are common.
So what I did was to change the atomtypes of sugar as H1S, OHS HOS for
sugar and H1, OH and HO for protein and I made the changes accordingly in
the pairtypes section for protein and sugar.
Is this a valid approach?
Any suggestion will be helpful.
To make things little more clear:
H1 H1 0.0000 0.0000 A 2.47135e-01 6.56888e-02
;originally obtained using amb2gmx.pl
H1S H1S 0.0000 0.0000 A 2.47135e-01 6.56888e-02 ;Glycam
Hydrogen of Sugar ( I changed this so that the common atom types be
separated)
In the ffnonbonded_complex_mod.itp:
;;using combination rule of 2
[ pairtypes ]
;;for protein
H1 H1 1 0.247135000 0.006568880 ;the epsilon is divided
by 10
;;for sugar
H1S H1S 1 0.24713500 0.010948133 ;the epsilon is divided
by 6
Thanks for your time,
Regards
Sai
On Mon, Sep 5, 2011 at 11:33 AM, ABEL Stephane 175950
<Stephane.ABEL at cea.fr>wrote:
> Dear All,
>
> Below a little update and results about the application of half double
> pair list method to scale properly the Coulombic 1-4 interactions in case
> of a system where the AMBER99SB (fudgeLJ=0.5 and fudgeLJ=0.8333) and
> GLYCAM06 (fudgeLJ=1.0 and fudgeLJ=1.0) force fields are combined.
>
> I have followed the 4 steps described in [1] and used the following
values
> in my forcefield.itp file
>
> [ defaults ]
> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> 1 2 yes 1.0 0.16666666666666
> #include "ffnonbonded_mod.itp"
> ;#include "ffnonbonded.itp"
> #include "ffbonded.itp"
>
> I used two different topology files for the glycolipid (bDM) and the
> peptide. One with (*_mod.itp) with pair list parameters duplicated 6
times
> (bDM) and 5 times (peptide) and with single pair list (*_no_mod.itp) as
> decribed in [1].
>
> TESTING:
>
> Three 3 different systems were examined:
>
> A. A first system containing 1 glycolipid (bDM) in water cubic box
> B. A second system with 1 peptide in TIP3P water
> C. And a third system with 1 peptide and 1 glycolipid in water cubic box
>
> To obtain the glycolipid and peptide energy pairs, I did one step of
MD in
> NVT ensemble with the *.mdp file given in [2] with different
energygrps and
> tc_grps.
> For 1. energygrps and tc_grps = bDM SOL
> For 2. energygrps and tc_grps = Protein SOL
> For 3. energygrps and = Protein bDM SOL
>
> bDM/water system
>
> Test_A1
>
> ## Control with GLYCAM force field fudgeLJ fudgeQQ parameters and the
> *_no_mod.itp file :
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 1.0 1.0
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> bDM-bDM -4.00855e+02 -3.71401e+01 2.03406e+03 2.79234e+02
>
> Test_A2
>
> #### with the topology *_mod.itp file and the directive
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 1.0
> 0.16666666666666
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> bDM-bDM -4.00855e+02 -3.71401e+01 2.03406e+03 2.79234e+02
>
> Test_A3
>
> #### with the topology *_no_mod.itp file and the directive
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 1.0
> 0.16666666666666
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> bDM-bDM -4.00855e+02 -3.71401e+01 3.39010e+02 2.79234e+02
>
> Coul-14 energy for bDM-bDM in Test_A1 = Coul-14 energy in Test_A2 -->
OK !
> Coul-14 energy for bDM-bDM Test_A3 is 6 smaller than Coul-14 energy in
> Test_A1 and Test 2 ---> OK !
>
> peptide/water system
>
> Test_B1
>
> #### Control with AMBER force field fudgeLJ fudgeQQ parameters, the
> *_no_mod.itp file
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 0.5 0.8333
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> Protein-Protein -1.49026e+03 -4.12114e+02 3.80551e+03
6.55321e+02
>
> Test_B2
>
> #### Control with the peptide *_no_mod.itp file and the directive
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 1.0 0.16666666666
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> Protein-Protein -1.49026e+03 -4.12114e+02 7.61132e+02
1.31064e+03
>
> Test_B3
>
> #### With the peptide *_mod.itp file and the directive
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 1.0
> 0.16666666666666
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> Protein-Protein -1.49026e+03 -4.12114e+02 3.80566e+03
6.55320e+02
>
> Coul-14 and LJ-14 energies for Protein-Protein in Test_B3 are 5 times
> Coul-14 and 2 times LJ-14 energies, respectively, than in Test_B1
---> OK !
> Coul-14 and LJ-14 energies for Protein-Protein in Test_B3 = Coul-14 and
> LJ-14 energies in Test_B1 ---> OK !
>
> peptide/bDM/water system
>
> Test_C1
>
> ## Control with the peptide and the bDM *_mod.itp topology files and the
> directive
> ##; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> ## 1 2 yes 1.0
0.16666666666666
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> Protein-Protein -1.48386e+03 -3.99146e+02 3.79532e+03
6.53553e+02
> bDM-bDM -3.69788e+02 -3.51804e+01 2.00228e+03
2.25026e+02
>
> Test_C2
>
> ## Control with the peptide and the bDM *no_mod.itp topology files and
> the directive
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> Protein-Protein -1.48386e+03 -3.99146e+02 7.59063e+02
1.30711e+03
> bDM-bDM -3.69788e+02 -3.51804e+01 3.33713e+02
2.25026e+02
>
> bDM-bDM Coul-14 energy in Test_C2 is 6 times smaller, respectively, than
> in the Test_C1 ---> OK !
> Protein-Protein Coul-14 and LJ1-4 energies are in Test_C2 is 5 and 2
times
> smaller, respectively, than in the Test_C1 ---> OK !
>
> Conclusions
>
> - The trick works, great and energy difference between extremely small
> (negligible) in my case.
>
> [1] http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html
> [2]
> http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html
>
> Let me know if you have any comments or correct me if I am wrong.
>
> Thank you again.
>
> St?phane
>
>
>
> Chris,
>
> Thank you for your confirmation. I did the changes. I am currently doing
> some tests, I will send you a feedback about the results off-list (if you
> want) shortly.
>
> A bient?t
>
> St?phane
> ________________________________________
>
>
>
> Message: 1
> Date: Thu, 01 Sep 2011 14:38:53 -0400
> From: chris.neale at utoronto.ca
> Subject: [gmx-users] Half double pair list method in GROMACS
> To: gmx-users at gromacs.org
> Message-ID: <20110901143853.537ln0dj94w4wc4c at webmail.utoronto.ca>
> Content-Type: text/plain; charset=ISO-8859-1; DelSp="Yes";
> format="flowed"
>
> It's a typo, but it's in the discussion and not in the "do this" part
> of the method so I decided not to mention it. I don't see another
> question in this post, so I hope that you have figured things out.
> Note that I have never tested the exact implementation that I
> suggested in that April email. It seems like it should work just fine,
> but it is numerically different than the OPLS/Berger combination so
> there is no way to be sure until you check the energies as I
> suggested. I'd be interested to have you report back on the energy
> matching once you have done the tests.
>
> Good luck,
> Chris.
>
> -- original message --
>
> Hi Chris,
>
> Sorry to repost the same question, but I have really tested your
method the
> last few weeks. My question about the gen-pairs directive come from the
> fact
> that I have read a message from you
>
> http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html
>
> Where you detailed how to use the Berger and OPLS force fields
together in
> the same system. By reading carefully the meaning of the gen-pairs
> directive, I found several errors in my force field.
>
> BYW in your previous message
>
> http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html
>
> In the step 3, I think there is a typo, it is "values/10" instead of
> "values/12". Am I right?
>
> Thank you again
>
> Stephane
>
>
>
> Dear Stephane:
>
> We discussed this in April:
>
> http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html
>
> At that time I also provided a method for you to verify your files
> (and the method in general).
>
> It is possible for you to answer your gen-pairs question by looking
> into the manual and reading about pairs, gen-pairs, and pairtypes. I
> think that this is one case where you would benefit more from fully
> understanding how these parts work than from a direct answer to your
> question.
>
> If you are having problems, please provide a whole bunch more
> information on the problems that you are seeing. If, on the other
> hand, you are just looking for somebody other than me to comment on
> the accuracy of the April post, then that is perfectly fine with me,
> but you should state that.
>
> Chris.
>
>
> -- original message --
>
> Dear All,
>
>
>
> I try to apply the half double pair list method for a system containing a
> glycolipid surfactant and a peptide modeled with the GLYCAM and AMBER99SB
> force field. Briefly what I did :
>
>
>
> 1- I have changed the forcefield.itp like this
>
>
>
> [ defaults ]
>
> ; gen_pairs set explicitly ---> gen-pairs = "no"
>
> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
>
> ;1 2 yes 0.5 0.8333 ---> for
> AMBER99SB only
>
> 1 2 no 1.0 0.16666666666666
> ---
> > for both GLYCAM et AMBER99SB
>
> ;1 2 yes 1.0 1.0 --> for
> GLYCAM only
>
>
>
> #include "ffnonbonded_mod.itp"
>
> #include "ffbonded.itp"
>
>
>
> 2- For the surfactant, I have calculated the pair-types interactions
> manually with the comb-rule = 2 and divided the values /6 and added these
> values in [pairtypes] section in the ffnonbonded_mod.itp files
>
>
>
> 3- For the peptide, I have calculated the pair-types interactions
manually
> comb-rule = 2 and divided the values /10 and added these values in
> [pairtypes] section in the ffnonbonded_mod.itp files.
>
>
>
> 4- In the surfactant topology file I have repeated 6 times the [pairs]
> directives 0.166666*6 ~=10
>
>
>
> 5 - In the peptide topology file I have repeated 5 times the [pairs]
> directives 0.166666*5 ~= 0.8333
>
>
>
> Is it correct ?
>
>
>
> However I have a little doubt about "gen-pairs" directive should I
have set
> it to "no" or "yes". in a previous message with a similar problem,
the gen
> directive was set "yes"
>
>
>
> http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html
>
>
>
> Thank you for your help
>
>
>
> Stephane
>
More information about the gromacs.org_gmx-users
mailing list