Error while undergoing uni-axial tensile deformation of a network of chains

Hello,

I am running the LAMMPS Dec 2018 version and studying a random network of chains undergoing tensile deformation. The length of each chain is ~ 600 Angstrom (~30 beads for one chain) and there are 4900 such chains in random orientation within a simulation box of size (67000 Ang by 67000 Ang by 1000 Ang). The data file was created in such a way that the distance between each chain will never be less than 110 Angstrom, which is basically a distance within the potential well of the force potential (lj 9/6) between the chains. Each chain is straight and has bonds, angles and dihedrals.

The initial energy minimization step looks good, in my opinion and converges near to zero for all x, y and z directions (log file posted below).

After that, I am trying to deform in the x direction by maintaining consistent pressure (~0) in the y and z direction to reproduce the Poisson’s effect of the lateral strain and make it close to the real scenario.

However, when I am trying to achieve this using:

"fix 9 all npt temp 300 300 100 y 1.0 1.0 500.0 z 1.0 1.0 500.0 drag 2.0

fix 10 all deform 1 x erate 0.00001 units box remap x"

I see the following error of "Non-numeric pressure - simulation unstable "

I tried several ways to address this as I know this error is a common occurrence in LAMMPS. Here are the following ways I tried:

(1) Increase the Pdamp for y and z to 500 ( more than ~1000 steps as timestep=0.5) and also to 1500. But I got the same error.

(2) Instead of NPT, I tried using the following commands along with “fix deform”

fix 8 all nvt temp 300 300 100
fix 9 all press/berendsen y 0.0 0.0 500.0 z 0.0 0.0 500.0

For this case, the simulation is running but with “nan” values of pxx, pyy and pzz.

(3) Reducing the erate and the time step (from 0.5 to 0.25) also are showing the same error for the approaches detailed in (1) and (2)

However, if I am using “fix nve” and “fix deform” together, the simulations run well. But I cannot capture the effect of lateral strain there as the box volume is fixed.

I am trying hard to understand the physics why I am failing to control the pressure. So if anyone can please provide me any suggestion, that will be extremely helpful.Thank you so much for your time and help.

The input files and log files of one sample run are posted below.

INPUT FILE

dimension 3
boundary p p p

units real
atom_style full
newton on

echo both

read_data data.random_0.6.TwoFibers

neighbor 100.0 nsq
neigh_modify every 10 delay 10 check yes page 100000

##################################SPECIFYING THE POTENTIAL
##NON BONDED POTENTIALS

pair_style hybrid lj/cut 500.0 lj96/cut 1500.0
pair_modify shift yes

##BETWEEN BEADS IN A SINGLE CHAIN

pair_coeff 1 1 lj/cut 1004.8715 201.66
pair_coeff 2 2 lj/cut 1004.8715 201.66
pair_coeff 3 3 lj/cut 1004.8715 201.66

.
.
.
.
.
pair_coeff 4895 4895 lj/cut 1004.8715 201.66
pair_coeff 4896 4896 lj/cut 1004.8715 201.66
pair_coeff 4897 4897 lj/cut 1004.8715 201.66
pair_coeff 4898 4898 lj/cut 1004.8715 201.66
pair_coeff 4899 4899 lj/cut 1004.8715 201.66
pair_coeff 4900 4900 lj/cut 1004.8715 201.66

#INTER POTENTIAL

pair_coeff 1 24900 lj96/cut 2128.4206 109.83
pair_coeff 2 3
4900 lj96/cut 2128.4206 109.83
pair_coeff 3 44900 lj96/cut 2128.4206 109.83
pair_coeff 4 5
4900 lj96/cut 2128.4206 109.83

.
.
.
.
pair_coeff 4895 48964900 lj96/cut 2128.4206 109.83
pair_coeff 4896 4897
4900 lj96/cut 2128.4206 109.83
pair_coeff 4897 48984900 lj96/cut 2128.4206 109.83
pair_coeff 4898 4899
4900 lj96/cut 2128.4206 109.83
pair_coeff 4899 4900 lj96/cut 2128.4206 109.83

#################BONDED POTENTIALS

bond_style harmonic
bond_coeff 1 629.72 203.04

angle_style harmonic
angle_coeff 1 3156851634.53 163.272

dihedral_style harmonic
dihedral_coeff 1 114307.2 1 1

special_bonds lj/coul 1.0 1.0 1.0

##################OUTPUT

thermo_style custom step temp press vol pe ke ebond eangle edihed etotal pxx pyy pzz lx ly lz
thermo 500
thermo_modify flush yes
timestep 0.5

##################################################DUMP

dump 1 all custom 1000000 Random.lammpstrj id type x y z
dump_modify 1 sort 1

restart 5000000 poly.restart

########################MINIMIZATION

displace_atoms all random 0.01 0.01 0.01 443322

min_style cg
min_modify line quadratic
minimize 1.0e-20 1.0e-15 100000000 1000000000

write_restart equil.restart2

You are getting this warning:
WARNING: Bond/angle/dihedral extent > half of periodic box length (…/domain.cpp:903)

Do you know why?

How much are the y,z box dimensions changing during
your simulation, since that is what you are barostatting?
If they are fluctuating improperly that would cause problems.
If you can run w/out a barotstat, and can guess what the
y,z dimension should do as x grows, try running with
deformation on y,z as well, to smoothly adjust it to a final
size you think is close to what you want. If that works,
then the problem must be coming from the barostat
and its behavior. Expecting the y,z barostat to “do the right
thing”, when you are pulling hard in x for a system that
is connected molecules is not always going to work.
It can’t always compensate for large changes to a system.

Steve