[lammps-users] A question about how to compute pressure tensor

Dear all,

  I want to compute pressure tensor. My test system is melt problem in examples folder. (all files is attached) I get the pressure tensor Pzz by two ways. First, from thermo compute pressure. Second, I write my own script to compute pressure tensor, but with same formula as lammps definition (http://lammps.sandia.gov/doc/compute_pressure.html). However, these results are very different. I checked my script many times and could found programming typo. Therefore, I don't know where is wrong? Could any one take a moment to help me figure it out?
  I attached lammps input, and my own python script(usage: ./press.py dump).
  Thanks all in advance

Best wishes,
Yangpeng OU

in.melt (511 Bytes)

press.py (697 Bytes)

Dear all,

  I used SVN version.
  The lammps thermo output is

  Step Press 3[1] 3[3] 3[6]
       0 -3.7033504 -3.7486388 -3.7093134 0.060536765
    1000 5.7734903 5.6886678 5.7903827 0.055700253
    2000 5.9071009 5.865803 5.7800424 0.31273757
    3000 5.8660135 5.9976298 5.5122119 -0.19276738
    4000 5.8100964 5.8223124 5.6771705 0.0938918
    5000 5.7839178 5.5126741 5.8726927 -0.025861736
    6000 5.828503 5.8504812 5.9011465 -0.00047498923
    7000 5.6922059 5.7534908 5.7893515 0.088534491
    8000 5.910558 5.7534054 5.8923415 -0.21761534
    9000 5.7641691 5.5967242 5.9330067 0.12196928
   10000 5.878123 5.9715603 5.6333512 -0.057131362
  3[1] is Pxx
  3[3] is Pzz
  dump.5000 is
      Kinetic part = 1.38922832997
      Potential part = -3.38518718179
      Pzz = -1.99595885182
  dump.6000 is
      Kinetic part = 1.3798937217
      Potential part = 0.288750544165
      Pzz = 1.66864426586
  Attachment is input file and script.

Best wishes,
Yangpeng Ou

in.melt (511 Bytes)

press.py (676 Bytes)

I suggest you look at the details of how LAMMPS computes a pressure - e.g.
pair_lj_cut.cpp and compute_presure.cpp and put in some print
statements if you need to. If you isolate the different terms, you should
be able to find what is the difference compared to your calculation.
E.g. do it for a 2-atom system.

Steve