Brownian Dynamics Computer Simulations of Quenched Lennard-Jones- like Fluids

LAMMPS LAMMPS Mailing List Mirror
1

Hi everyone,

I’m trying to perform brownian dynamics in lammps of quenched Lennard-Jones- like fluids. I read in the documentation page that this is done by using fix NVE and fix langevin together. I’d like to see gelification of this fluid.

I’m trying to reproduce the simulations performed in http://www.tandfonline.com/doi/abs/10.1080/08927029908022123

My system consisted of 864 spherical particles of diameter sigma (sigma = 1) in an implicit fluid. I chose the traditional time unit for colloidal liquids (a²/D), where D is the self diffusion coefficient. I simulated 82 a²/D , in colloidal time units. Being D = KbT/ coeff_friction, I set coeff_friction = 1, So D = T.

I equilibrated my system at T = 2.0 (homogeneous liquid) so I quenched to 0.7 and ran the simulation. I used the 12-6 Lennard-Jones potencial for pairwise interactions.

By changing the timestep I observe that the system evolves in different ways and this is not what I expect to happen.

If I chose very low timesteps the system seems to reach equilibrium in the earlier timesteps and I dont see any gelification.

My question is : How can I chose the correct timestep of this type of simulation?

I send you my code and graphics I obtained for RDF using different timesteps.

Thank you so much for your attention.

units lj
dimension 3
atom_style atomic
neigh_modify delay 0 every 1

#settings

variable N equal 864
variable T equal 0.7

variable phi equal 0.16

variable cuttof equal 2.5
variable epsilon equal 1
variable sigma equal 1
variable a equal {sigma}/2.0 variable cuttof equal 2.5*{sigma}

variable rho equal 6*{phi}/PI variable L equal ({N}/${rho})^(1.0/3.0)

variable D0 equal T variable coeff_fric equal 1.0 variable dt_equilibration equal 8e-4*{a}{a}/{D0}
variable damp_equi equal 100${dt_equilibration}

problem setup

region simbox block 0 $L 0 $L 0 $L
create_box 1 simbox
create_atoms 1 random $N 10 simbox

variable mass_par equal 1.0
mass * ${mass_par}
velocity all create 2.0 10 dist gaussian

pair_style mie/cut {cuttof} pair_coeff * * {epsilon} ${sigma} 12 6

minimize 1.0e-3 1.0e-3 1000 10000
reset_timestep 0

equilibration

fix 1 all nvt temp 2.0 2.0 ${damp_equi}

variable equil_steps equal 400000
timestep {dt_equilibration} thermo 1000 run {equil_steps}
reset_timestep 0

#Production
unfix 1
velocity all scale ${T}

variable damp equal {mass_par}/{coeff_fric}
variable dt equal 1.0e-5*{a}*{a}/${D0}

variable total_steps equal round(81.9*{a}*{a}/{D0}/{dt})
variable dump_interval equal round(3.2*{a}*{a}/{D0}/{dt})

fix NVE all nve
fix 1 all langevin {T} {T} ${damp} 48279 zero yes

dump 1 all custom ${dump_interval} dump_struc.lammpstrj id type x y z

thermo 1000
run ${total_steps}

RDf.pdf (53 KB)

2

Hi Fellipe,

Good job in performing a very thorough analysis. I think there is just
one oversight: You lower your time step size, but you do not increase
the number of time steps you perform. For your very low time step
size, this means that all simulations, no matter what the actual
damping time is, have barely changed from the minimised state at all.
Try to couple the number of equilibration steps to the time step size
so that each simulation run corresponds to an equal physical
simulation time. Then you should obtain results qualitatively more
similar to those of the larger time step size (and get some idea as to
how sensitive your results are to your time step size).

3

Hi Stefan,

First of all thank you for your attention to my problem.
I’m so sorry I didn’t understand why the number of time steps is not increasing when I decrease the time step size.

For the equilibration of all simulations I used dt_equilibration = 8e-4a²/D (a and D are constants) and I performed the same number of equilibration time steps = 400000 to attain T = 2.0 (this is my initial point, from this point the system initiates the gelation process) so I imagine that all simulations start from equivalent configurations.

For the production (this is the part of the simulation from which the RDF’s were calculated) I used dt = 0.2e-3 a²/D, dt = 1.0e-3 a²/D and dt = 1e-5 a²/D. In the following part of the code I set time step size and total number of time steps

variable dt equal 1.0e-5*{a}*{a}/{D0} # or dt equal 1.0e-3*{a}{a}/{D0} or dt equal 0.2e-3{a}*{a}/${D0}

variable total_steps equal round(81.9*{a}*{a}/{D0}/{dt})

In time units, all simulations have the same length = 81.9*{a}*{a}/{D0}. But for each one I have different number of time steps ( total_steps = 81.9*{a}*{a}/{D0}/${dt} ) because dt changes.

thank you so much!

4

Hi Stefan ,

I just realized I forgot to change dt in the production part.
thank you!

5

Yeah, that's what I meant. Your input is a bit complex so I don't know
if that is where the problem is, but it is easy to check. Good luck!