Skip to content

Fix vibrational/rotational energy sampling cutoff for many DOF (#64) #21

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
stanmoore1 wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Fix vibrational/rotational energy sampling cutoff for many DOF (#64) #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
bc95f93
Fix vibrational/rotational energy sampling cutoff for many DOF
claude Jun 26, 2026
37ca2bd
Apply evib/erot sampling cutoff fix to KOKKOS package
claude Jun 26, 2026
e5071ed
Rebless in.vibrate logs for the energy sampling cutoff fix
claude Jul 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 35 additions & 34 deletions examples/vibrate/log.12Jul24.mpi_1.vibrate → examples/vibrate/log.07Jul26.mpi_1.vibrate
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
SPARTA (7 Mar 2024)
SPARTA (24 Sep 2025)
Running on 1 MPI task(s)
################################################################################
# test of vibrational energy modes
Expand All @@ -22,8 +22,8 @@ create_box 0 1e-5 0 1e-5 0 1e-5
Created orthogonal box = (0 0 0) to (1e-05 1e-05 1e-05)
create_grid 2 2 2
Created 8 child grid cells
CPU time = 0.000811451 secs
create/ghost percent = 98.6133 1.38665
CPU time = 0.000814469 secs
create/ghost percent = 97.2833 2.71674

species co2.species N2 CO2 vibfile co2.species.vib

Expand All @@ -39,7 +39,7 @@ fix 1 vibmode

create_particles mix n 0 twopass
Created 1000 particles
CPU time = 0.00164516 secs
CPU time = 0.00144881 secs

variable collrate equal "ncoll*step*2/np"

Expand All @@ -61,63 +61,64 @@ stats_style step cpu np nattempt ncoll c_St c_Sr c_Srv c_Sv
stats 10

run 100
WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/runner/work/sparta/sparta/src/compute_grid.cpp:203)
WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/user/sparta/src/compute_grid.cpp:203)
Memory usage per proc in Mbytes:
particles (ave,min,max) = 1.9375 1.9375 1.9375
particles (ave,min,max) = 1.8125 1.8125 1.8125
grid (ave,min,max) = 1.51379 1.51379 1.51379
surf (ave,min,max) = 0 0 0
total (ave,min,max) = 3.45288 3.45288 3.45288
total (ave,min,max) = 3.32788 3.32788 3.32788
Step CPU Np Natt Ncoll c_St c_Sr c_Srv c_Sv
0 0 1000 0 0 20527.039 20857.18 0 0
10 0.002129962 1000 39 26 20205.949 20691.401 204.40929 661.99501
20 0.004260896 1000 41 32 19568.703 20782.232 415.73498 964.11516
30 0.006338579 1000 41 20 19231.273 20644.583 622.01036 1289.9613
40 0.008423717 1000 45 30 18809.903 20632.587 774.40678 1458.8993
50 0.010497413 1000 44 28 18453.821 20329.389 988.81379 1699.0558
60 0.012657381 1000 42 30 18346.728 20057.611 1138.0057 1876.238
70 0.014718163 1000 43 25 18087.449 19796.238 1276.5398 2033.2845
80 0.016783203 1000 43 28 18095.118 19426.564 1378.0915 2145.6484
90 0.018804391 1000 43 26 17644.746 19502.423 1565.1087 2347.3082
100 0.020832432 1000 43 31 17636.431 19349.211 1662.1014 2456.254
Loop time of 0.0208741 on 1 procs for 100 steps with 1000 particles
0 0 1000 0 0 20307.245 20991.439 0 0
10 0.001701039 1000 41 31 20136.929 20715.968 153.70002 553.49515
20 0.003346812 1000 42 27 19684.277 20809.904 300.58986 827.67411
30 0.005084473 1000 44 30 19189.196 20798.65 489.49675 1104.0232
40 0.006898945 1000 44 28 18700.169 20948.487 659.87707 1301.6076
50 0.008736291 1000 42 25 18590.564 20466.038 878.81439 1563.6939
60 0.010530764 1000 42 29 17949.78 20613.8 1046.2097 1750.0843
70 0.012353185 1000 43 28 17761.724 20111.423 1200.1754 1922.2281
80 0.014190234 1000 44 26 17863.5 20005.609 1289.8841 2023.4138
90 0.015823542 1000 43 25 17620.526 19846.638 1382.3751 2124.3913
100 0.017426353 1000 45 27 17533.167 19573.719 1453.4421 2205.8191
Loop time of 0.0174433 on 1 procs for 100 steps with 1000 particles
Performance: 5732.876 timesteps/s, 5.733 Mparticle-step/s

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Move | 0.018738 | 0.018738 | 0.018738 | 0.0 | 89.77
Coll | 0.001148 | 0.001148 | 0.001148 | 0.0 | 5.50
Sort | 0.00020238 | 0.00020238 | 0.00020238 | 0.0 | 0.97
Comm | 1.3567e-05 | 1.3567e-05 | 1.3567e-05 | 0.0 | 0.06
Move | 0.015794 | 0.015794 | 0.015794 | 0.0 | 90.54
Coll | 0.00097996 | 0.00097996 | 0.00097996 | 0.0 | 5.62
Sort | 0.00013298 | 0.00013298 | 0.00013298 | 0.0 | 0.76
Comm | 2.2031e-05 | 2.2031e-05 | 2.2031e-05 | 0.0 | 0.13
Modify | 0 | 0 | 0 | 0.0 | 0.00
Output | 0.0007668 | 0.0007668 | 0.0007668 | 0.0 | 3.67
Other | | 4.933e-06 | | | 0.02
Output | 0.00050734 | 0.00050734 | 0.00050734 | 0.0 | 2.91
Other | | 6.989e-06 | | | 0.04

Particle moves = 100000 (0.1M)
Cells touched = 332469 (0.332M)
Cells touched = 328596 (0.329M)
Particle comms = 0 (0K)
Boundary collides = 232448 (0.232M)
Boundary collides = 228741 (0.229M)
Boundary exits = 0 (0K)
SurfColl checks = 0 (0K)
SurfColl occurs = 0 (0K)
Surf reactions = 0 (0K)
Collide attempts = 4239 (4.24K)
Collide occurs = 2938 (2.94K)
Collide attempts = 4218 (4.22K)
Collide occurs = 2876 (2.88K)
Reactions = 0 (0K)
Particles stuck = 0
Axisymm bad moves = 0

Particle-moves/CPUsec/proc: 4.79062e+06
Particle-moves/CPUsec/proc: 5.73288e+06
Particle-moves/step: 1000
Cell-touches/particle/step: 3.32469
Cell-touches/particle/step: 3.28596
Particle comm iterations/step: 1
Particle fraction communicated: 0
Particle fraction colliding with boundary: 2.32448
Particle fraction colliding with boundary: 2.28741
Particle fraction exiting boundary: 0
Surface-checks/particle/step: 0
Surface-collisions/particle/step: 0
Surf-reactions/particle/step: 0
Collision-attempts/particle/step: 0.04239
Collisions/particle/step: 0.02938
Collision-attempts/particle/step: 0.04218
Collisions/particle/step: 0.02876
Reactions/particle/step: 0

Particles: 1000 ave 1000 max 1000 min
Expand Down
77 changes: 39 additions & 38 deletions examples/vibrate/log.12Jul24.mpi_4.vibrate → examples/vibrate/log.07Jul26.mpi_4.vibrate
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
SPARTA (7 Mar 2024)
SPARTA (24 Sep 2025)
Running on 4 MPI task(s)
################################################################################
# test of vibrational energy modes
Expand All @@ -22,8 +22,8 @@ create_box 0 1e-5 0 1e-5 0 1e-5
Created orthogonal box = (0 0 0) to (1e-05 1e-05 1e-05)
create_grid 2 2 2
Created 8 child grid cells
CPU time = 0.000860103 secs
create/ghost percent = 93.2626 6.73745
CPU time = 0.000879556 secs
create/ghost percent = 91.6968 8.30317

species co2.species N2 CO2 vibfile co2.species.vib

Expand All @@ -39,7 +39,7 @@ fix 1 vibmode

create_particles mix n 0 twopass
Created 1000 particles
CPU time = 0.0012523 secs
CPU time = 0.00118294 secs

variable collrate equal "ncoll*step*2/np"

Expand All @@ -61,67 +61,68 @@ stats_style step cpu np nattempt ncoll c_St c_Sr c_Srv c_Sv
stats 10

run 100
WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/runner/work/sparta/sparta/src/compute_grid.cpp:203)
WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/user/sparta/src/compute_grid.cpp:203)
Memory usage per proc in Mbytes:
particles (ave,min,max) = 1.9375 1.9375 1.9375
particles (ave,min,max) = 1.8125 1.8125 1.8125
grid (ave,min,max) = 1.51379 1.51379 1.51379
surf (ave,min,max) = 0 0 0
total (ave,min,max) = 3.45169 3.45169 3.45169
total (ave,min,max) = 3.32669 3.32669 3.32669
Step CPU Np Natt Ncoll c_St c_Sr c_Srv c_Sv
0 0 1000 0 0 19750.722 21250.869 0 0
10 0.001099912 1000 40 29 19584.106 20569.578 179.94805 581.31383
20 0.00204286 1000 45 35 19517.758 20359.032 333.16715 872.99551
30 0.003018749 1000 45 35 19023.886 20314.638 495.38299 1085.8964
40 0.00399013 1000 43 29 18478.529 20230.419 775.61845 1437.0413
50 0.004915165 1000 43 33 18189.411 20164.428 907.17702 1580.6091
60 0.005841441 1000 44 28 17823.085 20197.993 1059.5758 1775.4883
70 0.006737521 1000 44 32 17559.06 19893.231 1235.0641 1970.1942
80 0.007771981 1000 44 29 16907.053 20287.099 1438.5541 2191.7101
90 0.008669063 1000 43 25 16632.298 19978.181 1598.7441 2367.6259
100 0.009600359 1000 45 27 16568.79 19614.025 1741.1069 2529.6526
Loop time of 0.0096493 on 4 procs for 100 steps with 1000 particles
0 0 1000 0 0 20604.382 21301.272 0 0
10 0.001017648 1000 43 32 20499.65 20858.013 150.15022 522.48317
20 0.001881933 1000 41 32 20509.003 20359.192 322.99105 868.56297
30 0.002708559 1000 41 27 20109.663 20119.263 499.14626 1083.5941
40 0.003511703 1000 43 27 19613.153 20167.821 667.20383 1317.1638
50 0.004320494 1000 44 31 19352.856 20248.73 795.75074 1469.0043
60 0.005139096 1000 46 33 18972.622 20395.593 953.88246 1651.0698
70 0.005959587 1000 43 31 18602.262 20130.448 1107.9058 1823.2868
80 0.006810711 1000 43 31 18193.647 20309.316 1280.61 2018.1609
90 0.007636682 1000 45 30 17604.824 20072.144 1467.603 2220.0972
100 0.008497566 1000 43 25 17390.064 19762.296 1727.1455 2495.3107
Loop time of 0.00853214 on 4 procs for 100 steps with 1000 particles
Performance: 11720.384 timesteps/s, 11.720 Mparticle-step/s

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Move | 0.0045529 | 0.0045854 | 0.0045998 | 0.0 | 47.52
Coll | 0.00032341 | 0.00032793 | 0.00033003 | 0.0 | 3.40
Sort | 7.2893e-05 | 7.7937e-05 | 8.4618e-05 | 0.0 | 0.81
Comm | 0.0038536 | 0.0041274 | 0.0042462 | 0.2 | 42.77
Move | 0.0036794 | 0.0038398 | 0.0040027 | 0.2 | 45.00
Coll | 0.00023008 | 0.000236 | 0.00024157 | 0.0 | 2.77
Sort | 4.0392e-05 | 4.2848e-05 | 4.7076e-05 | 0.0 | 0.50
Comm | 0.0039121 | 0.0040888 | 0.0043192 | 0.2 | 47.92
Modify | 0 | 0 | 0 | 0.0 | 0.00
Output | 0.0003919 | 0.00049521 | 0.00078913 | 0.0 | 5.13
Other | | 3.536e-05 | | | 0.37
Output | 0.00022122 | 0.00029706 | 0.00051278 | 0.0 | 3.48
Other | | 2.768e-05 | | | 0.32

Particle moves = 100000 (0.1M)
Cells touched = 326495 (0.326M)
Particle comms = 74165 (74.2K)
Boundary collides = 226537 (0.227M)
Cells touched = 333151 (0.333M)
Particle comms = 74058 (74.1K)
Boundary collides = 233026 (0.233M)
Boundary exits = 0 (0K)
SurfColl checks = 0 (0K)
SurfColl occurs = 0 (0K)
Surf reactions = 0 (0K)
Collide attempts = 4218 (4.22K)
Collide occurs = 2938 (2.94K)
Collide attempts = 4287 (4.29K)
Collide occurs = 2954 (2.95K)
Reactions = 0 (0K)
Particles stuck = 0
Axisymm bad moves = 0

Particle-moves/CPUsec/proc: 2.59086e+06
Particle-moves/CPUsec/proc: 2.9301e+06
Particle-moves/step: 1000
Cell-touches/particle/step: 3.26495
Cell-touches/particle/step: 3.33151
Particle comm iterations/step: 1
Particle fraction communicated: 0.74165
Particle fraction colliding with boundary: 2.26537
Particle fraction communicated: 0.74058
Particle fraction colliding with boundary: 2.33026
Particle fraction exiting boundary: 0
Surface-checks/particle/step: 0
Surface-collisions/particle/step: 0
Surf-reactions/particle/step: 0
Collision-attempts/particle/step: 0.04218
Collisions/particle/step: 0.02938
Collision-attempts/particle/step: 0.04287
Collisions/particle/step: 0.02954
Reactions/particle/step: 0

Particles: 250 ave 267 max 230 min
Histogram: 1 0 1 0 0 0 0 0 0 2
Particles: 250 ave 272 max 206 min
Histogram: 1 0 0 0 0 0 0 1 1 1
Cells: 2 ave 2 max 2 min
Histogram: 4 0 0 0 0 0 0 0 0 0
GhostCell: 6 ave 6 max 6 min
Expand Down
14 changes: 10 additions & 4 deletions src/KOKKOS/particle_kokkos.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -262,9 +262,12 @@ double ParticleKokkos::erot(int isp, double temp_thermal, rand_type &erandom) co
eng = -log(erandom.drand()) * boltz * temp_thermal;
else {
a = 0.5*d_species[isp].rotdof-1.0;
// candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1,
// std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it
// at 10 kT, which is below the mean for large dof
double xmax = a + 1.0 + 9.0*sqrt(a+1.0);
while (1) {
// energy cut-off at 10 kT
erm = 10.0*erandom.drand();
erm = xmax*erandom.drand();
b = pow(erm/a,a) * exp(a-erm);
if (b > erandom.drand()) break;
}
Expand Down Expand Up @@ -297,9 +300,12 @@ double ParticleKokkos::evib(int isp, double temp_thermal, rand_type &erandom) co
eng = -log(erandom.drand()) * boltz * temp_thermal;
else if (d_species[isp].vibdof > 2) {
a = 0.5*d_species[isp].vibdof-1.;
// candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1,
// std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it
// at 10 kT, which is below the mean for large dof
double xmax = a + 1.0 + 9.0*sqrt(a+1.0);
while (1) {
// energy cut-off at 10 kT
erm = 10.0*erandom.drand();
erm = xmax*erandom.drand();
b = pow(erm/a,a) * exp(a-erm);
if (b > erandom.drand()) break;
}
Expand Down
14 changes: 10 additions & 4 deletions src/KOKKOS/surf_collide_diffuse_kokkos.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -330,9 +330,12 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse {
eng = -log(rand_gen.drand()) * boltz * temp_thermal;
} else {
a = 0.5*d_species[isp].rotdof-1.0;
// candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1,
// std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it
// at 10 kT, which is below the mean for large dof
double xmax = a + 1.0 + 9.0*sqrt(a+1.0);
while (1) {
// energy cut-off at 10 kT
erm = 10.0*rand_gen.drand();
erm = xmax*rand_gen.drand();
b = pow(erm/a,a) * exp(a-erm);
if (b > rand_gen.drand()) break;
}
Expand Down Expand Up @@ -370,9 +373,12 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse {
eng = -log(rand_gen.drand()) * boltz * temp_thermal;
else if (d_species[isp].vibdof > 2) {
a = 0.5*d_species[isp].vibdof-1.;
// candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1,
// std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it
// at 10 kT, which is below the mean for large dof
double xmax = a + 1.0 + 9.0*sqrt(a+1.0);
while (1) {
// energy cut-off at 10 kT
erm = 10.0*rand_gen.drand();
erm = xmax*rand_gen.drand();
b = pow(erm/a,a) * exp(a-erm);
if (b > rand_gen.drand()) break;
}
Expand Down
17 changes: 13 additions & 4 deletions src/particle.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1036,9 +1036,13 @@ double Particle::erot(int isp, double temp_thermal, RanKnuth *erandom)
eng = -log(erandom->uniform()) * update->boltz * temp_thermal;
} else {
a = 0.5*particle->species[isp].rotdof-1.0;
// sample E/kT from the equilibrium distribution x^a * exp(-x)
// candidate range must cover the high-energy tail: the mode is at
// x = a, the mean is a+1, the std dev is sqrt(a+1); use mean + ~9 std
// devs so the cut-off scales with rotdof rather than a fixed 10 kT
double xmax = a + 1.0 + 9.0*sqrt(a+1.0);
while (1) {
// energy cut-off at 10 kT
erm = 10.0*erandom->uniform();
erm = xmax*erandom->uniform();
b = pow(erm/a,a) * exp(a-erm);
if (b > erandom->uniform()) break;
}
Expand Down Expand Up @@ -1077,9 +1081,14 @@ double Particle::evib(int isp, double temp_thermal, RanKnuth *erandom)
eng = -log(erandom->uniform()) * update->boltz * temp_thermal;
else if (species[isp].vibdof > 2) {
a = 0.5*particle->species[isp].vibdof-1.;
// sample E/kT from the equilibrium distribution x^a * exp(-x)
// candidate range must cover the high-energy tail: the mode is at
// x = a, the mean is a+1, the std dev is sqrt(a+1); use mean + ~9 std
// devs so the cut-off scales with vibdof (a fixed 10 kT cut-off is
// below the mean for molecules with many vibrational modes)
double xmax = a + 1.0 + 9.0*sqrt(a+1.0);
while (1) {
// energy cut-off at 10 kT
erm = 10.0*erandom->uniform();
erm = xmax*erandom->uniform();
b = pow(erm/a,a) * exp(a-erm);
if (b > erandom->uniform()) break;
}
Expand Down
Loading