Solvent partition

GROMACS version: 2020.1
GROMACS modification: No
Hi all,

I am working with a water system, where some waters around the system’s COM will be flat-bottom restrained, while the rest will be unrestrained. I used CHARMM-GUI to create a water box, which was equilibrated during 1 ns. Then, following a previous suggestion (Setting spherical restraints), I modified the topology of the system to create two topologies for the restrained (TIP3) and unrestrained (TP3R, the name is also changed in the input .gro file) solvents, grouped in an index file.

When I run mdrun, get this error:

Fatal error:
The [molecules] section of your topology specifies more than one block of
a [moleculetype] with a [settles] block. Only one such is allowed.
If you are trying to partition your solvent into different *groups*
(e.g. for freezing, T-coupling, etc.), you are using the wrong approach. Index
files specify groups. Otherwise, you may wish to change the least-used
block of molecules with SETTLE constraints into 3 normal constraints.

If I comment the [ settles ] and [ exclusions ] directives in either molecule, then:

starting mdrun '10A spherical restrained water in a box'
1000000 steps,   1000.0 ps.
step 0
step 1: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
Segmentation fault (core dumped)

Bellow are the .itp files and the topology.

Unrestrained solvent:

[ moleculetype ]
; name	nrexcl
TIP3	     2

[ atoms ]
; nr	type	resnr	residu	atom	cgnr	charge	mass
     1         OT      1     TIP3    OH2      1    -0.834000    15.9994   ; qtot -0.834
     2         HT      1     TIP3     H1      2     0.417000     1.0080   ; qtot -0.417
     3         HT      1     TIP3     H2      3     0.417000     1.0080   ; qtot  0.000

[ settles ]
 OW	funct	doh	dhh
    1     1  9.572000e-02  1.513900e-01

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

Restrained solvent:

[ moleculetype ]
; name	nrexcl
TP3R	     2

[ atoms ]
; nr	type	resnr	residu	atom	cgnr	charge	mass
     1         OT      1     TIP3    OH2      1    -0.834000    15.9994   ; qtot -0.834
     2         HT      1     TIP3     H1      2     0.417000     1.0080   ; qtot -0.417
     3         HT      1     TIP3     H2      3     0.417000     1.0080   ; qtot  0.000

[ settles ]
 OW	funct	doh	dhh
  1     1  9.572000e-02  1.513900e-01

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

#ifdef POSRES_WATER_SPH
[ position_restraints ]
;	atom	type	g	r	k
	1	    2	    1	1	10000
#endif

“Mini” force field:

[ defaults ]
; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
	   1	        2	      yes	    1.0	    1.0

[ atomtypes ]
; name	at.num	mass	charge	ptype	sigma	epsilon	;	sigma_14	epsilon_14
      HT     1     1.0080      0.417     A    4.00013524445e-02    1.924640e-01 
      OT     8    15.9994     -0.834     A    3.15057422683e-01    6.363864e-01 

[ bondtypes ]
; i	j	func	b0	Kb
     HT      HT     1  1.513900e-01  0.000000e+00
     HT      OT     1  9.572000e-02  3.765600e+05

[ pairtypes ]
; i	j	func	sigma1-4	epsilon1-4
     HT      HT     1  4.00013524445e-02  1.92464000000e-01 
     HT      OT     1  1.77529387564e-01  3.49973530556e-01 
     OT      OT     1  3.15057422683e-01  6.36386400000e-01 

[ angletypes ]
; i	j	k	func	th0	Kth	s0	Kub
     HT      OT      HT     5  1.0452000e+02  4.6024000e+02  0.0000000e+00  0.0000000e+00

Topology:

; Include forcefield parameters
#include "charmm36-jul2020.ff/forcefield.itp"

; Include forcefield parameters for TIP3 water
#include "charmm36-jul2020.ff/TIP3.itp"

; Include forcefield parameters for TIP3 water restrained
#include "charmm36-jul2020.ff/TIP3_restrained.itp"

[ system ]
; Name
10A spherical restrained (10000 kJ-mol_nm2) water in a box

[ molecules ]
; Compound	#mols
TIP3  	    6704
TP3R		137

Attached is the .mdp file:
md_sphere_restrained.mdp (1.2 KB)

Am I missing something? Any help will be appreciated.

Best


The run is able to start, only if I comment the [ settles ] directive in TP3R (TIP3_restrained.itp), but the flat-bottom restrains are not applied.

You will have to delete [settles] from one of the topologies and replace it with a triangular set of constraints, e.g.

[ constraints ]
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139

Thanks a lot! It worked.

Hi @jalemkul, this thread is for flat-bottom restraint. I am using position restraints, can I still use

[ constraints ]
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139

to replace

[ settles ]
; OW	funct	doh	dhh
1       1       0.09572 0.15139

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

?