BD stops after some steps

User discussions
1

GROMACS version: 2020.3
GROMACS modification: No
I am performing Brownian Dynamics with hard spheres. But the simulation stops after a random number of steps, Gromacs keeps running but I don’t have evolution on the number of steps in any output file.
First I was using only the C12 part of the LJ potential and I thought this problem could be due to two particles getting to close together and the potential getting to big. So I added the C6 term, the problem took longer to appear, the simulation was able to run for a longer time, but it is still there.
I am putting below the mdp file I am using, the itp for the particles and the topology.
The system is 1000 particles in a box of 5 x 5 x 5 um3.
Any help is greatly appreciated! Thank you so much!

mdp file_
ld-seed = 55560 ; Use random seed
;################################# INTEGRATOR ##################################
integrator = bd ; Langevin thermostat
dt = 1000 ; Timestep (ps)
nsteps = 20000000 ; Simulation duration (timesteps)
nstcomm = 0 ; Center of mass motion removal interval
comm-mode = None ; Center of mass motion removal mode
bd_fric = 24000 ; The friction coefficient for each particle is calculated as mass/tau-t.
;################################## ENSEMBLE ###################################
ref_t = 295 ; System temperature (K)
tau_t = 2.0 ; Thermostat time constant (ps)
tc_grps = system ; Apply thermostat to complete system
;############################## IMPLICIT SOLVENT ###############################
implicit_solvent = no ; NO implicit solvent
;########################### NONBONDED INTERACTIONS ############################
cutoff_scheme = Verlet ; Method of managing neighbor lists
pbc = xyz ; Periodic Boundary Conditions in all directions
coulombtype = cut-off ; Calculate coulomb interactions using PME
rcoulomb = 10 ; Coulomb cutoff
vdwtype = cut-off ; Calculate van der Waals interactions using PME
rvdw = 10 ; Van der Waals cutoff
rlist = 10 ; Neighbor list cutoff
nstlist = 100 ;
;################################### OUTPUT ####################################
nstlog = 5000 ; Log output interval (timesteps)
nstenergy = 5000 ; Energy output interval (timesteps)
nstcalcenergy = 5000 ; Energy calculation interval (timesteps)
nstxout = 1000 ; Trajectory output interval (timesteps)
nstvout = 5000 ; Velocity outout interval (timesteps)
nstfout = 5000 ; Force output interval (timesteps)
nstxout-compressed = 0 ;

itp file_________
; Topology of solid particles.
; No include forcefield parameters

[ defaults ]
; nb-type combination-rule gen-pairs fudgeLJ fudgeQQ
1 1 no

[ atomtypes ]
;type mass charge ptype c6 c12
FA 4.8e+04 0.000 A 3.898e+04 6.090e+08
FB 3.0e+08 0.000 A 2.494e+12 2.494e+24
FC 1.03e+3 0.000 A 1.818e+03 1.326e+06

[nonbond_params]
FA FA 1 3.898e+04 6.090e+08
FB FB 1 2.494e+12 2.495e+24
FC FC 1 1.818e+03 1.326e+06

[ moleculetype ]
;name nrexcl
F005 1 ; sigma=05nm Epsilon= 0.6236 kJmol-1 mass= 48000 Da C12= 6.08984e+08 kJmol-1nm12

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 FA 1 F005 FA 0 0.00 4.8e+04

[ moleculetype ]
;name nrexcl
F100 1 ; sigma=100nm Epsilon= 0.6236 kJmol-1 mass= 3.0e+08 Da C12= 2.49440e+24 kJmol-1nm12

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 FB 1 F100 FB 0 0.00 3.0e+08

[ moleculetype ]
;name nrexcl
F003 1 ; sigma=03nm Epsilon= 0.6236 kJmol-1 mass= 1029 Da C12= 1.32563e+06 kJmol-1nm12

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 FC 1 F003 FC 0 0.00 1.02914e+03

top file___
; Topology of 1000 particles.

#include “solidParticlesC6.itp”

[ system ]
; systemname
1000 particles 005nm

[ molecules ]
; molecule number
F005 1000

2

I forgot to say I generated the initial configuration with:

gmx insert-molecules -ci 1particles5nm.gro -nmol 1000 -box 5000 5000 5000 -o 1000particles5nm.gro -radius 2.5

Thank again.

3

It is a bit difficult to judge what the values of the parameters should be with your different than usual sizes and masses, but I think the timestep is many orders of magnitude too large. My guess is that atoms move thousands of box length and pbc corrections take forever.

4

I get what you are saying but with BD ‘the system is assumed to be over-damped, large timesteps can be used’ as state in the Gromacs manual. So 1ns is too large?

5

What is relevant for BD is dt/gamma, so one can’t speak about dt in isolation. But I would guess that 1 ns is several orders of magnitude too large. Try and see.

1 Like
6

Thank you so much for your clarifications. I am learning a lot.

I have dt/gamma around 1e-1.
What would be reasonable values so I can adjust dt??

7

If you read the manual, you see that the displacement in a step is
force*dt/gamma plus a random term. So you need to come up with a
reasonable displacement and maximum force so you can estimate reasonable
parameter values. The proper check is to see that you distribution of
states is not affected by the time step.

1 Like
8

Thank you! I will do that.