Minus stress during deformation

LAMMPS LAMMPS General Discussion
1

Dear lammps users,

I am currently attempting to obtain the tensile strength (UTS) and Young’s modulus of my system by analyzing the stress-strain curve. However, I have noticed that the stress calculated at low strain is negative, which is unexpected since the stress should be zero at zero strain. Could you please advise me on what may be causing this issue and how I can obtain a proper stress-strain curve with zero stress and zero strain?

image
image930×639 25.5 KB

To provide some context, my system consists of both polymer (PCFF) and MOF (UFF), and I have equilibrated it in accordance with a reference paper. I then applied deformation in the z direction and measured the pressure using pzz. I have included my Lammps input below.

After analyzing my system, I have observed that the pzz value is consistently high after equilibration (around 1000 atm) even at NPT ensemble (fix 1 all npt temp 300 300 100 iso 1 1 1000). Could this be the cause of the issue? If so, what does a high pzz value indicate and how can I resolve this issue?

image
image1523×671 52.4 KB

Thank you for your assistance.

include "system.in.init"
read_data "system.data"
include "system.in.settings"

neigh_modify every 1 delay 0 check yes
timestep 1.0
 
dump my_dump all atom 1000 dump.lammpstrj
minimize 1.0e-7 1.0e-9 1000 10000

velocity all create 600 1231
thermo_style   custom step temp density press etotal lx ly lz
thermo 1000

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 5000
unfix 1
fix 1 all npt temp 300 300 100 z 19.7 19.7 1000
run 5000
unfix 1

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
fix 1 all npt temp 300 300 100 z 592 592 1000
run 5000
unfix 1

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
fix 1 all npt temp 300 300 100 z 986 986 1000
run 5000
unfix 1

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
fix 1 all npt temp 300 300 100 z 493 493 1000
run 5000
unfix 1

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
fix 1 all npt temp 300 300 100 z 98 98 1000
run 5000
unfix 1

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
fix 1 all npt temp 300 300 100 z 9.8 9.8 1000
run 5000
unfix 1

fix 1 all nvt temp 600 600 100 
run 5000
unfix 1
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
fix 1 all npt temp 300 300 100 z 0.98 0.98 1000
run 800000
unfix 1

fix 1 all nvt temp 300 300 100
run 1000000
unfix 1

variable tmp equal "lz"
variable L0 equal ${tmp}
variable strain equal "(lz - v_L0)/v_L0"
variable p1 equal "v_strain"
variable p2 equal "-pxx/10000*1.01325"
variable p3 equal "-pyy/10000*1.01325"
variable p4 equal "-pzz/10000*1.01325"
variable p5 equal "lx"
variable p6 equal "ly"
variable p7 equal "lz"
fix        1 all npt temp 300 300 50 x 0 0 1000 y 0 0 1000 drag 2
fix        2 all deform 1 z erate 1e-6 units box remap x
fix def2 all print 100 "${p1} ${p4}" file def2.txt screen no
timestep   1
reset_timestep 0
run        1000000
write_data deform.data
unfix 1
unfix 2

---

 system.data
LAMMPS Description

2440  atoms
2780  bonds
6312  angles
5820  dihedrals
2704  impropers

9  atom types
46  bond types
892  angle types
9  dihedral types
4  improper types

0 29.34 xlo xhi
0 25.409 ylo yhi
0.0 120 zlo zhi
-14.67 0 0 xy xz yz
---
 system.in.init
    atom_style full
    units           real
    pair_style      hybrid lj/cut/coul/long 10 lj/class2/coul/long 10
    bond_style      hybrid harmonic class2
    angle_style     hybrid fourier cosine/periodic class2
    dihedral_style  hybrid harmonic class2
    improper_style  hybrid fourier class2
    special_bonds   lj/coul 0.0 0.0 1.0
    pair_modify     tail yes mix arithmetic
    kspace_style    ewald 0.0001
  boundary p p p
2

It’s not really a LAMMPS question, but an MD one.

If you have a system equilibrated in some finite temperature, you must not use minimize as it will destroy the equilibration. Similarly, assigning new temperature with velocity command throws the system out of the equilibrium. Finally, if your system is anisotropic, using fix npt iso may lead to substantial stress in some directions, because with iso keyword all dimensions are scaled similarly.

There may also be other issues with your script. I advise you to find someone local who you can consult on this subject.

1 Like
3

During your equilibration you slowly decrease the pressure in only the z direction but for the final deformation you set the pressure in the x and y directions to 0. You should properly equilibrate all three diagonal components of the pressure virial before starting a production run. And as @mkanski says, these should be allowed to equilibrate independently with, e.g., npt aniso.

1 Like
4

Thank you for your help! I got more reasonable outcomes with fix npt aniso.

However, I have a question regarding the minimize and velocity commands. Typically, I optimize the initial structure of my system using the minimize command before starting a simulation because I believe it requires a geometric optimization process. In this case, would you still don’t recommend to use minimize option as well?
Additionally, I’ve heard that the velocity can take a long time to assign the correct velocity for a specific temperature, so it may be reasonable to assign velocities initially with a Gaussian distribution. Do you think it also makes sense?

Is it correct to understand that it is recommended not to use minimize and velocity in systems where equilibration is important?

5

Thank you for your response. I tried to control pressure in z direction because I followed the equilibration steps applied in the paper. However, npt aniso worked for me!