Hello LAMMPS users,
I am simulating tensile deformation of a calcio-silico crystal in LAMMPS, version:14MAY2016. I have equilibrated the system at 1K temperature by NVE+langevin and then by NPT. I am doing the MD run with fix NPT for time integration with fix deform for the deformation of the box. However, there are two issues I ran into:
-
When use fix deform to tension the box, at step 2, the temp rises to 30 K and then go back down.
-
The stress at step 1 is negative. It goes back up to positive and never go down again. Eventually, the material never breaks.
I am using CVFF-CLAYFF potential with morse potential for bond_style. Then I have used fix bond/break with delete_bonds commands but the stress kept on increasing. Note that I have used bond length a bit larger than the equilibrium bond length for fix bond/break. For example, for O-H bond, I used 1 \AA because the equil bond length is 0.96 \AA.
Can someone please tell me what am I doing wrong here? Is there any other way to break the bonds?
Here’s the input file:
Initialization
echo screen
log debug_9tclayffmix.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 13 10 # cutoff1 = 2.5*sigma
read_data data.9tnve
kspace_style pppm 1e-5
delete_atoms overlap 0.4 all all
neighbor 2.0 bin
neigh_modify delay 0 every 1 one 1000 check yes
pair_modify mix arithmetic tail no
#------------ define groups ------------------xxx
group st type 1
group ob type 2
group obts type 3
group oh type 4
group ho type 5
group cao type 6
region leftbdy block INF 2.65 INF INF INF INF units box
group leftbdy region leftbdy
region rightbdy block 17.31 INF INF INF INF INF units box
group rightbdy region rightbdy
group bdy union leftbdy rightbdy
group mobile subtract all bdy
Charge from Cygan et al. (2004)
set type 1 charge 2.1
set type 2 charge -1.05
set type 3 charge -1.1688
set type 4 charge -0.95
set type 5 charge 0.425
set type 6 charge 1.36
compute newtemp all temp
compute poteng all pe
compute bdytemp bdy temp
compute disc all centro/atom fcc
velocity all create 1 1234 temp newtemp dist gaussian mom yes rot yes
neigh_modify exclude type 4 5
run_style verlet
minimize 1.0e-20 1.0e-20 500000 1000000
min_style cg
min_modify dmax 0.1
fix 1 all nve
fix 2 all langevin 1 1 5 12345
thermo 100
thermo_style custom step temp pe etotal press vol c_bdytemp
run 10000
write_data data.9tnve
unfix 1
fix 3 mobile npt temp 1 1 5 iso 0 0 1000 drag 1
#fix 34 all ave/time 20 1 5 c_myRDF[*] file equil.rdf mode vector
dump 3 all xyz 100 t9_equil.xyz
run 100000
write_data data.9tnpt
unfix 2
unfix 3
undump 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]
compute atomid all property/atom id
compute fatom all pair/local dist eng fx fy fz
variable s1 equal c_st[1]
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_sxx1.01e-4/(22.6631.4820.90)
#fix leftfix leftbdy setforce 0 0 0
#fix rightfix rightbdy setforce 0 0 0
#fix leftmove leftbdy move linear -2e-4 0 0 units box
#fix rightmove rightbdy move linear 2e-4 0 0 units box
reset_timestep 0
fix 2 mobile langevin 1 1 5 12345
fix 5 mobile npt temp 1 1 5 y 0 0 1000 z 0 0 1000 drag 1
fix 6 all deform 1 x erate 1e-5 units box remap x
#fix 5 mobile nvt temp 1 1 5
fix 7 all ave/time 20 10 200 v_stress ave window 100 file stress.profile
thermo 100
thermo_style custom step temp v_stress vol press lx ly lz
fix break1 mobile bond/break 10 1 1.8 prob 1 12345
fix break2 mobile bond/break 10 2 1.8 prob 1 12345
fix break3 mobile bond/break 10 3 1.8 prob 1 12345
fix break4 mobile bond/break 10 4 1 prob 1 12345
delete_bonds mobile multi remove
dump 4 all xyz 10 t9_mix.xyz
dump disc all custom 100 disc.csp id type x y z c_disc
run 20000