OPC water from CHARMM-GUI

User discussions
1

GROMACS version: 2021.4
GROMACS modification: No

Dear gmx users,

The force field file for the solvated protein system was obtained from CHARMM-GUI’s Solution Builder.

protein = FF19SB
water = OPC
ligand = GAFF2

Comparing OPC.itp/forcefield.itp obtained from CHARMM-GUI with water files (spc.itp, tip3p.itp, tip4p.itp) of various forcefields bulit in GROMACS,

There is no information about [ bonds ] and [ angles ] in the OPC.itp file, and [ angletypes ] in the forcefield.itp file.

Is MD simulation correct without these information?

Also, when I obtained force field files for tip3p water instead of opc water, there is no information about [ bonds ], [ angles ], and [ angletypes ].

How can this be solved?

Thanks for your help.

--------------- OPC.itp ---------------
[ moleculetype ]
; name nrexcl
OPC 1

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 OW 1 OPC O 1 0.000000 16.0000 ; qtot 0.000
2 HW 1 OPC H1 2 0.679142 1.0080 ; qtot 0.679
3 HW 1 OPC H2 3 0.679142 1.0080 ; qtot 1.358
4 EP 1 OPC EPW 4 -1.358284 0.0000 ; qtot 0.000

[ settles ]
; OW funct doh dhh
1 1 8.724331e-02 1.371205e-01

[ virtual_sites3 ]
; Vsite from func a b
4 1 2 3 1 1.477224e-01 1.477224e-01

[ exclusions ]
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3

--------------- forcefield.itp ---------------
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.500000 0.833333

[ atomtypes ]
; name at.num mass charge ptype sigma epsilon ; sigma_14 epsilon_14
EP 0 0.0000 -1.358 A 1.78179743628e-01 0.000000e+00
HW 1 1.0080 0.679 A 0.00000000000e+00 0.000000e+00
OW 8 16.0000 0.000 A 3.16655208196e-01 8.903586e-01

[ bondtypes ]
; i j func b0 Kb
HW HW 1 1.371205e-01 4.627504e+05
OW EP 1 1.593983e-02 4.627504e+05
OW HW 1 8.724331e-02 4.627504e+05

[ pairtypes ]
; i j func sigma1-4 epsilon1-4

2

Hi,

the [ settles ] directive means that the OPC water molecule is constrained using the SETTLE algorithm, which is an analytical solution to constraining water molecule atoms. Thus no bonds need to be specified. Same goes for TIP3P, SPC/E, etc.

See this page for more information, including a link to the paper introducing the algorithm.

Regards,
Petter

1 Like
3

It is still a difficult problem for me.

Q1.
No need to set #ifndef FLEXIBLE and #else in itp file for water?
If [ settles ] is set in #ifndef FLEXIBLE, it seems that information about [ bonds ] and [ angles ] should be included in the #else part.

Q2.
#ifndef FLEXIBLE
[ settles ]
Does this mean that OPC water is processed by the settle algorithm when -DFLEXIBLE?

#else
Does this part mean when -DPOSRES condition?

4

All the water models commonly used in GROMACS are intended to be rigid, via SETTLE. Making them flexible is rarely used and in the case of CHARMM-GUI, they are assuming you are using conventional approaches. In this topology, -DFLEXIBLE will do nothing, and it has nothing to do with restraints, which are topologically different from constraints (SETTLE).

5

In the minimization mdp file obtained from CHARMM-GUI, define is set to -DPOSRES/-DPOSRES_FC_SC/-DPOSRES_FC_BB instead of -FLEXIBLE.

Also, the following error seems to occur because there is no setting such as #ifndef FLEXIBLE in the water molecule file.
.
.
.
If define = -DFLEXIBLE in mdp file,

WARNING 1 [file topol.top, line 20]:
The following macros were defined in the ‘define’ mdp field with the -D
prefix, but were not used in the topology:
FLEXIBLE
If you haven’t made a spelling error, either use the macro you defined,
or don’t define the macro
.
.
If define = -DPOSRES,

WARNING 1 [file topol.top, line 20]:
The following macros were defined in the ‘define’ mdp field with the -D
prefix, but were not used in the topology:
POSRES_FC_SC
POSRES_FC_BB
POSRES
If you haven’t made a spelling error, either use the macro you defined,
or don’t define the macro
.
.
.
If I delete the define part of the mdp file, grompp proceeds without error.
Is it ok to delete this?

Then, how to solve the error about -DPOSRES during equilibration?

6

This is because restraints and constraints are separate concepts, and as noted before, CHARMM-GUI does not assume that the user will every need/want flexible water. So they use their conventions for typical use.

This suggests a problem with the input file and you should notify the CHARMM-GUI developers that they are providing inputs that generate warnings. The typical workflow from CHARMM-GUI is a progressive relaxation of restraints on different parts of the structure, but those may or may not be relevant or useful in all cases.

If you don’t need those restraints, sure.

You’ll need to generate a suitable restraint .itp file and enclose it within an #ifdef POSRES...#endif block in the topology.

7

So far, in the minimization and production run stages,
I set define = -DFLEXIBLE in the mdp file, is there no need for that?

I saw the word ‘flexible water’ in the manual and thought that I should set it to -DFLEXIBLE.

However, from the answers here, I found out that rigid water is generally applied.

So, in general, there is no need to set define in the minimization and production run stages, and is it correct to set -DPOSRES only in the equilibration stage?

Thank you

8

Yes, the conventional approach to equilibration is to restrain the solute and allow the solvent to relax around it.

9

I am a bit confused about how to correctly implement OPC4 water in gromacs

The OPC4.itp file I have is:

“”"
[ atomtypes ]
OW OW 0.0000 0.0000 A 3.16655e-01 8.903586e-01
HW HW 0.0000 0.0000 A 0.00000e+00 0.00000e+00
MW MW 0.0000 0.0000 A 0.00000e+00 0.00000e+00

[ moleculetype ]
; molname nrexcl
SOL 2

[ atoms ]
; id at type res nr res name at name cg nr charge mass
1 OW 1 SOL OW 1 0 16.00000
2 HW 1 SOL HW1 1 0.67914 1.00800
3 HW 1 SOL HW2 1 0.67914 1.00800
4 MW 1 SOL MW 1 -1.35828 0.00000

#ifndef FLEXIBLE

[ settles ]
; i funct doh dhh
1 1 0.08724 0.13712

#else

[ bonds ]
; i j funct length force.c.
1 2 1 0.08724 502416.0 0.08724 502416.0
1 3 1 0.08724 502416.0 0.08724 502416.0

[ angles ]
; i j k funct angle force.c.
2 1 3 1 103.6 628.02 103.6 628.02

#endif

[ virtual_sites3 ]
; Vsite from funct a b
4 1 2 3 1 0.147722363 0.147722363

[ exclusions ]
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3

; The position of the virtual site is computed as follows:
;
; O
;
; V
;
; H H
;
; Ewald OPC:
; const = distance (OV) / [ cos (angle(VOH)) * distance (OH) ]
; 0.01594 nm / [ cos (51.8 deg) * 0.0872433 nm ]
; then a = b = 0.5 * const = 0.147722363
;
; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1)
“”"

OP seems to have an OPC file without the bonds and angles description, and in the OPC4 paper it says they didn’t impose any geometry constraints on the model other than the symmetry and optimized the distribution of point charges (https://pubs.acs.org/doi/10.1021/jz501780a).

Should I be using the [settles] section of this .itp file or the [bonds] [angles] section?

10

OPC is rigid, like all other popular simple water models,so you should use SETTLE. In GROMACS, bonds and angles are only there to support energy minimization which do not support constraints.