Hi,I need help please. I have been trying to simulate perfect crystal hexagonal ice using TIP4P under uniaxial loading. I have the unit cell working properly but when I try to replicate it it generates this error of hydrogen missing on proc 0,… I have changed the box but it generated illogical results.
Appreciate you help!
Here is my input file: unit cell atatched below
set the random seed for velocity initialization
variable seed equal 7545
------ System initialization --------
set periodic boundary conditions in 3 directions
boundary p p p
unit in real where energy in kcal/mol and distance in angstroms
units real
include all atom info (charge, atom, type)
atom_style full
bond_style harmonic
angle_style harmonic
adjust atoms to be stored internally as array and sorted
atom_modify map array sort 0 0
read atomic configs
read_data Ih.in
Define bond coefficients for O-H bond
bond_coeff 1 450.0 0.9572 # k = 450 kcal/mol/ ^e , r0 = 0.9572 ^e
Define angle coefficients for H-O-H angle
angle_coeff 1 55.0 104.52 # k = 55.0 kcal/mol/rad
#pair_coeff 1 1 0.1554 3.15365 # Oxygen-Oxygen LJ parameters for TIP4P units #$
Use TIP4P water model
pair_style lj/cut/tip4p/long 2 1 1 1 0.1577 12.0
pair_coeff 2 2 0.1852 3.1589
pair_coeff 1 1 0.0 0.0
pair_coeff 2 1 0 0
#fix shake all shake 0.0001 100 0 b 1 a 1
change_box all triclinic
change_box all x final 0 10 y final 0 20 z final 0 10
replicate 6 6 6
fix shake all shake 0.0001 100 0 b 1 a 1
#delete_atoms overlap 0.5 all all
set type 2 charge -1.1794
set type 1 charge 0.5897
-----Potential and forcefield definitions-----
#pair_modify mix arithmetic
Set up long-range electrostatics
kspace_style pppm/tip4p 1.0e-4 1.0e-4
#delete_atoms overlap 2.8 all all mol yes
reset atom ids and sort them
#reset_atoms id sort yes
----- Variables definition -------
temperature in K
variable T equal 200
pressure in atmospheres
variable P equal 1.0
----- Initialize -----
velocity all create {T} {seed}
set neighbor list parameters
#neighbor 2.0 nsq
neighbor 3.0 bin
neigh_modify delay 0 every 1 check yes
adjust neighbor list updates
#neigh_modify delay 0 every 1 check yes
----- Computes mean square displacement of the center of mass-----
compute msd all msd com yes
----- Log/output -----
thermo 1000
thermo_style custom step atoms etotal ke pe temp press lx ly lz enthalpy density c_msd[4]
thermo_modify norm yes flush yes format float %.6f
----- Output Trajectory -----
dump positions all custom 1000 output_ML-BOP.lammpstrj id type x y z
----- Minimize energy with specified tolerances & iterations -----
minimize 1e-4 1e-6 1000 10000
use conjugate gradient method for minimization
min_style cg
----- Temperature/Pressure Controls -----
fix integrate all npt temp {T} {T} 500 iso {P} {P} 2500
timestep 1 # 5 femtoseconds
run 100000 # no. of steps
unfix integrate # removes the NPT integration fix
undump positions
#fix shake all shake 0.0001 100 0 b 1 a 1
----- Uniaxial loading -----
reset_timestep 0
lz box dimension in z-direction and L0 initial length before stress
variable tmp equal “lz”
variable L0 equal {tmp}
print "Initial Length, L0: {L0}"
variable strain equal “(lz - v_L0)/v_L0”
— captures the strain value
variable p1 equal “v_strain”
— normal stresses
variable p2 equal “-pxx/9869.23266716” # normal stress in x-face
variable p3 equal “-pyy/9869.23266716” # normal stress in y-face
variable p4 equal “-pzz/9869.23266716” # normal stress in z-face
— shear stresses
variable p5 equal “-pxy/9869.23266716”
variable p6 equal “-pxz/9869.23266716”
variable p7 equal “-pyz/9869.23266716”
— captures the step
variable p8 equal "step "
----- Temperature/Pressure Controls for uniaxial loading -----
fix 1 all npt temp {T} {T} 500 x 0 0 2500 y 0 0 2500
----- fix 2a applies uniaxial deformation in z-direction -----
fix 2a all deform 1 z erate 0.00000005 units box remap x
— print to Final.dat
fix def1 all print 200 “{p8} {p1} {p2} {p3} {p4} {p5} {p6} {p7}” file Final.dat screen no
dump 1 all custom 10000 dump.EQ2.*.cfg id type xs ys zs
run 10000000
undump 1
unfix 1