Pressure issue in a piston simulation

LAMMPS
,
1

Hello everyone.

I am trying to perform a piston simulation with silica glass (LAMMPS version: 2Aug2023).
The idea is to launch a shock wave inside silica which contains a nanopore in the center and analyze the pressure distribution around the pore and the subsequent pore’s evolution (collapse, shrinkage, elongation, expantion, etc.).

My strategy was this:

  • Amorphization of quartz by melt-quench (NPT)
  • Replication into bulk, 30nm, 30nm, 100nm (NPT)
  • Creation of nanopore by removing spherical number of atoms from center (radius = 10nm) (NPT)
  • Adding vacuum zones before left Z and after right Z boundaries (NVT) and switching from ppp to ppf boundaries
  • Define a piston region of 3 nm thickness at left Z part of material and use “fix move” to linearly move it in positive Z direction

After each of these steps a relaxation was performed as well.

The script I have used for the simulation with piston is the following:

units           metal
dimension       3
boundary        p p f
atom_style      atomic

read_data		d3.lmp

pair_style      vashishta
pair_coeff      * * SiO.1994.vashishta Si O

neighbor        1.0 bin
neigh_modify    every 1 delay 0 check yes

timestep        0.001

thermo          100
thermo_style    custom step temp press pxx pyy pzz pxy pxz pyz pe ke etotal density
thermo_modify	lost	warn

run 0

# ------------------------------------------------------
variable	Nevery		equal		100						
variable	Nrepeat		equal		10						
variable	Nfreq		equal		1000					

variable	deltax		equal		10										
variable	deltay		equal		10
variable	deltaz		equal		10
# ------------------------------------------------------

compute		spa_vir			all		stress/atom		NULL virial
compute		spa_tot			all		stress/atom		NULL

compute		volvoronoi		all		voronoi/atom
compute		KE				all		ke/atom

variable	temp		atom	c_KE/(1.5*8.617333262e-5)

variable	p_vir		atom	-(c_spa_vir[1]+c_spa_vir[2]+c_spa_vir[3])/(3*c_volvoronoi[1]*10000)
variable	p_tot		atom	-(c_spa_tot[1]+c_spa_tot[2]+c_spa_tot[3])/(3*c_volvoronoi[1]*10000)



compute		chunkYZ		all		chunk/atom		bin/3d  x  lower  ${deltax}	y  lower  ${deltay}	z  lower  ${deltaz}		units box

compute		zmatmin		all		reduce min z

run 0

# --------------------------------------------------------------------------

fix		temp_profile	all		ave/chunk	${Nevery}  ${Nrepeat}  ${Nfreq}  chunkYZ  v_temp							file  Temp.out
fix		press_vir		all		ave/chunk	${Nevery}  ${Nrepeat}  ${Nfreq}  chunkYZ  v_p_vir 	c_volvoronoi[1] 		file  Pressure_virial.out
fix		press_tot		all		ave/chunk	${Nevery}  ${Nrepeat}  ${Nfreq}  chunkYZ  v_p_tot 	c_volvoronoi[1] 		file  Pressure_total.out
fix     density_X 		all		ave/chunk	${Nevery}  ${Nrepeat}  ${Nfreq}  chunkYZ  density/number					file  density_X.out
fix     density_m 		all		ave/chunk	${Nevery}  ${Nrepeat}  ${Nfreq}  chunkYZ  density/mass						file  density_m.out
# --------------------------------------------------------------------------

variable        zstart equal c_zmatmin
variable        piston_thick equal 30.0
variable        zend equal v_zstart+v_piston_thick

region          piston block INF INF INF INF ${zstart} ${zend} units box
group           piston region piston
group           mobile subtract all piston

variable        vis_type atom type+2*gmask(piston)

fix             fmob mobile nve
fix             move_piston piston move linear 0.0 0.0 5.0 units box

dump d all custom 1000 random_glass.dump type v_vis_type x y z v_p_vir v_p_tot

run             20000

unfix temp_profile
unfix press_vir
unfix press_tot
unfix density_X
unfix density_m
unfix fmob
unfix move_piston
undump d

Here are the central YZ cut (2 nm thickness) snapshots of the simulation at 0, 2, 4, 6 and 8 ps (piston atoms are just colored differently to detect its movement):

image
image896×3648 2.71 MB
image
image896×3648 2.68 MB
image
image896×3648 2.65 MB
image
image896×3648 2.64 MB
image
image896×3648 2.63 MB

Here are the thermo output from log.lammps at these timesteps:

Step Temp Press Pxx Pyy Pzz Pxy Pxz Pyz PotEng KinEng TotEng Density

0 299.84861 -92.655724 -426.91877 169.46818 -20.516584 -342.48716 70.301908 159.14002 -47832730 247101.44 -47585628 1.6265079

2000 331.5874 5538.2084 2464.0346 3138.1702 11012.42 -305.9637 111.07692 116.43121 -47806557 273256.97 -47533300 1.6265079

4000 366.24102 11070.564 5363.9006 6131.929 21715.862 -279.90045 217.32051 84.681754 -47780596 301814.58 -47478782 1.6265079

6000 407.80975 15270.319 7475.1635 8344.2999 29991.492 -253.7481 251.0739 28.545036 -47760354 336070.85 -47424283 1.6265079

8000 449.65348 18919.753 9245.5801 10055.532 37458.146 -245.73758 308.07806 19.312801 -47740332 370553.74 -47369779 1.6265079

So the average pressure in the system is rising (as I assume due to the compression coming from the piston movement).

However, when I analyze “Pressure_total.out” file, all the chunks have negative ‘meanpress’.
For example, this is a small part from the file at timestep = 4000:

4000 134144 6375417.999999948
  1 5 5 5 0 0 0
  2 5 5 15 0 0 0
  3 5 5 25 0 0 0
  4 5 5 35 0 0 0
  5 5 5 45 0 0 0
  6 5 5 55 0 0 0
  7 5 5 65 0 0 0
  8 5 5 75 0 0 0
  9 5 5 85 0 0 0
  10 5 5 95 0 0 0
  11 5 5 105 0 0 0
  12 5 5 115 16.1 -18.6587 322
  13 5 5 125 64.4 -26.4133 23.9566
  14 5 5 135 60.5 -33.6368 16.1615
  15 5 5 145 65.8 -20.1408 15.7471
  16 5 5 155 66.3 -23.2984 14.3699
  17 5 5 165 69.8 -20.2591 13.7678
  18 5 5 175 74.1 -23.9721 14.3519
  19 5 5 185 67.2 -24.5764 14.9116
  20 5 5 195 63 -26.8098 15.4573
  21 5 5 205 67 -16.5735 15.1864
  22 5 5 215 67.4 -23.0683 14.1267
  23 5 5 225 68.9 -21.6004 13.7351
  24 5 5 235 74.1 -23.9937 14.3833
  25 5 5 245 67.1 -23.8472 14.9961
  26 5 5 255 65.4 -26.6838 15.3292
  27 5 5 265 66.1 -16.0587 15.1527
  28 5 5 275 67.5 -21.9131 14.1692
  29 5 5 285 70.1 -22.4122 13.7528
  30 5 5 295 73.6 -24.5818 14.3574
  31 5 5 305 66.6 -22.668 15.085
  32 5 5 315 67.7 -26.6164 15.2829
  33 5 5 325 63.7 -17.0795 15.1284
  34 5 5 335 67.9 -22.2271 14.1556
  35 5 5 345 70.4 -21.4663 13.843
  36 5 5 355 73.5 -23.9674 14.4152
  37 5 5 365 67.7 -22.2385 15.1415
  38 5 5 375 66.4 -26.0545 15.3424
  39 5 5 385 64.4 -19.3647 15.1383
  40 5 5 395 66.1 -23.8582 14.3493
  41 5 5 405 67.5 -25.1646 14.2503
  42 5 5 415 71.3 -28.1724 14.8306
  43 5 5 425 64.7 -25.1509 15.7212
  44 5 5 435 60.6 -32.4566 15.9678
  45 5 5 445 59.8 -18.4704 16.3612
  46 5 5 455 67.2 -30.6717 14.9568
  47 5 5 465 64 -26.2135 14.5911
  48 5 5 475 68.4 -30.004 15.1633
  49 5 5 485 61.1 -28.4228 15.8285
  50 5 5 495 60.5 -31.7638 16.4574
  51 5 5 505 67.7 -21.5281 16.3187
  52 5 5 515 60.1 -29.3998 15.3211
  53 5 5 525 65.2 -29.3541 14.4579
  54 5 5 535 69.5 -26.7551 15.4428
  55 5 5 545 64.4 -29.1904 15.5219
  56 5 5 555 61.6 -31.6584 16.3276
  57 5 5 565 60.7 -25.2053 16.2075
  58 5 5 575 61.1 -29.433 15.489
  59 5 5 585 64.8 -30.8082 14.7659
  60 5 5 595 64.5 -25.72 15.1077
  61 5 5 605 71.8 -26.3013 15.6224
  62 5 5 615 60.3 -26.3214 16.2304
  63 5 5 625 60.2 -34.3637 16.2053
  64 5 5 635 60.8 -21.4473 16.3988
  65 5 5 645 64.2 -32.2817 14.8753
  66 5 5 655 65.8 -27.6793 14.6071
  67 5 5 665 67 -27.5511 15.3893
  68 5 5 675 59.6 -31.8983 15.4188
  69 5 5 685 63.4 -30.3952 16.9703
  70 5 5 695 64.7 -22.1961 16.012
  71 5 5 705 63.5 -26.7214 15.4914
  72 5 5 715 63.7 -29.0048 14.2471
  73 5 5 725 70 -28.4169 15.3666
  74 5 5 735 65.8 -30.2928 15.6217
  75 5 5 745 59.1 -27.7829 16.5703
  76 5 5 755 62.2 -24.7824 16.2147
  77 5 5 765 58.6 -28.4565 15.5325
  78 5 5 775 66.4 -29.8506 14.8174
  79 5 5 785 65.8 -25.82 15.1544
  80 5 5 795 68.7 -29.6663 15.3219
  81 5 5 805 57.5 -28.6536 16.0155
  82 5 5 815 65.3 -28.2892 16.7128
  83 5 5 825 58.5 -26.1946 15.9455
  84 5 5 835 61.9 -30.2598 15.0312
  85 5 5 845 67 -24.2526 14.8995
  86 5 5 855 68.2 -26.8086 15.4431
  87 5 5 865 60.5 -26.8519 15.8144
  88 5 5 875 58.7 -33.4644 16.6482
  89 5 5 885 66 -28.8987 15.6627
  90 5 5 895 59.7 -24.8619 15.6189
  91 5 5 905 64.8 -29.1338 14.2987
  92 5 5 915 68.2 -25.4512 15.5853
  93 5 5 925 62.8 -31.0277 15.3115
  94 5 5 935 56.9 -27.5638 16.6404
  95 5 5 945 67.3 -27.1446 16.0847
  96 5 5 955 58.9 -26.2478 15.5754
  97 5 5 965 62.6 -27.147 14.8946
  98 5 5 975 65.8 -24.6685 15.2491
  99 5 5 985 64.8 -30.3299 15.2923
  100 5 5 995 61.7 -32.1287 15.8133
  101 5 5 1005 65.6 -26.5368 16.6382
  102 5 5 1015 56.6 -27.3331 15.7665
  103 5 5 1025 62.6 -26.6478 15.1784
  104 5 5 1035 63 -27.7997 14.9498
  105 5 5 1045 66.3 -25.8861 15.5312
  106 5 5 1055 63.9 -30.7878 15.4814
  107 5 5 1065 57.4 -28.9938 16.9569
  108 5 5 1075 60.1 -27.9559 15.8332
  109 5 5 1085 63.7 -27.4459 15.353
  110 5 5 1095 63.6 -30.5535 14.3517
  111 5 5 1105 66.7 -16.956 42.5608
  112 5 5 1115 0 0 0
  113 5 5 1125 0 0 0
  114 5 5 1135 0 0 0
  115 5 5 1145 0 0 0
  116 5 5 1155 0 0 0
  117 5 5 1165 0 0 0
  118 5 5 1175 0 0 0
  119 5 5 1185 0 0 0
  120 5 5 1195 0 0 0
  121 5 5 1205 0 0 0
  122 5 5 1215 0 0 0
  123 5 5 1225 0 0 0
  124 5 5 1235 0 0 0
  125 5 5 1245 0 0 0
  126 5 5 1255 0 0 0
  127 5 5 1265 0 0 0
  128 5 5 1275 0 0 0
  129 5 5 1285 0 0 0
  130 5 5 1295 0 0 0
  131 5 5 1305 0 0 0

So the pressure is negative everywhere, while in a compression simulation like this it should normally be positive, just like in log.lammps.

I also tried to output virial pressure only (“Pressure_virial.out“), but the values are similar to total pressure values:

4000 134144 6375417.999999948
  1 5 5 5 0 0 0
  2 5 5 15 0 0 0
  3 5 5 25 0 0 0
  4 5 5 35 0 0 0
  5 5 5 45 0 0 0
  6 5 5 55 0 0 0
  7 5 5 65 0 0 0
  8 5 5 75 0 0 0
  9 5 5 85 0 0 0
  10 5 5 95 0 0 0
  11 5 5 105 0 0 0
  12 5 5 115 16.1 -18.8151 322
  13 5 5 125 64.4 -26.6516 23.9566
  14 5 5 135 60.5 -33.8927 16.1615
  15 5 5 145 65.8 -20.5277 15.7471
  16 5 5 155 66.3 -23.9435 14.3699
  17 5 5 165 69.8 -20.9463 13.7678
  18 5 5 175 74.1 -24.6593 14.3519
  19 5 5 185 67.2 -25.2087 14.9116
  20 5 5 195 63 -27.4663 15.4573
  21 5 5 205 67 -17.2128 15.1864
  22 5 5 215 67.4 -23.7263 14.1267
  23 5 5 225 68.9 -22.247 13.7351
  24 5 5 235 74.1 -24.6244 14.3833
  25 5 5 245 67.1 -24.4694 14.9961
  26 5 5 255 65.4 -27.3284 15.3292
  27 5 5 265 66.1 -16.6797 15.1527
  28 5 5 275 67.5 -22.5297 14.1692
  29 5 5 285 70.1 -23.0344 13.7528
  30 5 5 295 73.6 -25.1721 14.3574
  31 5 5 305 66.6 -23.2353 15.085
  32 5 5 315 67.7 -27.2089 15.2829
  33 5 5 325 63.7 -17.6702 15.1284
  34 5 5 335 67.9 -22.8275 14.1556
  35 5 5 345 70.4 -22.1259 13.843
  36 5 5 355 73.5 -24.5714 14.4152
  37 5 5 365 67.7 -22.7985 15.1415
  38 5 5 375 66.4 -26.574 15.3424
  39 5 5 385 64.4 -19.8766 15.1383
  40 5 5 395 66.1 -24.3822 14.3493
  41 5 5 405 67.5 -25.6305 14.2503
  42 5 5 415 71.3 -28.5792 14.8306
  43 5 5 425 64.7 -25.5421 15.7212
  44 5 5 435 60.6 -32.8332 15.9678
  45 5 5 445 59.8 -18.8051 16.3612
  46 5 5 455 67.2 -31.0291 14.9568
  47 5 5 465 64 -26.5454 14.5911
  48 5 5 475 68.4 -30.3354 15.1633
  49 5 5 485 61.1 -28.7468 15.8285
  50 5 5 495 60.5 -32.0769 16.4574
  51 5 5 505 67.7 -21.8422 16.3187
  52 5 5 515 60.1 -29.7168 15.3211
  53 5 5 525 65.2 -29.6978 14.4579
  54 5 5 535 69.5 -27.0868 15.4428
  55 5 5 545 64.4 -29.5022 15.5219
  56 5 5 555 61.6 -31.988 16.3276
  57 5 5 565 60.7 -25.5344 16.2075
  58 5 5 575 61.1 -29.7773 15.489
  59 5 5 585 64.8 -31.1621 14.7659
  60 5 5 595 64.5 -26.0302 15.1077
  61 5 5 605 71.8 -26.6365 15.6224
  62 5 5 615 60.3 -26.651 16.2304
  63 5 5 625 60.2 -34.6938 16.2053
  64 5 5 635 60.8 -21.7543 16.3988
  65 5 5 645 64.2 -32.6475 14.8753
  66 5 5 655 65.8 -28.0862 14.6071
  67 5 5 665 67 -27.9065 15.3893
  68 5 5 675 59.6 -32.2555 15.4188
  69 5 5 685 63.4 -30.7086 16.9703
  70 5 5 695 64.7 -22.4987 16.012
  71 5 5 705 63.5 -27.0716 15.4914
  72 5 5 715 63.7 -29.3341 14.2471
  73 5 5 725 70 -28.7386 15.3666
  74 5 5 735 65.8 -30.6259 15.6217
  75 5 5 745 59.1 -28.0838 16.5703
  76 5 5 755 62.2 -25.1137 16.2147
  77 5 5 765 58.6 -28.8029 15.5325
  78 5 5 775 66.4 -30.1658 14.8174
  79 5 5 785 65.8 -26.1662 15.1544
  80 5 5 795 68.7 -29.9912 15.3219
  81 5 5 805 57.5 -28.9792 16.0155
  82 5 5 815 65.3 -28.5747 16.7128
  83 5 5 825 58.5 -26.5126 15.9455
  84 5 5 835 61.9 -30.5906 15.0312
  85 5 5 845 67 -24.5599 14.8995
  86 5 5 855 68.2 -27.1512 15.4431
  87 5 5 865 60.5 -27.174 15.8144
  88 5 5 875 58.7 -33.7997 16.6482
  89 5 5 885 66 -29.2109 15.6627
  90 5 5 895 59.7 -25.1594 15.6189
  91 5 5 905 64.8 -29.4537 14.2987
  92 5 5 915 68.2 -25.7601 15.5853
  93 5 5 925 62.8 -31.3555 15.3115
  94 5 5 935 56.9 -27.8679 16.6404
  95 5 5 945 67.3 -27.472 16.0847
  96 5 5 955 58.9 -26.561 15.5754
  97 5 5 965 62.6 -27.4718 14.8946
  98 5 5 975 65.8 -24.9788 15.2491
  99 5 5 985 64.8 -30.6377 15.2923
  100 5 5 995 61.7 -32.4441 15.8133
  101 5 5 1005 65.6 -26.8644 16.6382
  102 5 5 1015 56.6 -27.6824 15.7665
  103 5 5 1025 62.6 -27.016 15.1784
  104 5 5 1035 63 -28.149 14.9498
  105 5 5 1045 66.3 -26.2152 15.5312
  106 5 5 1055 63.9 -31.1527 15.4814
  107 5 5 1065 57.4 -29.3604 16.9569
  108 5 5 1075 60.1 -28.2977 15.8332
  109 5 5 1085 63.7 -27.7808 15.353
  110 5 5 1095 63.6 -30.9023 14.3517
  111 5 5 1105 66.7 -17.2999 42.5608
  112 5 5 1115 0 0 0
  113 5 5 1125 0 0 0
  114 5 5 1135 0 0 0
  115 5 5 1145 0 0 0
  116 5 5 1155 0 0 0
  117 5 5 1165 0 0 0
  118 5 5 1175 0 0 0
  119 5 5 1185 0 0 0
  120 5 5 1195 0 0 0
  121 5 5 1205 0 0 0
  122 5 5 1215 0 0 0
  123 5 5 1225 0 0 0
  124 5 5 1235 0 0 0
  125 5 5 1245 0 0 0
  126 5 5 1255 0 0 0
  127 5 5 1265 0 0 0
  128 5 5 1275 0 0 0
  129 5 5 1285 0 0 0
  130 5 5 1295 0 0 0
  131 5 5 1305 0 0 0

I also live here the total pressure colormap along Z chunks for central X and Y chunks:

Screenshot_58
Screenshot_581408×671 87.8 KB

Of course, I understand that I have reflection from the rear boundary, but that’s at later times, in the beginning I expect the pressure to be positive like in log.lammps.

So I have three questions:

  1. Why the pressure is negative
  2. What is the correct way to treat the rear surface to model shock wave propagation in bulk silica without reflection? should I dump the velocity of atoms, using fixes like “fix setforce“ or “fix viscous“, or put a second piston at rear surface and fix it’s position or allow it to move dynamically?
  3. What is the overall best way to model shock wave propagation in porous silica? Is piston the best approach, or something like “fix msst“ or maybe other thing I don’t know exists?

Thanks a lot in advance.

2

@davies97 This is a well crafted forum post, but your chances to here anything better than what follows are slim. You need to discuss with an expert in this kind of this kind of simulations and they are rarely the kind that monitor a public forum. We can usually help with general problems or technical details with LAMMPS but questions that are mostly about the science are outside the expertise for almost anybody here.

What follows is a mix of AI-assisted analysis (Claude Opus 4.8) and curation by
a human who does not have specific expertise in shock physics or much experience
in modeling silica, but does have many years of reading and
answering mailing-list and forum posts. Please weigh the suggestions
accordingly, check them against the documentation, and verify with your own
tests rather than taking any of it on faith. Since the thread is still fresh
and hadn’t drawn a detailed answer yet, I figured a careful starting point was
better than none.

Thanks for posting the actual input — that makes this answerable. There are
really two layers here: (A) why your per-chunk pressure comes out negative, which
is a postprocessing issue, and (B) a couple of bigger-picture concerns about
whether the model and cell are set up to answer the question you have in mind.
I’ll take them in that order, then address your rear-boundary / method questions.

A. Why the local pressure is negative

1. The per-chunk estimator is the likely culprit (average-of-ratios)

You build a per-atom pressure first and then let ave/chunk number-average it:

variable p_tot atom -(c_spa_tot[1]+c_spa_tot[2]+c_spa_tot[3])/(3*c_volvoronoi[1]*10000)
fix press_tot all ave/chunk ... chunkYZ v_p_tot c_volvoronoi[1]

That computes the mean of Sᵢ/Vᵢ over the atoms in a bin, which is not the
local pressure. The local virial pressure of a chunk is the ratio of sums, not
the mean of ratios:

P_chunk = -( Σ Sxx_i + Σ Syy_i + Σ Szz_i ) / ( 3 · Σ V_i )

The two agree only if every atom has the same Voronoi volume. They don’t — and
pre-dividing by V_i per atom strongly over-weights atoms with small Voronoi
cells. A handful of those with an opposite-sign virial can dominate the bin
average and flip both its magnitude and its sign. This is the most likely reason
you see a uniformly negative result whose magnitude doesn’t track the bulk
pressure.

The fix needs no new computes — just average the extensive stress components
and the volume separately and combine them afterwards:

fix press_tot all ave/chunk ${Nevery} ${Nrepeat} ${Nfreq} chunkYZ &
    c_spa_tot[1] c_spa_tot[2] c_spa_tot[3] c_volvoronoi[1]

then form -(Sxx+Syy+Szz)/(3*Vsum) in postprocessing. (Because the atom count
cancels in the ratio, using the chunk averages of S and V is equivalent and
avoids the per-atom division entirely.)

2. A decisive cross-check

Before going further, sum the per-atom stress over all atoms and divide by the
total box volume, and compare to thermo Press / Pzz. You already have both
stress/atom NULL (spa_tot, virial+KE) and stress/atom NULL virial
(spa_vir, virial only), so this is essentially free:

  • If the global per-atom sum matches thermo but the chunks don’t → it’s the
    binning/volume estimator in point 1.
  • If even the global doesn’t match thermo → the problem is in the volume
    normalization or the tally, not the binning.
  • Comparing p_tot vs p_vir chunk-by-chunk isolates the kinetic contribution.

3. Remove the streaming velocity from the kinetic term

stress/atom NULL uses lab-frame velocities, so behind the piston the bulk
drift enters the kinetic part of the stress. For a shocked/streaming system the
kinetic stress must be built from the peculiar (thermal) velocities. Pass a
binned, bias-removing temperature compute as the first argument instead of NULL:

compute ctmp all temp/profile 1 0 0 1 z ${nbins}   # or compute temp/chunk ... com yes
compute spa_tot all stress/atom ctmp

This term is positive, so it isn’t what makes your result negative — but as
written the “thermal” part of p_tot is contaminated by the piston motion and
won’t be physically meaningful.

4. Voronoi volumes are ill-defined next to vacuum

With a 10 nm pore and a non-periodic p p f boundary, atoms at the pore wall and
the free z-surfaces have unbounded / clipped Voronoi cells. Those are exactly the
small-or-huge-volume outliers that wreck an average-of-ratios. Restrict the
analysis to a bulk group (exclude a surface/pore skin), or at least flag those
bins.

5. thermo Press is not comparable to a local chunk value

Press in the log is a whole-box average — early in the run that is mostly
unshocked material, plus the box-wide kinetic term and the total volume in the
denominator. A magnitude mismatch versus a local shocked bin is expected; only the
sign flip is a genuine bug.

B. Two concerns about the model and the cell

These are more general and reflect experience rather than anything specific to your
script, but I think they matter more than the pressure bug:

The geometry is a high-porosity periodic lattice, not “bulk + a pore”

A 10 nm-radius pore is 20 nm across in a 30 nm transverse cell, so the wall
(ligament) between a pore and its periodic image is only 30 − 20 = 10 nm, half
the pore diameter, and the areal porosity is about 35%. The shock then couples
to a periodic array of pores with thin walls: release waves off one pore reach
its neighbors across the ligament, and the 10 nm walls deform as ligaments, not as
bulk. If the intent is “an isolated pore in bulk silica,” the transverse cell needs
to be several pore-diameters wide so the ligament is much larger than the pore. As
drawn, the result describes a 35%-porous periodic solid.

Vashishta is least reliable exactly where this problem lives — at the surfaces

Vashishta SiO₂ is a fixed-charge model with an exponentially-screened,
truncated Coulomb term and a three-body bond-bending term, all parameterized for
the bulk environment. A high-porosity cell is mostly surface, and every one of
those assumptions breaks at a pore wall or free surface:

  • No charge equilibration — under-coordinated surface Si/O can’t redistribute
    charge, so the polar surface and any silanol/water passivation are simply absent;
    surface energies and surface stresses are unphysical.
  • The screened Coulomb assumes a surrounding dielectric — at a vacuum interface
    there’s no screening medium, but the bulk screening length is applied anyway.
  • The 3-body term references bulk angles, so under-coordinated wall atoms carry
    spurious bond-bending stress.
  • No way to saturate dangling bonds with H / OH or to add water — the actual
    chemistry of a real silica pore isn’t represented.

Note how this compounds part A: the atoms whose per-atom virial is least
trustworthy (pore rim, free surfaces) are the same ones the average-of-ratios
estimator over-weights — two independent reasons to distrust the surface-region
chunk pressures.

If the science depends on the pore/surface response, a fixed-charge bulk potential
is probably the wrong tool: a reactive / charge-equilibration model (ReaxFF, COMB3)
or a surface-trained machine-learning potential, with hydroxylated/wetted walls,
would be needed — at a substantial cost increase. If you only need bulk shock
response, dropping the pore and using a homogeneous cell keeps Vashishta in its
valid regime (and, conveniently, makes MSST applicable — see below).

C. Rear boundary and overall method (your Q2 / Q3)

  • Avoiding rear-surface reflection (Q2). LAMMPS has no true absorbing boundary.
    Practical options: (a) make the sample long enough that the production window ends
    before the reflected wave returns; (b) add a damped “sponge” layer — fix viscous
    or fix langevin on a rear region — to soak up the wave; (c) avoid the rear
    boundary altogether with MSST (fix msst), which evolves a steady shock in a
    comoving frame with a small cell.

  • Best approach (Q3). For an explicit moving piston, fix wall/piston is the
    purpose-built tool (it provides a proper momentum-mirror reflecting wall), rather
    than fix move on a region group, which is a rigid pusher with no momentum
    mirror. For Hugoniot end-states there’s fix msst and fix nphug. The
    examples/shock directory in the distribution has working setups. Caveat for
    your case:     MSST and nphug assume a homogeneous material and a planar, steady
    shock — with a 10 nm pore the response is inherently non-steady and non-planar,
    so for the porous problem an explicit NEMD piston (with the rear damping-region
    trick) is the appropriate choice and MSST is not. MSST becomes the right tool
    only if you switch to the homogeneous bulk cell mentioned above.

Hope that helps — I’d start with the average-of-ratios fix and the global
cross-check, since those directly settle the sign question, and treat the model /
geometry points as the more important longer-term questions.