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.ne...@utoronto.ca Subject: [gmx-users] Half double pair list method in GROMACS To: gmx-users@gromacs.org Message-ID: <20110901143853.537ln0dj94w4w...@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 -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists