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.