Unstable simulation

[krlitros87]

Hi,
I'm trying to parametrize a small molecule.
I was able to create both, the .itp and .gro files without issues. After this I first minimize in vaccum the molecule and then solvate 10 copies in a 10x10x10 solvent box, where finally I I minimize the system (10 molecules in a non-polarizable water box).
Sadly, I'm having problems equilibrating the system, where the simulation crash (Too many LINCS warnings (1002))
Do you think you can send me a copy of an .mdp file you used in your simulations?
Here are my .itp and .mdp files:
;;;; GENERATED WITH auto-martini
; INPUT SMILES: FC(F)(F)c1ccc(cc1)OC(CCNC)c2ccccc2
; Tristan Bereau (2014)

[moleculetype]
; molname nrexcl
GUA 2

[atoms]
; id type resnr residu atom cgnr charge smiles
1 N0 1 GUA N01 1 0 ; FCF
2 SC5 1 GUA S01 2 0 ; [c]1[c][c][c][c][c]1
3 SNa 1 GUA S02 3 0 ; [O]c1[c][c][c][c][c]1
4 SC5 1 GUA S03 4 0 ; [c]1[c][c][c][c][c]1
5 C5 1 GUA C01 5 0 ; [C][C]
6 P1 1 GUA P01 6 0 ; [C][N][C]
7 SC5 1 GUA S04 7 0 ; [c]1[c][c][c][c][c]1
8 SC5 1 GUA S05 8 0 ; [c]1[c][c][c][c][c]1
9 SC5 1 GUA S06 9 0 ; [c]1[c][c][c][c][c]1

[bonds]
; i j funct length force.c.
1 4 1 0.25 1250
3 5 1 0.25 1250
5 6 1 0.26 1250
5 7 1 0.25 1250

[constraints]
; i j funct length
2 3 1 0.24
2 4 1 0.24
3 4 1 0.24
7 8 1 0.24
7 9 1 0.24
8 9 1 0.24

[angles]
; i j k funct angle force.c.
1 4 2 2 62.0 25.0
1 4 3 2 122.1 25.0
2 3 5 2 136.9 25.0
3 5 6 2 131.1 25.0
3 5 7 2 114.6 25.0
4 3 5 2 143.4 25.0
5 7 8 2 122.0 25.0
5 7 9 2 62.0 25.0
6 5 7 2 97.7 25.0

[dihedrals]
; i j k l funct angle force.c.
1 2 3 4 2 -0.7 10.0
1 2 4 3 2 179.3 10.0
1 4 2 3 2 -179.3 10.0
1 4 3 2 2 0.7 10.0
2 1 4 3 2 -0.7 10.0
2 4 3 5 2 -129.5 10.0
2 4 3 9 2 -112.0 10.0
3 2 1 4 2 0.7 10.0
3 9 7 8 2 143.7 10.0
3 9 8 7 2 -63.8 10.0
4 2 3 5 2 137.7 10.0
4 2 3 9 2 81.9 10.0
5 7 8 9 2 0.1 10.0
5 7 9 8 2 -179.9 10.0
5 9 7 8 2 179.9 10.0
5 9 8 7 2 -0.1 10.0
7 5 9 8 2 0.1 10.0
8 7 5 9 2 -0.1 10.0

.mdp
integrator = md
dt = 0.01
nsteps = 500000
nstcomm = 10
comm-grps =

nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 10000
nstenergy = 1000
nstxtcout = 1000
xtc_precision = 100
xtc-grps =
energygrps =

nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.4

cutoff-scheme = group
coulombtype = Shift ;Reaction_field (for use with Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)

tcoupl = Berendsen
tc-grps = Solvent
tau_t = 1.0
ref_t = 300
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 1.0 1.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0

gen_vel = no
gen_temp = 320
gen_seed = 473529

constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30

Where Solvent : molecules + water

[tbereau]

Hi,

Couple tips:

1. The molecule you're parametrizing has many rings. The code tends to generate too many dihedrals (it's a bug), which can easily make the integration unstable. Have a quick look as to whether you could remove a few that may be redundant. Alternatively you can try with the alternative branch refactor

git checkout refactor

we may have fixed it there, so you would get fewer dihedrals. Alternatively, you can try to replace a few constraints by harmonic springs. I would then check that the molecule behaves by running a quick simulation in the gas phase (with sd integrator).

2. You may want to equilibrate your simulation first without pressure coupling. That would help ensure stability. Once you've equilibrated it a bit, turn on the pressure coupling.

Best,
Tristan

[krlitros87]

Hi Tristan,
First of all, thanks a lot for the quick reply!
I'll try with your tips ASAP and I'll let you know how this goes.
Thanks again,
Carlos

[ljmartin]

Hi @tbereau , could you please provide a link to this alternative branch? Does it still exist? I'm having a similar issue with dihedrals leading to LINCS failures
git checkout refactor returns: fatal: Not a git repository (or any of the parent directories): .git

[tbereau]

Hi @ljmartin, it's still here. My guess is you didn't checkout the git repo, but instead downloaded a zip of the code. Make sure to first git clone [XXX] with [XXX] being the address of the repo, and then checkout refactor.