Tensile stress calculation in LAMMPS

Hello everyone,

I want to do the tensile simulation of tobermorite crystal in LAMMPS. I need clarification on the stress calculation. I used compute stress/atom and compute reduce sum to calculate the stress. However, from the doc we know that the stress is in P*V form. So, (1) should we divide the stress by the volume of the simulation box or by volume of the actual crystal? Note that, I have some vacuum in the box.

(2) I calculated the hydrostatic pressure from the stress/atom calculation and that totally matches with the pressure that LAMMPS is calculating.

But somehow, I am not getting correct stress-strain relationship curve. It seems my structure is not experiencing fracture at all. I am using morse potential for fracturing the material. Note that I am also using CVFF-CLAYFF potential with some non-bonded parameters from CSHFF potential developed by Shahsavari et al. I used fix deform to perform the tensile stress while keeping 0 pressure in Y and Z direction. I have seen that the box volume was increasing during simulation.

I need to know if my stress calculation understanding is correct or not. Can anyone help me please?

Here is my input script

Tobermorite-14A

Uniaxial strain in x-direction

Creating and observing dislocation

Ca/Si ratio= 120/144 = 0.833

#---------Define Input Variables------------------#
variable timestep equal 1
variable Nsamplefreq equal 10
variable Nsamplesize equal 100
variable Nblocksize equal v_Nsamplefreqv_Nsamplesize
variable Nblock equal 70
variable Nsteps equal v_Nblocksize
v_Nblock

Initialization

echo screen
log debug_14t.log
newton on

units real
atom_style full
boundary p p p

Bonded Interaction

bond_style morse
angle_style harmonic
dihedral_style harmonic
improper_style cvff

Non-Bonded Interaction

pair_style lj/cut/coul/long 10.0 10.0 # cutoff1 = 2.5*sigma
read_data data.14tcshproject

kspace_style ewald 1e-6

pair_coeff 1 2 1.84e-6 3.7064
pair_coeff 1 3 0.00146 4.898
pair_coeff 1 4 0.00146 4.898
pair_coeff 1 5 0.00146 4.898
pair_coeff 1 6 0.00603 5.0168
pair_coeff 1 7 0.0 3.12
pair_coeff 1 8 0.0 3.12

pair_coeff 2 3 0.000563 3.6716
pair_coeff 2 4 0.000563 3.6716
pair_coeff 2 5 0.000595 3.6627
pair_coeff 2 6 0.000530 3.6298
pair_coeff 2 7 0.0 1.8532
pair_coeff 2 8 0.0 1.8532

pair_coeff 3 4 1.2433 3.0687
pair_coeff 3 5 0.0455 4.0654
pair_coeff 3 6 0.00525 4.7557
pair_coeff 3 7 0.0 1.7766
pair_coeff 3 8 0.0 1.7766

pair_coeff 4 5 0.0455 4.0654
pair_coeff 4 6 0.0455 4.0654
pair_coeff 4 7 0.0 1.7766
pair_coeff 4 8 0.0 1.7766

pair_coeff 5 6 0.8705 3.2513
pair_coeff 5 7 0.0 1.7766
pair_coeff 5 8 0.0 1.7766

pair_coeff 6 7 0.0 0.0
pair_coeff 6 8 0.0 0.0
pair_coeff 7 8 0.0 0.0

neighbor 3.0 bin
neigh_modify delay 0 every 1 one 10000 check yes
comm_modify mode single cutoff 3.0
#pair_modify mix arithmetic tail no

#------------ define groups ------------------xxx
group cao type 1
group st type 2
group ob type 3
group obts type 4
group oh type 5
group o* type 6
group h* type 7
group ho type 8
group water type 6 7
group solid subtract all water

#set type 1 charge 1.7
#set type 2 charge 1.72
#set type 3 charge -1.14
#set type 4 charge -1.14
#set type 5 charge 1.0
#set type 6 charge -0.82
#set type 7 charge 0.41
#set type 8 charge 0.29

region leftbdy block INF -18.13 INF INF INF INF units box
group leftbdy region leftbdy
region rightbdy block 6.24 INF INF INF INF INF units box
group rightbdy region rightbdy
group mobile subtract all leftbdy rightbdy

compute newtemp all temp
velocity all create 300 1234 temp newtemp dist gaussian mom yes rot no

minimize 1.0e-4 1.0e-6 100 1000
min_style cg
min_modify dmax 0.1

fix 1 mobile nve
fix 2 all langevin 300 300 100 1234
thermo 100
run 10000
unfix 1
fix 3 mobile npt temp 300 300 100 iso 0 0 1000 drag 1
run 10000

unfix 3

compute peratom all stress/atom newtemp
compute stp all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]

variable sxx equal c_stp[1]
variable syy equal c_stp[2]
variable szz equal c_stp[3]
variable p equal -(v_sxx+v_syy+v_szz)/(3vol)
variable stress equal v_sxx
0.000101325/vol

fix 4 all deform 1 x erate 1e-5 units box remap x
fix 5 mobile npt temp 300 300 100 y 0 0 1000 z 0 0 1000 drag 1

fix 6 all ave/time {Nsamplefreq} {Nsamplesize} {Nblocksize} v_stress ave window {Nblock} file stress.profile

thermo 100
thermo_style custom step temp v_stress vol press v_p
run 70000

Regards,

Baig

The divide by volume is just a bookkeeping
normalization. If you don’t have a periodic
filll-the-box system, it’s up to you
what to divide by. LAMMPS doesn’t have a
better idea than you, what that should be.