I have been trying to preheat a simulation box containing BCC crystal upto 2000K and then perform some studies on the spatial profile of stress. I have been using the following script and the obtained output files are attached. First NPT is ran to bring system to 0 bar and temperature 2000K. Then NVE is ran for some time to hold this temperature. And it accurately does so, as can be seen from log.lammps.
boundary p p p
velocity all create 300 12345 rot yes dist gaussian
fix equilibration all npt temp 300 2100 0.1 iso 0 0 1 drag 1
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol
fix nve all nve
And then to check the spatial profile of stress , box has been binned in Z-direction and centre of mass velocity of each bin has been subtracted to calculate Kinetic energy contribution to stress plus interaction energy contribution as shown, (This has been done since later on I want to check the stress due to impact in z direction). For this following lammps commands has been added to the script.
compute mytemp all temp/partial 1 1 0
compute peratom all stress/atom mytemp ke pair
compute vorvol all voronoi/atom
compute zavg all chunk/atom bin/1d z lower 10 units lattice
variable meanpress atom -(c_peratom+c_peratom+c_peratom)/3/(c_vorvol+1e-99)
variable meanpresszz atom -(c_peratom)/(c_vorvol+1e-99)
variable meanpressxx atom -(c_peratom)/(c_vorvol+1e-99)
fix pres_prof all ave/chunk 100 5 500 zavg v_meanpress file pressure.file
fix presz_prof all ave/chunk 100 5 500 zavg v_meanpresszz file pressurezz.file
fix presx_profile all ave/chunk 100 5 500 zavg v_meanpressxx file pressurexzchunk.file
My question is even if from log.lammps file we can see that the pzz value comes down to around -109 bar which I consider approximately equal to 0 bar. Why is the value of stress so high in z direction when it has been printed binwise. i.e. around 12000 bar. Check pressurezz.file for the chunkwise printing of stress in each bin.
What is that I can be doing wrong or may be understanding wrong! I feel the problem lies in stress/atom calculation. But I can’t figure it out.
Thanks and Regards,
pressurezz.file (4.33 KB)
log.lammps (55.4 KB)
pressurexzchunk.file (4.17 KB)
pressure.file (4.34 KB)