Hello Lammps Developers/Users,
I want to compute the (short- & long-term) self-diffusion coefficients of different colloidal dispersions with an implicit solvent. I have already written an input script according to Bolintineanu, Grest, Lechman et al. 2014 and the example from the documentation "DIFFUSE" (see code below) and checked it with some mailing entries like https://sourceforge.net/p/lammps/mailman/message/34081599/ ("Mysterious attraction between hard spheres...").
Now, when I run a test for a diluted dispersion (e.g. volume fraction 0.06 %), I do not obtain an approximate value for the diffusion coefficient in infinite dilution (given by the Stokes-Einstein relation D_0 = k_B T /_6πμ_0a). Even if I change the radius of the colloid particles, the obtained coefficients do not scale accordingly.
I'm running the pre-built version of LAMMPS for Ubuntu.
Here the input script:
////////////////////////////////////////////////////////////////////////////////////
variable d1 equal 20.0 # 2*a -- diameter
variable rc equal 30.0 # 3*a -- cut-off radius
variable irc equal 20.002 # 2.0002*a -- inner cut-off
variable eta0 equal 1.0 # viscosity
variable Tempt equal 300.0 # target temperature
# ~ variable shearR equal 1e-12 #shear rate only for LubricateU (non-inertial)
variable AVdW equal 0.1 # Hamaker
variable Sigmalj equal 0.25 # solvent/constituent size
units nano
dimension 3
atom_style sphere
boundary p p p
newton on
lattice fcc 300.0
region box block 0 5 0 5 0 5
create_box 1 box
create_atoms 1 box
comm_modify vel yes # ghost atoms needed for hydrodynamics interactions
pair_style hybrid/overlay colloid \{rc\} brownian {eta0} 1 1 \{irc\} {rc} \{Tempt\} 5878567 1 1 lubricate {eta0} 1 1 \{irc\} {rc} 1 1
# I have also tried the following combinations:
# ~ pair_style hybrid/overlay brownian \{eta0\} 1 1 {irc} \{rc\} {Tempt} 5878567 1 1 lubricate \{eta0\} 1 1 {irc} \{rc\} 1 1 \# no colloid pair style: same results since dilute
\# \~ pair\_style hybrid/overlay brownian {eta0} 1 1 \{irc\} {rc} \{Tempt\} 5878567 1 1 lubricateU {eta0} 1 \{irc\} {rc} ${shearR} 1 1 # this non-intertial version was unstable
pair_coeff * * colloid \{AVdW\} {Sigmalj} \{d1\} {d1} ${rc}
pair_coeff * * brownian
pair_coeff * * lubricate
velocity all create ${Tempt} 4928459 dist gaussian
# ~ fix 1 all nve/noforce
fix 1 all nve/sphere
timestep 0.001
compute msd all msd com yes
variable twopoint equal c_msd[4]/6/(step*dt+1.0e-6)
fix 9 all vector 10000 c_msd[4]
variable fitslope equal slope(f_9)/6/(10000*dt)
thermo_style custom step temp c_msd[4] v_twopoint v_fitslope
thermo 10000
run 100000000
Am I missing something?
Thanks in advance.
Best,
Hector Rusinque