[gmx-users] Using the CHARMM TIP3P water model in Gromacs

[Mark Abraham]

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