diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index e502480ce5..488de848bf 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -1250,10 +1250,10 @@ tabulate tool .. versionadded:: 22Dec2022 -The ``tabulate`` folder contains Python scripts scripts to generate tabulated -potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module -in the ``tabulate.py`` file. Some example files demonstrating its use are -included. See the README file for more information. +The ``tabulate`` folder contains Python scripts scripts to generate and +visualize tabulated potential files for LAMMPS. The bulk of the code is in the +``tabulate`` module in the ``tabulate.py`` file. Some example files +demonstrating its use are included. See the README file for more information. ---------- @@ -1276,7 +1276,7 @@ Those scripts were written by Steve Plimpton sjplimp at gmail.com valgrind tool ------------- -The ``valgrind`` folder contains additional suppressions fur LAMMPS when +The ``valgrind`` folder contains additional suppressions for LAMMPS when using `valgrind's `_ ` `memcheck tool `_ to search for memory access violation and memory leaks. These suppressions are automatically diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index cfbd01dd99..3eb5b0c61d 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -215,6 +215,9 @@ for an overview of LAMMPS output options. The vector or array will be floating point values that correspond to the specified attribute. +Any settings with the *store/local* option are not saved to a restart +file and must be redefined. + The single() function of this bond style returns 0.0 for the energy of a bonded interaction, since energy is not conserved in these dissipative potentials. It also returns only the normal component of diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 27962c81a1..fb40e5164b 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -215,6 +215,9 @@ for an overview of LAMMPS output options. The vector or array will be floating point values that correspond to the specified attribute. +Any settings with the *store/local* option are not saved to a restart +file and must be redefined. + The potential energy and the single() function of this bond style return :math:`k (r - r_0)^2 / 2` as a proxy of the energy of a bonded interaction, ignoring any volumetric/smoothing factors or dissipative forces. The single() diff --git a/doc/src/compute_rheo_property_atom.rst b/doc/src/compute_rheo_property_atom.rst index 8686a0dec2..fdf1ec8a1f 100644 --- a/doc/src/compute_rheo_property_atom.rst +++ b/doc/src/compute_rheo_property_atom.rst @@ -88,7 +88,7 @@ The *phase* property indicates whether the particle is in a fluid state, a value of 0, or a solid state, a value of 1. The *surface* property indicates the surface designation produced by -the *interface/reconstruct* option of :doc:`fix rheo `. Bulk +the *surface/detection* option of :doc:`fix rheo `. Bulk particles have a value of 0, surface particles have a value of 1, and splash particles have a value of 2. The *surface/r* property is the distance from the surface, up to the kernel cutoff length. Surface particles diff --git a/examples/granular/data.particles b/examples/granular/data.particles deleted file mode 100644 index c9f3bd7a9c..0000000000 --- a/examples/granular/data.particles +++ /dev/null @@ -1,18 +0,0 @@ -Python generated LAMMPS data file - -2 atoms -1 atom types - -0 0.08 xlo xhi -0 0.04 ylo yhi -0 0.08 zlo zhi - -Atoms # sphere - -1 1 0.004 2500 0.04 0.02 0.04 0 0 0 -2 1 0.004 2500 0.04 0.02 0.04416 0 0 0 - -Velocities - -1 0.0 0.0 1 0 0 0 -2 0.0 0.0 -1 0 0 0 diff --git a/examples/granular/in.restitution b/examples/granular/in.restitution index e441ed67a7..e1959c9db3 100644 --- a/examples/granular/in.restitution +++ b/examples/granular/in.restitution @@ -1,24 +1,43 @@ units si atom_style sphere +comm_modify vel yes +boundary p p p -boundary p p f -region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4 +region box block -0.01 0.01 -0.01 0.01 0 0.08 create_box 2 box -read_data data.particles add append -group mb type 1 +create_atoms 1 single 0.0 0.0 0.04 +create_atoms 1 single 0.0 0.0 0.04416 +set group all diameter 0.004 density 2500 -pair_style granular -pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution -# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution -comm_modify vel yes +pair_style granular +pair_coeff * * hertz/material 1e6 0.8 0.4 & + tangential mindlin NULL 0.0 0.5 & + damping coeff_restitution +#pair_coeff * * hooke 1e6 0.5 & +# tangential mindlin 1 1.0 0.0 & +# damping coeff_restitution timestep 1e-9 fix 1 all nve/sphere -compute s all stress/atom NULL pair -#dump 1 all custom 2000000 op.dump id x y z vx vy vz -#dump_modify 1 pad 8 -thermo_style custom step ke -run_style verlet -run 10000000 +group a1 id 1 +group a2 id 2 + +velocity a1 set 0 0 1.0 +velocity a2 set 0 0 -1.0 + +compute z1 a1 reduce sum z +compute z2 a2 reduce sum z +compute v1 a1 reduce sum vz +compute v2 a2 reduce sum vz +compute f1 a1 reduce sum fz +compute f2 a2 reduce sum fz +variable dz equal c_z1 - c_z2 + +# dump 1 all custom 2000000 op.dump id x y z vx vy vz + +thermo 10000 +thermo_style custom step ke v_dz c_v1 c_v2 c_f1 c_f2 + +run 1000000 diff --git a/examples/granular/log.13May23.restitution.g++.1 b/examples/granular/log.13May23.restitution.g++.1 deleted file mode 100644 index e51709d10d..0000000000 --- a/examples/granular/log.13May23.restitution.g++.1 +++ /dev/null @@ -1,80 +0,0 @@ -LAMMPS (17 Apr 2024 - Development - patch_17Apr2024-93-g4e7bddaa0b) - using 1 OpenMP thread(s) per MPI task -units si -atom_style sphere - -boundary p p f -region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4 -create_box 2 box -Created orthogonal box = (0 0 0) to (0.08 0.04 0.08) - 1 by 1 by 1 MPI processor grid - -read_data data.particles add append -Reading data file ... - orthogonal box = (0 0 0) to (0.08 0.04 0.08) - 1 by 1 by 1 MPI processor grid - reading atoms ... - 2 atoms - reading velocities ... - 2 velocities - read_data CPU = 0.002 seconds -group mb type 1 -2 atoms in group mb - -pair_style granular -pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution -# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution -comm_modify vel yes - -timestep 1e-9 -fix 1 all nve/sphere -compute s all stress/atom NULL pair - -#dump 1 all custom 2000000 op.dump id x y z vx vy vz -#dump_modify 1 pad 8 -thermo_style custom step ke -run_style verlet -run 10000000 -Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule -Neighbor list info ... - update: every = 1 steps, delay = 0 steps, check = yes - max neighbors/atom: 2000, page size: 100000 - master list distance cutoff = 0.005 - ghost atom cutoff = 0.005 - binsize = 0.0025, bins = 32 16 32 - 1 neighbor lists, perpetual/occasional/extra = 1 0 0 - (1) pair granular, perpetual - attributes: half, newton on, size, history - pair build: half/size/bin/atomonly/newton - stencil: half/bin/3d - bin: standard -Per MPI rank memory allocation (min/avg/max) = 10.1 | 10.1 | 10.1 Mbytes - Step KinEng - 0 8.3775804e-05 - 10000000 5.3616513e-05 -Loop time of 5.99782 on 1 procs for 10000000 steps with 2 atoms - -77.9% CPU use with 1 MPI tasks x 1 OpenMP threads - -MPI task timing breakdown: -Section | min time | avg time | max time |%varavg| %total ---------------------------------------------------------------- -Pair | 0.60235 | 0.60235 | 0.60235 | 0.0 | 10.04 -Neigh | 0.00021965 | 0.00021965 | 0.00021965 | 0.0 | 0.00 -Comm | 1.7939 | 1.7939 | 1.7939 | 0.0 | 29.91 -Output | 2.5955e-05 | 2.5955e-05 | 2.5955e-05 | 0.0 | 0.00 -Modify | 1.7622 | 1.7622 | 1.7622 | 0.0 | 29.38 -Other | | 1.839 | | | 30.66 - -Nlocal: 2 ave 2 max 2 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Nghost: 0 ave 0 max 0 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Neighs: 0 ave 0 max 0 min -Histogram: 1 0 0 0 0 0 0 0 0 0 - -Total # of neighbors = 0 -Ave neighs/atom = 0 -Neighbor list builds = 14 -Dangerous builds = 0 -Total wall time: 0:00:06 diff --git a/examples/granular/log.4Feb25.restitution.g++.1 b/examples/granular/log.4Feb25.restitution.g++.1 new file mode 100644 index 0000000000..ea2f463b84 --- /dev/null +++ b/examples/granular/log.4Feb25.restitution.g++.1 @@ -0,0 +1,195 @@ +LAMMPS (4 Feb 2025) +units si +atom_style sphere +comm_modify vel yes +boundary p p p + +region box block -0.01 0.01 -0.01 0.01 0 0.08 +create_box 2 box +Created orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08) + 1 by 1 by 1 MPI processor grid + +create_atoms 1 single 0.0 0.0 0.04 +Created 1 atoms + using lattice units in orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08) + create_atoms CPU = 0.000 seconds +create_atoms 1 single 0.0 0.0 0.04416 +Created 1 atoms + using lattice units in orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08) + create_atoms CPU = 0.000 seconds +set group all diameter 0.004 density 2500 +Setting atom values ... + 2 settings made for diameter + 2 settings made for density + +pair_style granular +pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution +#pair_coeff * * hooke 1e6 0.5 # tangential mindlin 1 1.0 0.0 # damping coeff_restitution + +timestep 1e-9 +fix 1 all nve/sphere + +group a1 id 1 +1 atoms in group a1 +group a2 id 2 +1 atoms in group a2 + +velocity a1 set 0 0 1.0 +velocity a2 set 0 0 -1.0 + +compute z1 a1 reduce sum z +compute z2 a2 reduce sum z +compute v1 a1 reduce sum vz +compute v2 a2 reduce sum vz +compute f1 a1 reduce sum fz +compute f2 a2 reduce sum fz +variable dz equal c_z1 - c_z2 + +# dump 1 all custom 2000000 op.dump id x y z vx vy vz + +thermo 10000 +thermo_style custom step ke v_dz c_v1 c_v2 c_f1 c_f2 + +run 1000000 +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 0.005 + ghost atom cutoff = 0.005 + binsize = 0.0025, bins = 8 8 32 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair granular, perpetual + attributes: half, newton on, size, history + pair build: half/size/bin/atomonly/newton + stencil: half/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 9.998 | 9.998 | 9.998 Mbytes + Step KinEng v_dz c_v1 c_v2 c_f1 c_f2 + 0 8.3775804e-05 -0.00416 1 -1 0 0 + 10000 8.3775804e-05 -0.00414 1 -1 0 0 + 20000 8.3775804e-05 -0.00412 1 -1 0 0 + 30000 8.3775804e-05 -0.0041 1 -1 0 0 + 40000 8.3775804e-05 -0.00408 1 -1 0 0 + 50000 8.3775804e-05 -0.00406 1 -1 0 0 + 60000 8.3775804e-05 -0.00404 1 -1 0 0 + 70000 8.3775804e-05 -0.00402 1 -1 0 0 + 80000 8.3775804e-05 -0.004 1 -1 0 0 + 90000 8.3411065e-05 -0.003980019 0.99782075 -0.99782075 -0.023914747 0.023914747 + 100000 8.2852688e-05 -0.0039600945 0.9944753 -0.9944753 -0.032005131 0.032005131 + 110000 8.2139641e-05 -0.0039402463 0.99018672 -0.99018672 -0.039875404 0.039875404 + 120000 8.1272296e-05 -0.0039204934 0.98494496 -0.98494496 -0.048007824 0.048007824 + 130000 8.0246788e-05 -0.0039008551 0.97871113 -0.97871113 -0.056503872 0.056503872 + 140000 7.9058986e-05 -0.0038813518 0.97144075 -0.97144075 -0.065373554 0.065373554 + 150000 7.7705654e-05 -0.0038620047 0.9630903 -0.9630903 -0.074595879 0.074595879 + 160000 7.6184906e-05 -0.0038428357 0.95361959 -0.95361959 -0.084137355 0.084137355 + 170000 7.4496418e-05 -0.0038238676 0.94299284 -0.94299284 -0.093958893 0.093958893 + 180000 7.2641536e-05 -0.0038051239 0.93117907 -0.93117907 -0.10401872 0.10401872 + 190000 7.0623328e-05 -0.0037866285 0.91815243 -0.91815243 -0.11427372 0.11427372 + 200000 6.8446602e-05 -0.003768406 0.90389221 -0.90389221 -0.12468011 0.12468011 + 210000 6.6117901e-05 -0.0037504812 0.88838298 -0.88838298 -0.13519381 0.13519381 + 220000 6.3645478e-05 -0.0037328791 0.87161455 -0.87161455 -0.14577066 0.14577066 + 230000 6.1039243e-05 -0.003715625 0.85358204 -0.85358204 -0.1563666 0.1563666 + 240000 5.8310702e-05 -0.0036987442 0.83428576 -0.83428576 -0.16693776 0.16693776 + 250000 5.5472871e-05 -0.003682262 0.8137313 -0.8137313 -0.17744066 0.17744066 + 260000 5.2540172e-05 -0.0036662033 0.79192936 -0.79192936 -0.18783225 0.18783225 + 270000 4.9528314e-05 -0.003650593 0.76889577 -0.76889577 -0.19807012 0.19807012 + 280000 4.6454158e-05 -0.0036354555 0.74465137 -0.74465137 -0.20811258 0.20811258 + 290000 4.3335566e-05 -0.0036208149 0.71922195 -0.71922195 -0.21791884 0.21791884 + 300000 4.0191232e-05 -0.0036066944 0.69263807 -0.69263807 -0.22744912 0.22744912 + 310000 3.7040511e-05 -0.0035931168 0.66493499 -0.66493499 -0.23666485 0.23666485 + 320000 3.390323e-05 -0.0035801042 0.63615249 -0.63615249 -0.24552878 0.24552878 + 330000 3.0799488e-05 -0.0035676776 0.60633473 -0.60633473 -0.25400513 0.25400513 + 340000 2.7749462e-05 -0.0035558573 0.57553002 -0.57553002 -0.26205979 0.26205979 + 350000 2.4773197e-05 -0.0035446626 0.54379063 -0.54379063 -0.2696604 0.2696604 + 360000 2.1890403e-05 -0.0035341116 0.51117262 -0.51117262 -0.27677654 0.27677654 + 370000 1.9120254e-05 -0.0035242212 0.47773551 -0.47773551 -0.28337983 0.28337983 + 380000 1.6481181e-05 -0.0035150072 0.44354212 -0.44354212 -0.28944409 0.28944409 + 390000 1.3990689e-05 -0.0035064841 0.40865822 -0.40865822 -0.29494544 0.29494544 + 400000 1.1665166e-05 -0.003498665 0.37315233 -0.37315233 -0.2998624 0.2998624 + 410000 9.5197195e-06 -0.0034915617 0.33709536 -0.33709536 -0.304176 0.304176 + 420000 7.5680136e-06 -0.0034851844 0.30056032 -0.30056032 -0.30786987 0.30786987 + 430000 5.8221324e-06 -0.003479542 0.26362205 -0.26362205 -0.31093026 0.31093026 + 440000 4.2924559e-06 -0.0034746417 0.22635684 -0.22635684 -0.31334618 0.31334618 + 450000 2.9875585e-06 -0.0034704894 0.18884214 -0.18884214 -0.31510936 0.31510936 + 460000 1.9141264e-06 -0.0034670891 0.15115621 -0.15115621 -0.31621432 0.31621432 + 470000 1.0768988e-06 -0.0034644437 0.11337783 -0.11337783 -0.31665837 0.31665837 + 480000 4.786302e-07 -0.0034625541 0.075585893 -0.075585893 -0.31644161 0.31644161 + 490000 1.2007709e-07 -0.0034614198 0.037859142 -0.037859142 -0.31556689 0.31556689 + 500000 6.3727744e-12 -0.0034610388 0.0002758068 -0.0002758068 -0.31403982 0.31403982 + 510000 1.1522726e-07 -0.0034614073 -0.03708671 0.03708671 -0.31186864 0.31186864 + 520000 4.6064472e-07 -0.0034625203 -0.074152149 0.074152149 -0.30906426 0.30906426 + 530000 1.029334e-06 -0.0034643709 -0.1108457 0.1108457 -0.30564009 0.30564009 + 540000 1.812635e-06 -0.0034669512 -0.14709431 0.14709431 -0.30161199 0.30161199 + 550000 2.8002645e-06 -0.0034702513 -0.18282695 0.18282695 -0.29699817 0.29699817 + 560000 3.9804448e-06 -0.0034742603 -0.21797491 0.21797491 -0.29181905 0.29181905 + 570000 5.3400475e-06 -0.0034789659 -0.25247202 0.25247202 -0.28609717 0.28609717 + 580000 6.8647484e-06 -0.0034843545 -0.28625495 0.28625495 -0.27985701 0.27985701 + 590000 8.5391931e-06 -0.003490411 -0.31926339 0.31926339 -0.2731249 0.2731249 + 600000 1.0347171e-05 -0.0034971194 -0.35144026 0.35144026 -0.2659288 0.2659288 + 610000 1.2271792e-05 -0.0035044627 -0.38273192 0.38273192 -0.25829822 0.25829822 + 620000 1.429567e-05 -0.0035124225 -0.41308835 0.41308835 -0.25026402 0.25026402 + 630000 1.6401103e-05 -0.0035209797 -0.44246327 0.44246327 -0.24185823 0.24185823 + 640000 1.8570255e-05 -0.0035301142 -0.47081428 0.47081428 -0.23311393 0.23311393 + 650000 2.078533e-05 -0.0035398052 -0.498103 0.498103 -0.22406507 0.22406507 + 660000 2.3028745e-05 -0.003550031 -0.52429514 0.52429514 -0.21474628 0.21474628 + 670000 2.5283288e-05 -0.0035607695 -0.54936056 0.54936056 -0.20519274 0.20519274 + 680000 2.7532278e-05 -0.0035719977 -0.57327337 0.57327337 -0.19544002 0.19544002 + 690000 2.9759697e-05 -0.0035836926 -0.59601193 0.59601193 -0.18552392 0.18552392 + 700000 3.1950329e-05 -0.0035958303 -0.61755887 0.61755887 -0.17548034 0.17548034 + 710000 3.408987e-05 -0.0036083869 -0.63790112 0.63790112 -0.16534511 0.16534511 + 720000 3.6165032e-05 -0.0036213382 -0.65702988 0.65702988 -0.15515392 0.15515392 + 730000 3.8163631e-05 -0.00363466 -0.67494058 0.67494058 -0.14494218 0.14494218 + 740000 4.0074659e-05 -0.0036483277 -0.69163285 0.69163285 -0.13474489 0.13474489 + 750000 4.1888343e-05 -0.0036623172 -0.7071105 0.7071105 -0.1245966 0.1245966 + 760000 4.3596185e-05 -0.0036766041 -0.7213814 0.7213814 -0.11453133 0.11453133 + 770000 4.5190996e-05 -0.0036911645 -0.73445747 0.73445747 -0.1045825 0.1045825 + 780000 4.6666906e-05 -0.0037059745 -0.74635458 0.74635458 -0.094782939 0.094782939 + 790000 4.8019368e-05 -0.0037210109 -0.75709246 0.75709246 -0.085164857 0.085164857 + 800000 4.9245153e-05 -0.0037362507 -0.76669467 0.76669467 -0.075759891 0.075759891 + 810000 5.0342325e-05 -0.0037516713 -0.77518852 0.77518852 -0.066599178 0.066599178 + 820000 5.1310215e-05 -0.003767251 -0.78260499 0.78260499 -0.057713474 0.057713474 + 830000 5.2149388e-05 -0.0037829686 -0.78897875 0.78897875 -0.049133335 0.049133335 + 840000 5.2861601e-05 -0.0037988035 -0.79434809 0.79434809 -0.040889384 0.040889384 + 850000 5.3449761e-05 -0.0038147361 -0.79875498 0.79875498 -0.03301269 0.03301269 + 860000 5.3917883e-05 -0.0038307476 -0.80224517 0.80224517 -0.025535302 0.025535302 + 870000 5.4271056e-05 -0.0038468201 -0.80486832 0.80486832 -0.018491033 0.018491033 + 880000 5.4515415e-05 -0.0038629369 -0.80667827 0.80667827 -0.011916591 0.011916591 + 890000 5.4658137e-05 -0.0038790822 -0.80773352 0.80773352 -0.00585333 0.00585333 + 900000 5.4707463e-05 -0.0038952416 -0.80809791 0.80809791 -0.00035001143 0.00035001143 + 910000 5.4672786e-05 -0.003911402 -0.80784176 0.80784176 0.0045325009 -0.0045325009 + 920000 5.4564826e-05 -0.0039275517 -0.80704376 0.80704376 0.0087126344 -0.0087126344 + 930000 5.4396015e-05 -0.0039436807 -0.80579439 0.80579439 0.012070422 -0.012070422 + 940000 5.4181322e-05 -0.0039597811 -0.80420264 0.80420264 0.014404618 -0.014404618 + 950000 5.3940325e-05 -0.0039758475 -0.80241211 0.80241211 0.015295834 -0.015295834 + 960000 5.3704817e-05 -0.0039918778 -0.8006585 0.8006585 0.013305187 -0.013305187 + 970000 5.3616513e-05 -0.0040078808 -0.79999999 0.79999999 0 0 + 980000 5.3616513e-05 -0.0040238808 -0.79999999 0.79999999 0 0 + 990000 5.3616513e-05 -0.0040398808 -0.79999999 0.79999999 0 0 + 1000000 5.3616513e-05 -0.0040558808 -0.79999999 0.79999999 0 0 +Loop time of 0.572025 on 1 procs for 1000000 steps with 2 atoms + +99.3% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0.1248 | 0.1248 | 0.1248 | 0.0 | 21.82 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0.21166 | 0.21166 | 0.21166 | 0.0 | 37.00 +Output | 0.00057419 | 0.00057419 | 0.00057419 | 0.0 | 0.10 +Modify | 0.12695 | 0.12695 | 0.12695 | 0.0 | 22.19 +Other | | 0.108 | | | 18.89 + +Nlocal: 2 ave 2 max 2 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 1 ave 1 max 1 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 1 +Ave neighs/atom = 0.5 +Neighbor list builds = 0 +Dangerous builds = 0 +Total wall time: 0:00:01 diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index db6b96cf05..284626464e 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -397,6 +397,7 @@ double BondBPM::equilibrium_distance(int /*i*/) void BondBPM::write_restart(FILE *fp) { fwrite(&overlay_flag, sizeof(int), 1, fp); + fwrite(&break_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -405,8 +406,12 @@ void BondBPM::write_restart(FILE *fp) void BondBPM::read_restart(FILE *fp) { - if (comm->me == 0) utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error); + if (comm->me == 0) { + utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &break_flag, sizeof(int), 1, fp, nullptr, error); + } MPI_Bcast(&overlay_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&break_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- */ diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index 7cc04955c6..b1e6163d97 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -28,6 +28,7 @@ #include "memory.h" #include "modify.h" #include "neighbor.h" +#include "update.h" #include #include @@ -483,6 +484,7 @@ void BondBPMRotational::compute(int eflag, int vflag) int newton_bond = force->newton_bond; double **bondstore = fix_bond_history->bondstore; + const bool allow_breaks = (update->setupflag == 0) && break_flag; for (n = 0; n < nbondlist; n++) { @@ -527,7 +529,7 @@ void BondBPMRotational::compute(int eflag, int vflag) breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2, torque1on2, torque2on1); - if ((breaking >= 1.0) && break_flag) { + if ((breaking >= 1.0) && allow_breaks) { bondlist[n][2] = 0; process_broken(i1, i2); continue; @@ -763,6 +765,7 @@ void BondBPMRotational::read_restart(FILE *fp) void BondBPMRotational::write_restart_settings(FILE *fp) { fwrite(&smooth_flag, sizeof(int), 1, fp); + fwrite(&normalize_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -771,8 +774,12 @@ void BondBPMRotational::write_restart_settings(FILE *fp) void BondBPMRotational::read_restart_settings(FILE *fp) { - if (comm->me == 0) utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error); + if (comm->me == 0) { + utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &normalize_flag, sizeof(int), 1, fp, nullptr, error); + } MPI_Bcast(&smooth_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&normalize_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- */ diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index 873b5010b4..a5c31ea1b2 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/src/BPM/bond_bpm_spring.cpp @@ -26,6 +26,7 @@ #include "memory.h" #include "modify.h" #include "neighbor.h" +#include "update.h" #include #include @@ -218,6 +219,7 @@ void BondBPMSpring::compute(int eflag, int vflag) double invdim = 1.0 / dim; double **bondstore = fix_bond_history->bondstore; + const bool allow_breaks = (update->setupflag == 0) && break_flag; for (n = 0; n < nbondlist; n++) { @@ -249,7 +251,7 @@ void BondBPMSpring::compute(int eflag, int vflag) r = sqrt(rsq); e = (r - r0) / r0; - if ((fabs(e) > ecrit[type]) && break_flag) { + if ((fabs(e) > ecrit[type]) && allow_breaks) { bondlist[n][2] = 0; process_broken(i1, i2); diff --git a/src/BPM/bond_bpm_spring_plastic.cpp b/src/BPM/bond_bpm_spring_plastic.cpp index 7059449f3b..21ef0f12fb 100644 --- a/src/BPM/bond_bpm_spring_plastic.cpp +++ b/src/BPM/bond_bpm_spring_plastic.cpp @@ -26,6 +26,7 @@ #include "force.h" #include "memory.h" #include "neighbor.h" +#include "update.h" #include #include @@ -185,6 +186,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag) int newton_bond = force->newton_bond; double **bondstore = fix_bond_history->bondstore; + const bool allow_breaks = (update->setupflag == 0) && break_flag; for (n = 0; n < nbondlist; n++) { @@ -217,7 +219,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag) r = sqrt(rsq); e = (r - r0) / r0; - if ((fabs(e) > ecrit[type]) && break_flag) { + if ((fabs(e) > ecrit[type]) && allow_breaks) { bondlist[n][2] = 0; process_broken(i1, i2); continue; diff --git a/src/BPM/fix_update_special_bonds.cpp b/src/BPM/fix_update_special_bonds.cpp index 3c37769fee..b05f5196cc 100644 --- a/src/BPM/fix_update_special_bonds.cpp +++ b/src/BPM/fix_update_special_bonds.cpp @@ -26,6 +26,7 @@ #include "neigh_list.h" #include "neighbor.h" +#include #include using namespace LAMMPS_NS; @@ -267,18 +268,40 @@ void FixUpdateSpecialBonds::post_run() void FixUpdateSpecialBonds::add_broken_bond(int i, int j) { - auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]); + tagint *tag = atom->tag; + int mintag = MIN(tag[i], tag[j]); + int maxtag = MAX(tag[i], tag[j]); + auto tag_pair = std::make_pair(mintag, maxtag); new_broken_pairs.push_back(tag_pair); - broken_pairs.push_back(tag_pair); + + // cancel out if created->destroyed before nlist rebuild + // however, still cannot break + create in the same timestep + auto const &it = created_pairs.find(tag_pair); + if (it != created_pairs.end()) { + created_pairs.erase(it); + } else { + broken_pairs.insert(tag_pair); + } } /* ---------------------------------------------------------------------- */ void FixUpdateSpecialBonds::add_created_bond(int i, int j) { - auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]); + tagint *tag = atom->tag; + int mintag = MIN(tag[i], tag[j]); + int maxtag = MAX(tag[i], tag[j]); + auto tag_pair = std::make_pair(mintag, maxtag); new_created_pairs.push_back(tag_pair); - created_pairs.push_back(tag_pair); + + // cancel out if destroyed->created before nlist rebuild + // however, still cannot break + create in the same timestep + auto const &it = broken_pairs.find(tag_pair); + if (it != broken_pairs.end()) { + broken_pairs.erase(it); + } else { + created_pairs.insert(tag_pair); + } } /* ---------------------------------------------------------------------- diff --git a/src/BPM/fix_update_special_bonds.h b/src/BPM/fix_update_special_bonds.h index b244c21a15..12e775176e 100644 --- a/src/BPM/fix_update_special_bonds.h +++ b/src/BPM/fix_update_special_bonds.h @@ -22,6 +22,7 @@ FixStyle(UPDATE_SPECIAL_BONDS,FixUpdateSpecialBonds); #include "fix.h" +#include #include namespace LAMMPS_NS { @@ -39,15 +40,15 @@ class FixUpdateSpecialBonds : public Fix { void write_restart(FILE *) override; protected: - // Create two arrays to store bonds broken this timestep (new) - // and since the last neighbor list build + // Create array to store bonds broken this timestep (new) + // and a set for those broken since the last neighbor list build std::vector> new_broken_pairs; - std::vector> broken_pairs; + std::set> broken_pairs; - // Create two arrays to store newly created this timestep (new) - // and since the last neighbor list build + // Create arrays to store newly created this timestep (new) + // and a set for those created since the last neighbor list build std::vector> new_created_pairs; - std::vector> created_pairs; + std::set> created_pairs; }; } // namespace LAMMPS_NS diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index f7b74dbf70..d5b2dab1b0 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -44,7 +44,7 @@ using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : - BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr), + BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), id_fix(nullptr), compute_surface(nullptr) { partial_flag = 1; @@ -71,7 +71,6 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix)); index_nb = atom->find_custom("shell_nbond", tmp1, tmp2); } - nbond = atom->ivector[index_nb]; //Store non-persistent per atom quantities, intermediate @@ -181,6 +180,7 @@ void BondRHEOShell::compute(int eflag, int vflag) double **v = atom->v; double **f = atom->f; tagint *tag = atom->tag; + int *nbond = atom->ivector[index_nb]; int *status = atom->rheo_status; int **bondlist = neighbor->bondlist; int nbondlist = neighbor->nbondlist; diff --git a/src/RHEO/bond_rheo_shell.h b/src/RHEO/bond_rheo_shell.h index 828f693ea3..6640cab2ad 100644 --- a/src/RHEO/bond_rheo_shell.h +++ b/src/RHEO/bond_rheo_shell.h @@ -45,7 +45,7 @@ class BondRHEOShell : public BondBPM { double *k, *ecrit, *gamma; double tform, rmax, rsurf; - int *dbond, *nbond; + int *dbond; int index_nb, nmax_store; char *id_fix; diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 9c92e21166..0cd6e46417 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -39,12 +39,10 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace MathExtra; -static constexpr double EPSILON = 1e-1; - /* ---------------------------------------------------------------------- */ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr), rho0(nullptr), + Compute(lmp, narg, arg), chi(nullptr), fix_rheo(nullptr), rho0(nullptr), norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), compute_kernel(nullptr), fix_pressure(nullptr) @@ -64,13 +62,12 @@ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) : // between timesteps (fix property atom will handle callbacks) int tmp1, tmp2; - int index = atom->find_custom("fp_store", tmp1, tmp2); - if (index == -1) { + index_fp_store = atom->find_custom("fp_store", tmp1, tmp2); + if (index_fp_store == -1) { id_fix_pa = utils::strdup(id + std::string("_fix_property_atom")); modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa)); - index = atom->find_custom("fp_store", tmp1, tmp2); + index_fp_store = atom->find_custom("fp_store", tmp1, tmp2); } - fp_store = atom->darray[index]; } /* ---------------------------------------------------------------------- */ @@ -114,7 +111,7 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr) void ComputeRHEOInterface::compute_peratom() { - int a, i, j, ii, jj, jnum, itype, fluidi, fluidj, status_match; + int a, i, j, ii, jj, jnum, fluidi, fluidj, status_match; double xtmp, ytmp, ztmp, rsq, w, dot, dx[3]; int inum, *ilist, *jlist, *numneigh, **firstneigh; @@ -126,6 +123,7 @@ void ComputeRHEOInterface::compute_peratom() int newton = force->newton; int *status = atom->rheo_status; double *rho = atom->rho; + double **fp_store = atom->darray[index_fp_store]; inum = list->inum; ilist = list->ilist; @@ -151,7 +149,6 @@ void ComputeRHEOInterface::compute_peratom() xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; - itype = type[i]; fluidi = !(status[i] & PHASECHECK); jlist = firstneigh[i]; jnum = numneigh[i]; @@ -212,9 +209,10 @@ void ComputeRHEOInterface::compute_peratom() if (status[i] & PHASECHECK) { if (normwf[i] != 0.0) { // Stores rho for solid particles 1+Pw in Adami Adams 2012 - rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], i)); + // cap out at a tenth of equilibrium + rho[i] = MAX(0.1 * rho0[type[i]], fix_pressure->calc_rho(rho[i] / normwf[i], i)); } else { - rho[i] = rho0[itype]; + rho[i] = rho0[type[i]]; } } } @@ -230,6 +228,7 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int / { int m = 0; double *rho = atom->rho; + double **fp_store = atom->darray[index_fp_store]; for (int i = 0; i < n; i++) { int j = list[i]; @@ -250,6 +249,8 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int / void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf) { double *rho = atom->rho; + double **fp_store = atom->darray[index_fp_store]; + int m = 0; int last = first + n; for (int i = first; i < last; i++) { @@ -335,6 +336,7 @@ void ComputeRHEOInterface::store_forces() double *mass = atom->mass; double *rmass = atom->rmass; double **f = atom->f; + double **fp_store = atom->darray[index_fp_store]; // When this is called, fp_store stores the pressure force // After this method, fp_store instead stores non-pressure forces diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h index 2f88eca50f..c890f70595 100644 --- a/src/RHEO/compute_rheo_interface.h +++ b/src/RHEO/compute_rheo_interface.h @@ -40,7 +40,8 @@ class ComputeRHEOInterface : public Compute { double correct_rho(int); void store_forces(); - double *chi, **fp_store; + double *chi; + int index_fp_store; class FixRHEO *fix_rheo; private: diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 3a6629bfef..61aad6a9b8 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -220,6 +220,8 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji); if (kernel_style == QUINTIC) return calc_dw_quintic(delx, dely, delz, r, dWij, dWji); + if (kernel_style == RK0) + return calc_dw_quintic(delx, dely, delz, r, dWij, dWji); double wp; int corrections_i = check_corrections(i); diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 6d122a2ec1..a6226dbb8d 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -502,7 +502,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n) void ComputeRHEOPropertyAtom::pack_nbond_shell(int n) { - int *nbond = fix_oxidation->nbond; + int *nbond = atom->ivector[fix_oxidation->index_nb]; int *mask = atom->mask; int nlocal = atom->nlocal; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 56738c7a0c..c2ff2d8b09 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -39,6 +39,7 @@ using namespace RHEO_NS; using namespace FixConst; static const char cite_rheo[] = + "RHEO package: doi:10.1063/5.0228823\n\n" "@article{Palermo2024,\n" " journal = {Physics of Fluids},\n" " title = {Reproducing hydrodynamics and elastic objects: A hybrid mesh-free model framework for dynamic multi-phase flows},\n" @@ -117,17 +118,17 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : shift_flag = 1; memory->create(shift_type, n + 1, "rheo:shift_type"); for (i = 1; i <= n; i++) shift_type[i] = 1; - while (iarg < narg) { // optional sub-arguments - if (strcmp(arg[iarg], "scale/cross/type") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error); + while (iarg + 1 < narg) { // optional sub-arguments + if (strcmp(arg[iarg + 1], "scale/cross/type") == 0) { + if (iarg + 4 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error); shift_cross_type_flag = 1; - shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); - shift_cmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - shift_wmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + shift_scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + shift_wmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg], "exclude/type") == 0) { - if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); - utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error); + } else if (strcmp(arg[iarg + 1], "exclude/type") == 0) { + if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); + utils::bounds(FLERR, arg[iarg + 2], 1, n, nlo, nhi, error); for (i = nlo; i <= nhi; i++) shift_type[i] = 0; iarg += 1; } else { diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp index ab21307323..5c76094435 100644 --- a/src/RHEO/fix_rheo_oxidation.cpp +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -38,6 +38,7 @@ using namespace FixConst; enum { NONE, CONSTANT }; static const char cite_rheo_oxide[] = + "RHEO oxidation: doi:10.1016/j.apm.2024.02.027\n\n" "@article{ApplMathModel.130.310,\n" " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n" " journal = {Applied Mathematical Modelling},\n" @@ -109,7 +110,6 @@ void FixRHEOOxidation::init() int tmp1, tmp2; index_nb = atom->find_custom("shell_nbond", tmp1, tmp2); if (index_nb == -1) error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation"); - nbond = atom->ivector[index_nb]; // need a half neighbor list auto req = neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/RHEO/fix_rheo_oxidation.h b/src/RHEO/fix_rheo_oxidation.h index 6dddea867a..f845e1c064 100644 --- a/src/RHEO/fix_rheo_oxidation.h +++ b/src/RHEO/fix_rheo_oxidation.h @@ -39,11 +39,11 @@ class FixRHEOOxidation : public Fix { void post_force(int) override; int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; - int *nbond; double rsurf, cut; + int index_nb; private: - int btype, index_nb; + int btype; double cutsq; class NeighList *list; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index f8449ead70..e9842e19e2 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -221,11 +221,6 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : FixRHEOThermal::~FixRHEOThermal() { - // Remove custom property if it exists - int tmp1, tmp2, index; - index = atom->find_custom("rheo_conductivity", tmp1, tmp2); - if (index != -1) atom->remove_custom(index, 1, 0); - memory->destroy(cv_style); memory->destroy(Tc_style); memory->destroy(kappa_style); diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 0fdd711ba5..db83b01539 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -115,7 +115,7 @@ void PairRHEO::compute(int eflag, int vflag) double **fp_store, *chi; if (compute_interface) { - fp_store = compute_interface->fp_store; + fp_store = atom->darray[compute_interface->index_fp_store]; chi = compute_interface->chi; for (i = 0; i < atom->nmax; i++) { @@ -539,7 +539,7 @@ double PairRHEO::init_one(int i, int j) int PairRHEO::pack_reverse_comm(int n, int first, double *buf) { - double **fp_store = compute_interface->fp_store; + double **fp_store = atom->darray[compute_interface->index_fp_store]; int m = 0; int last = first + n; for (int i = first; i < last; i++) { @@ -555,7 +555,7 @@ int PairRHEO::pack_reverse_comm(int n, int first, double *buf) void PairRHEO::unpack_reverse_comm(int n, int *list, double *buf) { - double **fp_store = compute_interface->fp_store; + double **fp_store = atom->darray[compute_interface->index_fp_store]; int m = 0; for (int i = 0; i < n; i++) { int j = list[i]; diff --git a/tools/tabulate/README.md b/tools/tabulate/README.md index d3b4df0738..995bb8abdd 100644 --- a/tools/tabulate/README.md +++ b/tools/tabulate/README.md @@ -27,6 +27,7 @@ Please see the individual tabulation scripts in this folder for examples: | wall_harmonic_tabulate.py | creates a table for fix wall/table with a simple repulsive harmonic potential | | wall_multi_tabulate.py | creates a table for fix wall/table with multiple tables | | pair_bi_tabulate.py | creates a table from radial distribution file using Boltzmann Inversion | +| plot_forces.py | plots and extracts tabulated forces from table files | Common command line flags: diff --git a/tools/tabulate/plot_forces.py b/tools/tabulate/plot_forces.py new file mode 100755 index 0000000000..ee55a93be6 --- /dev/null +++ b/tools/tabulate/plot_forces.py @@ -0,0 +1,278 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# Author: Germain Clavier (Unicaen), germain.clavier at unicaen.fr + +""" +Plot LAMMPS tabulated forces. +""" + +import argparse +import numpy as np +import os +import logging +import sys +from matplotlib import pyplot as plt + +logger = logging.getLogger(__name__) + + +units = { + "lj": { + "Distance": "Reduced units", + "Energy": "Reduced units", + "Force": "Reduced units", + "kb": 1, + }, + "real": { + "Distance": "[A]", + "Energy": "[kcal/mol]", + "Force": "[kcal/mol/A]", + "kb": 0.001985875, + }, + "metal": { + "Distance": "[A]", + "Energy": "[eV]", + "Force": "[eV/A]", + "kb": 8.6173332e-5, + }, + "si": {"Distance": "[m]", "Energy": "[J]", "Force": "[N]", "kb": 1.380649e-23}, + "cgs": { + "Distance": "[cm]", + "Energy": "[ergs]", + "Force": "[dynes]", + "kb": 1.3806504e-16, + }, + "electron": { + "Distance": "[Bohr]", + "Energy": "[Hartrees]", + "Force": "[Hartree/Bohr]", + "kb": 3.16681534e-6, + }, + "micro": { + "Distance": "[um]", + "Energy": "[pg·um^2/us^2]", + "Force": "[pg·um/us^2]", + "kb": 1.3806504e-8, + }, + "nano": { + "Distance": "[nm]", + "Energy": "[ag·nm^2/ns^2]", + "Force": "[ag·nm/ns^2]", + "kb": 0.013806504, + }, +} + + +def compute_energy(tp): + r = tp[0] + fo = tp[2] + e = np.zeros(r.shape) + for i, (ri, fi) in enumerate(zip(r, fo)): + if i == 0: + continue + dr = ri - r[i - 1] + e[i] = e[i - 1] - dr * fo[i - 1] + e -= e[-1] + return e + + +def main(): + + parser = argparse.ArgumentParser( + description=""" + Plots LAMMPS tabulated forces. This script takes a table + file as an input and plots all the tabulated forces inside with their + corresponding energy. The forces label is the token used to name the + force in the file. It can be used to output all the forces in separate + files and/or recompute the energy from forces through finite difference + (assuming e(rc)=0). This script requires the matplotlib and numpy + Python libraries. Bitmap format is not supported. + """ + ) + parser.add_argument( + "-u", + "--units", + dest="units", + default="real", + help="Units of the file (LAMMPS units system)", + ) + parser.add_argument( + "-f", + "--file", + dest="infile", + default="", + help="File to read", + ) + parser.add_argument( + "-x", + dest="xrange", + default="", + help="xrange separated by : (for negative values use the '=' sign: -x=-3:10)", + ) + parser.add_argument( + "-y", + dest="yrange", + default="", + help="yrange separated by :", + ) + parser.add_argument( + "-t", + dest="temp", + default=None, + type=float, + help="temperature for KbT plot [default none]", + ) + parser.add_argument( + "-d", + "--diff-num", + dest="recompute", + action="store_true", + help="Recompute the energies from forces and distances through finite differences", + ) + parser.add_argument( + "-e", + dest="extract", + action="store_true", + help="Extract the forces in separate files", + ) + args = parser.parse_args() + logging.basicConfig(level=logging.INFO) + + ########## + # Manage arguments + + udistance = units[args.units]["Distance"] + uenergy = units[args.units]["Energy"] + uforce = units[args.units]["Force"] + kb = units[args.units]["kb"] + rlabel = " ".join(["Rij", udistance]) + elabel = " ".join(["E", uenergy]) + flabel = " ".join(["F", uforce]) + etitle = "Energy" + ftitle = "Force" + font = "DejaVu Sans" + fontsize = 30 + + infile = args.infile + if not os.path.isfile(infile): + logger.error("Input file not found") + sys.exit(1) + + toplot = [] + with open(infile, "r") as f: + lines = iter(f.readlines()) + while True: + try: + r = [] + force = [] + ener = [] + tok = [] + while not tok: + tok = next(lines).partition("#")[0].rstrip() + logger.info("Found {} token".format(tok)) + infos = next(lines).split() + npoints = int(infos[1]) + next(lines) + if "bitmap" in infos: + logger.info("Unsupported bitmap format for token {:s}".format(tok)) + for _ in range(npoints): + continue + else: + for i in range(npoints): + line = next(lines).split() + r.append(float(line[1])) + ener.append(float(line[2])) + force.append(float(line[3])) + r = np.array(r) + ener = np.array(ener) + force = np.array(force) + toplot.append([r, ener, force, tok]) + tok = [] + next(lines) + except StopIteration: + break + if args.recompute: + etitle = "Estimated energy" + for tp in toplot: + tp[1] = compute_energy(tp) + + fig, axes = plt.subplots(1, 2) + + for tp in toplot: + axes[0].plot(tp[0], tp[1], label=tp[3], linewidth=3) + axes[1].plot(tp[0], tp[2], label=tp[3], linewidth=3) + hmin, hmax = axes[1].get_xlim() + axes[1].hlines(0, hmin, hmax, color="black", linewidth=3, linestyles="dashdot") + + if args.temp: + if args.temp > 0: + hmin, hmax = axes[0].get_xlim() + axes[0].hlines( + kb * args.temp, + hmin, + hmax, + color="orange", + label=r"$k_BT$", + linewidth=3, + linestyles="dashdot", + ) + axes[0].text(hmax / 2.0, kb * args.temp, "KbT", fontsize=0.7 * fontsize) + logger.info("KbT value= {:e} {:s}".format(kb * args.temp, uenergy)) + else: + logger.info("Invalid temperature value: {:e}".format(args.temp)) + + if args.xrange: + xmin, xmax = list(map(float, args.xrange.split(":"))) + axes[0].set_xlim(xmin, xmax) + axes[1].set_xlim(xmin, xmax) + if args.yrange: + ymin, ymax = list(map(float, args.yrange.split(":"))) + axes[0].set_ylim(ymin, ymax) + axes[1].set_ylim(ymin, ymax) + + # Setting axes 0 + axes[0].set_title(etitle, fontsize=fontsize) + axes[0].set_xlabel( + rlabel, fontname=font, fontsize=fontsize + ) # xlabel name, size 30pts + axes[0].set_ylabel( + elabel, fontname=font, fontsize=fontsize + ) # ylabel name, size 30pts + axes[0].tick_params( + axis="both", which="major", labelsize=fontsize + ) # Biggers ticks, bigger tick labels! + + # Setting axes 1 + axes[1].set_title(ftitle, fontsize=fontsize) + axes[1].legend(frameon=False, fontsize=fontsize) # Fat font, no frame + axes[1].set_xlabel( + rlabel, fontname=font, fontsize=fontsize + ) + axes[1].set_ylabel( + flabel, fontname=font, fontsize=fontsize + ) + axes[1].tick_params( + axis="both", which="major", labelsize=fontsize + ) + + figManager = plt.get_current_fig_manager() + figManager.window.showMaximized() + plt.show() + + if args.extract: + for tp in toplot: + outfile = "".join([tp[3], ".plot"]) + logger.info("Writing file {}".format(outfile)) + with open(outfile, "w") as f: + f.write("# {} force extracted from {}\n".format(tp[3], infile)) + f.write("# {:^20} {:^20} {:^20}\n".format("r", "energy", "force")) + for a, b, c in zip(tp[0], tp[1], tp[2]): + f.write("{:>18.16e} {:>18.16e} {:>18.16e}\n".format(a, b, c)) + return + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + raise SystemExit("User interruption.")