Compressing a hydrogel with moving plates or sphere for mechanic moduli (using DPD)

I try to calculate the bulk modulus K and the elastic modulus E for spherical hydrogels using Dissipative Particle Dynamics. To calculate the bulk modulus I used a semipermeable spherical region with a repulsive harmonic potential (only repulsive for network beads, not for solvent). The radius of the sphere shrinks with every timestep (with a constant radial velocity). The calculation of the bulk modulus relies on the diagonal components of the stress tensor. The isotropic compression seems to work but if I try to generate a uniaxial compression (for the calculation of the elastic modulus) with two parallel plates, similar to the isotropic compression, no compression occurs. So I tried different versions using plane and box regions with the “move” parameter or compression by deformation of the simulation box - nothing worked.

The script for the isotropic compression:

units lj
atom_style angle

boundary p p p

angle_style cosine
bond_style harmonic

read_data /data/user/data.dat

group net type 1

variable dt equal 0.01

variable t equal step*${dt}

variable expandv equal 0.1

variable r0 equal 27.0

variable r equal ${r0}-${expandv}*v_t

region mySphere sphere 88.6 84.0 85.6 v_r side in
comm_modify vel yes

pair_style hybrid dpd 1 1.0 12345 srp 0.8 * min exclude yes
pair_coeff * * dpd 25.0 4.5 1.0
pair_coeff * * srp 100.0 0.5
pair_coeff 1 1 dpd 25.0 4.5 1.0
pair_coeff 2 2 dpd 25.0 4.5 1.0
pair_coeff 1 2 dpd 25 4.5 1.0
pair_coeff 1 3 none
pair_coeff 2 3 none
pair_coeff 3 3 srp 100.0 0.5

comm_modify cutoff 3.5

fix 1 all nve
fix wall net wall/region mySphere harmonic 25.0 0.0 2.5

compute mystress net stress/atom NULL
compute s net reduce sum c_mystress[1] c_mystress[2] c_mystress[3]

variable sigma equal -(c_s[1]+c_s[2]+c_s[3])

variable sigmapV equal -(c_s[1]+c_s[2]+c_s[3])/(4*v_r*v_r*v_r*PI)

thermo 100
thermo_style custom step temp epair emol press vol density

fix stress all print 1 "$t ${r} ${sigma} " file thermo_output.dat screen no

fix sigma all print 1 "${r} ${sigma} " file sigma.dat screen no

fix sigmapV all print 1 "${r} ${sigmapV} " file sigmapV.dat screen no

timestep 0.01
run 1000

Here are some snapshots:





Is there a way to compress the gelsphere uniaxial with two parallel plates with a script similar to the script for the isotropic compression with a sphere?

Hi @drOetker,

This is a tricky question. Since dpd potentials are very soft, it is very likely that, while compressing using only two planes, your gel can still strain and reorganize easily along the 2 orthogonal directions. Hence you are not really doing a uniaxial strain with regard to your gel.

Maybe you could consider constraining your gel particles in a cylindrical (or rectangular) region, which will maintain your cross-section surface constant along the compression, and then use the walls for compression of your gel. This would require a combined use of cylinder region (or block) and plane region. You could make a union of both and change the shape as desired as you are already doing with the sphere.

1 Like

It would be helpful to see an input that is not working as expected instead of one that does.

Thanks for your response Dr. Kohlmeyer. Since you are so familiar with all LAMMPS related topics, presenting you all the variations of combinations of (maybe topic related) input- script commands and parameters I could find in the command reference seemed like a waste of time for both of us.

  1. version, region plane:

units lj
atom_style angle

boundary p p p

angle_style cosine
bond_style harmonic

read_data /data/data

group net type 1

variable dt equal 0.1
variable t equal step*${dt}
variable dx equal 1.0
variable mdx equal -1.0
region upper_plane plane 1.0 85.0 85.0 1.0 0.0 0.0 move v_dx NULL NULL
region lower_plane plane 169.0 85.0 85.0 -1.0 0.0 0.0 move v_mdx NULL NULL

comm_modify vel yes

pair_style hybrid dpd 1 1.0 12345 srp 0.8 * min exclude yes
pair_coeff * * dpd 25.0 4.5 1.0
pair_coeff * * srp 100.0 0.5
pair_coeff 1 1 dpd 25.0 4.5 1.0
pair_coeff 2 2 dpd 25.0 4.5 1.0
pair_coeff 1 2 dpd 25 4.5 1.0
pair_coeff 1 3 none
pair_coeff 2 3 none
pair_coeff 3 3 srp 100.0 0.5

comm_modify cutoff 3.5

fix 1 all nve
fix upper_wall net wall/region upper_plane harmonic 25.0 0.0 2.5
fix lower_wall net wall/region lower_plane harmonic 25.0 0.0 2.5

compute mystress net stress/atom NULL
compute s net reduce sum c_mystress[1] #only xx component
variable sigmaxx equal -c_s

thermo 100
thermo_style custom step temp epair emol press vol density

fix stress all print 1 "$t ${sigmaxx} " file sigmaxx.dat screen no

timestep 0.01

run 1000

  1. version, region block:

units lj
atom_style angle

boundary p p p

angle_style cosine
bond_style harmonic

read_data /data/data

group net type 1

variable dt equal 0.01
variable t equal step*${dt}
variable velx equal 0.1

variable lowerx0 equal 1.0
variable upperx0 equal 168.0

variable lxlo equal ${lowerx0}+${velx}*v_t
variable lxhi equal ${lowerx0}+1.0+${velx}*v_t
variable lylo equal 1.0
variable lyhi equal 169.0
variable lzlo equal 1.0
variable lzhi equal 169.0

variable uxlo equal ${upperx0}+${velx}*v_t
variable uxhi equal ${upperx0}+1.0+${velx}*v_t
variable uylo equal 1.0
variable uyhi equal 169.0
variable uzlo equal 1.0
variable uzhi equal 169.0

variable dx equal 1.0

region lower_box block 1.0 2.0 1.0 169.0 1.0 169.0 side out move v_dx NULL NULL
region upper_box block 168.0 169.0 1.0 169.0 1.0 169.0 side out

comm_modify vel yes

pair_style hybrid dpd 1 1.0 12345 srp 0.8 * min exclude yes
pair_coeff * * dpd 25.0 4.5 1.0
pair_coeff * * srp 100.0 0.5
pair_coeff 1 1 dpd 25.0 4.5 1.0
pair_coeff 2 2 dpd 25.0 4.5 1.0
pair_coeff 1 2 dpd 20 4.5 1.0
pair_coeff 1 3 none
pair_coeff 2 3 none
pair_coeff 3 3 srp 100.0 0.5

comm_modify cutoff 3.5

fix 1 all nve

fix upper_box net wall/region upper_box harmonic 25.0 0.0 2.5
fix lower_box net wall/region lower_box harmonic 25.0 0.0 2.5

compute mystress net stress/atom NULL
compute s net reduce sum c_mystress[1] #only xx component
variable sigmaxx equal -c_s
thermo 100
thermo_style custom step temp epair emol press vol density

fix stress all print 1 "$t ${sigmaxx} " file sigmaxx.dat screen no

timestep 0.01

run 1000

  1. version, fix move:

units lj
atom_style angle

boundary p p p

angle_style cosine
bond_style harmonic

read_data /data/data

group net type 1
group gel type 1 2

variable dt equal 0.1
variable t equal step*${dt}
variable dx equal 0.01
variable mdx equal -0.01
region upper_plane plane 0.5 5.0 5.0 1.0 0.0 0.0
region lower_plane plane 169.0 5.0 5.0 -1.0 0.0 0.0

group uplaneg region upper_plane
group lplaneg region lower_plane

comm_modify vel yes

pair_style hybrid dpd 1 1.0 12345 srp 0.8 * min exclude yes
pair_coeff * * dpd 25.0 4.5 1.0
pair_coeff * * srp 100.0 0.5
pair_coeff 1 1 dpd 25.0 4.5 1.0
pair_coeff 2 2 dpd 25.0 4.5 1.0
pair_coeff 1 2 dpd 25 4.5 1.0
pair_coeff 1 3 none
pair_coeff 2 3 none
pair_coeff 3 3 srp 100.0 0.5

comm_modify cutoff 3.5

fix upper_wall net wall/region upper_plane harmonic 25.0 0.0 2.5
fix lower_wall net wall/region lower_plane harmonic 25.0 0.0 2.5
fix ufix uplaneg move linear ${dx} NULL NULL
fix lfix lplaneg move linear ${mdx} NULL NULL
fix nve gel nve

compute
.
.
.
run 1000

shouldn’t one of the planes be “side in” and the other “side out” for this to work?

here is a simple example using fix indent, which is about the easiest way to have a “traveling repulsive plane” in LAMMPS. For simplicity this is using a simple LJ droplet that is been prepared in the first part of the run and non-periodic boundaries.

units lj
lattice  sc 1.0

boundary m m m
region droplet sphere 0.0 0.0 0.0 10.0
create_box 1 droplet
create_atoms 1 region droplet

pair_style lj/cut 4.0
pair_coeff * *  1.0 1.0
mass * 1.0
fix 1 all nve

run 100 post no

minimize 0.0 0.0 1000 10000

reset_timestep 0

thermo 500

run 2500

variable umove equal vdisplace(10.0,-0.1)
variable lmove equal vdisplace(-10.0,0.1)
fix 2 all indent 1.0 plane y v_umove hi
fix 3 all indent 1.0 plane y v_lmove lo

dump 1 all movie 50 movie.mp4 type type size 1024 960 # ssao yes 3561234 0.6

run 20000

1 Like

I tried to adapt the fix indent example:

units lj
atom_style angle

boundary p p p
angle_style cosine
bond_style harmonic

read_data /data/data

variable umove equal vdisplace(10.0,-0.01) # the simulationbox in the data file is defined 0.0 to 20.0 for each dimension
variable lmove equal vdisplace(-10.0,0.01)

comm_modify vel yes

pair_style hybrid dpd 1 1.0 12345 srp 0.8 * min exclude yes
pair_coeff * * dpd 25.0 4.5 1.0
pair_coeff * * srp 100.0 0.5
pair_coeff 1 1 dpd 25.0 4.5 1.0
pair_coeff 2 2 dpd 25.0 4.5 1.0
pair_coeff 1 2 dpd 25 4.5 1.0
pair_coeff 1 3 none
pair_coeff 2 3 none
pair_coeff 3 3 srp 100.0 0.5

comm_modify cutoff 3.5

fix 1 all nve
fix 2 all indent 5.0 plane x v_umove hi
fix 3 all indent 5.0 plane x v_lmove lo

compute …

.
.
.
run …

changing lo and hi with each other:

blue particles are unbound solvent, red is the gelnetwork (bonded), yellow represnts the particles, generated by the srp potential.

I don’t think you want to apply the indenter to all particles. In my demo was no explicit solvent.

If your box is positioned from 0.0 to 20, then you have to adapt these variables accordingly, so that the indenter planes are placed properly. You need to study the manual.

Another item worth mentioning from viewing your movies is that you need to take finite size effects into account. Your “target object” is rather small, but curvature and surface tension will have a larger impact on your results the smaller the object is, so you will see that your results will change when using larger objects.