Nicolas SAPAY wrote:
Hello all gmx users,
I'm doing some tests on the CHARMM TIP3P water model. After a several
consultations of the GROMACS mailing-list and doc, I still have some
questions.
First of all, I have used a tip3p.itp file were the TIP3 rigid water model
is defined :
[ moleculetype ]
;name nrexcl
TIP3 1
[ atoms ]
;id atom_type res residu_name atom_name charge_group
charge
1 OT 1 TIP3 OH2 1
-0.834
2 HT 1 TIP3 H1 1
+0.417
3 HT 1 TIP3 H2 1
+0.417
[ settles ]
; the SETTLE ALGORITHM is an analytical solution of SHAKE specifically for
water
;oxygen func d(O,H) d(H,H)
1 1 0.09572 0.15139
[ exclusions ]
; for non bonded interactions exclusion
1 2 3
2 1 3
3 1 2
Then, I have defined a TIP3 in a ffcharmm.rtp
[ TIP3 ]
[ atoms ]
;name type charge charge_group
OH2 OT -0.8340 1
H1 HT +0.4170 1
H2 HT +0.4170 1
And finally, I have defined the LJ parameters in th ffcharmmnb.itp
HT 1.0080 0.000 A 3.153970e-09 1.292128e-17
OT 15.9994 0.000 A 2.489540e-03 2.434767e-06
That all seems to be in order. My LJ parameters have only 3 sig figs,
however.
This parameters are consistent with the implementation of Mark Abraham.
Everything works fine when I build the .top topology with pdb2gmxfile
except that the TIP3 residues in my pdb are not recognized as water and
are thus listed in the [atoms] section (the tip3p.itp is included by an
#include
command). I can fix that by creating a HOH or a SOL residue with same
property than TIP3 but I'm not sure it is a good idea.
All I did was hack tip3p.itp to include an #ifdef for CHARMM that had
the appropriate atom types and names and a SOL residue name.
As far as I have
understood, GROMACS uses some specific code to improve the non-bonded
interactions computation of water. This code is designed for water models
without LJ parameters and I don't want to use it (since TIP3P has got LJ
paramters in CHARMM).
Mostly correct - non-CHARMM TIP3P has LJ on the oxygen. However, the
decision about which code to use is made on-the-fly by mdrun. If the
molecule matches the pattern of a water molecule and only has LJ
parameters on the oxygen then the optimized inner loops are used for
water-nonwater and water-water interactions. The optimization,
obviously, omits calculating LJ interactions with hydrogens.
You can see what is happening by looking at the accounting at the end of
.log files... lines like "Coul(T) + LJ [W3-W3]" indicate the time spent
calling the innerloop routine that does coulomb lookups on a table,
computes LJ and does so for water-water interactions. If the description
of CHARMM TIP3P water is correct, you should see no lines with the [W3]
indicators, and thus the LJ terms are being treated correctly.
As an aside, I did develop assembly innerloops for CHARMM TIP3P water
for x86 some time back and posted here. My main compute platform isn't
x86 though, so I've not used them in anger. See discussions back in Feb
and March this year on gmx-developers.
But in this case, I not sure that my TIP3 residues
will be considered as rigid water model and then that the H1...H2
distance will be constrained (I don't use the -DFLEXIBLE flag in my .mdp
file and use the -I flag ot include the directory where the CHARMM force
field is defined).
Again, you can see in the .log file whether Settle is being used
appropriately. I can't tell you whether your setup above will do this
correctly, since I've not done it. My approach seems to work.
Mark
_______________________________________________
gmx-users mailing list [email protected]
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [EMAIL PROTECTED]
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
