Shearing microgels

Hello everyone,

I am trying to shear my microgel and I would like to calculate its gyration tensor components along the way as well as the intrinsic viscosity and I can’t seem to get it to work:

LAMMPS script for polymer gel under shear flow

Box and units (LJ units and periodic boundaries)

units lj
atom_style bond
boundary p p p

Neighbor list settings

neighbor 0.3 bin
neigh_modify every 1 delay 1 check yes

RESTART file generation

restart 100000 mgel.restart.a mgel.restart.b

READ initial configuration

read_data data_final_alfa0
reset_timestep 0

Convert box to triclinic for shear deformation

change_box all triclinic

Define groups

group all type 1

Dump configurations

dump mgelconf all custom 50000 mgel.dump.* id type x y z vx vy vz ix iy iz
dump_modify mgelconf format line “%d d .10f .10f .10f .10f .10f %.10f %d %d %d”
dump_modify mgelconf pbc yes
dump_modify mgelconf sort id

Pair interactions

pair_style lj/cut 1.122462048309
pair_modify shift yes
pair_coeff * * 1.0 1.0

Bond interactions (FENE)

bond_style fene
special_bonds fene
bond_coeff 1 30 1.5 1 1

Shear flow setup using SLLOD algorithm

Shear rate = 0.001 in LJ units

variable shear_rate equal 0.001

Fixes for SLLOD integration

First: deform the box at constant shear rate

fix deform all deform 1 xy erate ${shear_rate} remap v

Second: Use SLLOD equations with NVT thermostat

The “temp” thermostat acts on thermal velocities relative to shear flow

fix sllod all nvt/sllod temp 1.0 1.0 0.2 tchain 30

Compute for stress tensor components

compute stress all stress/atom NULL virial
compute p all pressure NULL temp c_stress
compute p all pressure temp c_stress
compute p all pressure NULL

Compute gyration tensor and related properties

compute gyrtensor all gyration tensor
compute rg all gyration

Extract individual components of gyration tensor

variable Rg2 equal c_rg
variable Rgxx equal c_gyrtensor[1]
variable Rgyy equal c_gyrtensor[2]
variable Rgzz equal c_gyrtensor[3]
variable Rgxy equal c_gyrtensor[4]

Compute intrinsic viscosity

σ_xy = -1/V * ∑(R_ix * F_iy)

Intrinsic viscosity = σ_xy / (shear_rate)

variable box_vol equal vol
variable sigma_xy equal -(c_p[4]) # c_p[4] is the xy component of stress tensor
variable intrinsic_viscosity equal v_sigma_xy/${shear_rate}

Additional computes for monitoring

compute ke all ke
compute pe all pe
compute temp_all all temp

Density profile calculation

compute densprof all chunk/atom bin/sphere 0 0 0 0.0 120.0 120
fix fordensprof all ave/chunk 2 20 10000 densprof density/number file densprofmgel

Recentering fix

fix recentermgel all recenter 0 0 0 shift all units box

Thermodynamic output with enhanced observables

thermo 1000
thermo_style custom step temp etotal pe ke press c_rg v_Rgxx v_Rgyy v_Rgzz v_Rgxy v_sigma_xy v_intrinsic_viscosity

Timestep

timestep 0.002

Minimization

minimize 1.0e-4 1.0e-6 1000 1000

Production run

run 50000000

Write final configuration

write_data data_final_sheared

I am happy for any feedback, thank you in advance!!

Leah

PS.: I also later on would like to shear a bulk of tighly packed microgels, do I simply remove the recentering command and lammps takes proper care of it (Lees-Edwards-boundary conditions, …) or is it not possible this way. Right now with one microgel I can recenter it and no boundaries are in the way but I am aware that this is not the case in bulk simulations.

Please have a good luck at this post and follow the suggestions: Please Read This First: Guidelines and Suggestions for posting LAMMPS questions

Otherwise you are not likely to get much help since we won’t be able to understand what your problem is and reproduce it to make suggestions.

Hello,

Thank you for the fast reply, I apologise for I was too inconcrete…
I got the problems with the variables solved and now the code is running… I shear like this:

Shear flow setup using SLLOD algorithm

Shear rate = 0.001 in LJ units

variable shear_rate equal 0.001

Fixes for SLLOD integration

First: deform the box at constant shear rate

fix deform all deform 1 xy erate ${shear_rate} remap v

Second: Use SLLOD equations with NVT thermostat

The “temp” thermostat acts on thermal velocities relative to shear flow

fix sllod all nvt/sllod temp 1.0 1.0 0.2 tchain 30

################################################
#which now runs but I also hopes will produce useful results
#AND then write output the quantities I am interested like so
################################################

Compute for stress tensor components

compute stress all stress/atom NULL virial
compute p all pressure NULL temp c_stress

Sum the xy components over all atoms

compute p all reduce sum c_stress[4]

Convert to macroscopic stress

variable sigma_xy equal -c_p/(vol)

Compute intrinsic viscosity

σ_xy = -1/V * ∑(R_ix * F_iy)

Intrinsic viscosity = σ_xy / (shear_rate)

variable box_vol equal vol
#variable sigma_xy equal -(c_p[4]) # c_p[4] is the xy component of stress tensor
variable intrinsic_viscosity equal v_sigma_xy/${shear_rate}

Additional computes for monitoring

compute ke all ke
compute pe all pe
compute temp_all all temp

Recentering fix

fix recentermgel all recenter 0 0 0 shift all units box

Thermodynamic output with enhanced observables

thermo 1000
thermo_style custom step temp etotal pe ke press v_sigma_xy v_intrinsic_viscosity

… (any feedback appreciated)

Now it also would be cool if you could tell me: If I move to a bigger system (here I have one microgel kept CENTERED in a very large box with no contact to any boundaries) with many microgels tightly packed, is lammps then able to account for the needed boundary conditions? And would I be therefore getting physical simulations and also data without any kinks after stopping and restarting using the restart file.

Best,
Leah

My previous statement still applies.