Merge branch 'develop' into collected-small-fixes
This commit is contained in:
Axel Kohlmeyer
@ -1250,10 +1250,10 @@ tabulate tool
|
|||||||
|
|
||||||
.. versionadded:: 22Dec2022
|
.. versionadded:: 22Dec2022
|
||||||
|
|
||||||
The ``tabulate`` folder contains Python scripts scripts to generate tabulated
|
The ``tabulate`` folder contains Python scripts scripts to generate and
|
||||||
potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module
|
visualize tabulated potential files for LAMMPS. The bulk of the code is in the
|
||||||
in the ``tabulate.py`` file. Some example files demonstrating its use are
|
``tabulate`` module in the ``tabulate.py`` file. Some example files
|
||||||
included. See the README file for more information.
|
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
|
valgrind tool
|
||||||
-------------
|
-------------
|
||||||
|
|
||||||
The ``valgrind`` folder contains additional suppressions fur LAMMPS when
|
The ``valgrind`` folder contains additional suppressions for LAMMPS when
|
||||||
using `valgrind's <https://valgrind.org/>`_ ` `memcheck tool
|
using `valgrind's <https://valgrind.org/>`_ ` `memcheck tool
|
||||||
<https://valgrind.org/info/tools.html#memcheck>`_ to search for memory
|
<https://valgrind.org/info/tools.html#memcheck>`_ to search for memory
|
||||||
access violation and memory leaks. These suppressions are automatically
|
access violation and memory leaks. These suppressions are automatically
|
||||||
|
|||||||
@ -215,6 +215,9 @@ for an overview of LAMMPS output options.
|
|||||||
The vector or array will be floating point values that correspond to
|
The vector or array will be floating point values that correspond to
|
||||||
the specified attribute.
|
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
|
The single() function of this bond style returns 0.0 for the energy
|
||||||
of a bonded interaction, since energy is not conserved in these
|
of a bonded interaction, since energy is not conserved in these
|
||||||
dissipative potentials. It also returns only the normal component of
|
dissipative potentials. It also returns only the normal component of
|
||||||
|
|||||||
@ -215,6 +215,9 @@ for an overview of LAMMPS output options.
|
|||||||
The vector or array will be floating point values that correspond to
|
The vector or array will be floating point values that correspond to
|
||||||
the specified attribute.
|
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
|
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,
|
: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()
|
ignoring any volumetric/smoothing factors or dissipative forces. The single()
|
||||||
|
|||||||
@ -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.
|
a value of 0, or a solid state, a value of 1.
|
||||||
|
|
||||||
The *surface* property indicates the surface designation produced by
|
The *surface* property indicates the surface designation produced by
|
||||||
the *interface/reconstruct* option of :doc:`fix rheo <fix_rheo>`. Bulk
|
the *surface/detection* option of :doc:`fix rheo <fix_rheo>`. Bulk
|
||||||
particles have a value of 0, surface particles have a value of 1, and
|
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
|
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
|
distance from the surface, up to the kernel cutoff length. Surface particles
|
||||||
|
|||||||
@ -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
|
|
||||||
@ -1,24 +1,43 @@
|
|||||||
units si
|
units si
|
||||||
atom_style sphere
|
atom_style sphere
|
||||||
|
comm_modify vel yes
|
||||||
|
boundary p p p
|
||||||
|
|
||||||
boundary p p f
|
region box block -0.01 0.01 -0.01 0.01 0 0.08
|
||||||
region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
|
|
||||||
create_box 2 box
|
create_box 2 box
|
||||||
|
|
||||||
read_data data.particles add append
|
create_atoms 1 single 0.0 0.0 0.04
|
||||||
group mb type 1
|
create_atoms 1 single 0.0 0.0 0.04416
|
||||||
|
set group all diameter 0.004 density 2500
|
||||||
|
|
||||||
pair_style granular
|
pair_style granular
|
||||||
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
|
pair_coeff * * hertz/material 1e6 0.8 0.4 &
|
||||||
# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution
|
tangential mindlin NULL 0.0 0.5 &
|
||||||
comm_modify vel yes
|
damping coeff_restitution
|
||||||
|
#pair_coeff * * hooke 1e6 0.5 &
|
||||||
|
# tangential mindlin 1 1.0 0.0 &
|
||||||
|
# damping coeff_restitution
|
||||||
|
|
||||||
timestep 1e-9
|
timestep 1e-9
|
||||||
fix 1 all nve/sphere
|
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
|
group a1 id 1
|
||||||
#dump_modify 1 pad 8
|
group a2 id 2
|
||||||
thermo_style custom step ke
|
|
||||||
run_style verlet
|
velocity a1 set 0 0 1.0
|
||||||
run 10000000
|
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
|
||||||
|
|||||||
@ -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
|
|
||||||
195
examples/granular/log.4Feb25.restitution.g++.1
Normal file
195
examples/granular/log.4Feb25.restitution.g++.1
Normal file
@ -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
|
||||||
@ -397,6 +397,7 @@ double BondBPM::equilibrium_distance(int /*i*/)
|
|||||||
void BondBPM::write_restart(FILE *fp)
|
void BondBPM::write_restart(FILE *fp)
|
||||||
{
|
{
|
||||||
fwrite(&overlay_flag, sizeof(int), 1, 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)
|
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(&overlay_flag, 1, MPI_INT, 0, world);
|
||||||
|
MPI_Bcast(&break_flag, 1, MPI_INT, 0, world);
|
||||||
}
|
}
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|||||||
@ -28,6 +28,7 @@
|
|||||||
#include "memory.h"
|
#include "memory.h"
|
||||||
#include "modify.h"
|
#include "modify.h"
|
||||||
#include "neighbor.h"
|
#include "neighbor.h"
|
||||||
|
#include "update.h"
|
||||||
|
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
@ -483,6 +484,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
|
|||||||
int newton_bond = force->newton_bond;
|
int newton_bond = force->newton_bond;
|
||||||
|
|
||||||
double **bondstore = fix_bond_history->bondstore;
|
double **bondstore = fix_bond_history->bondstore;
|
||||||
|
const bool allow_breaks = (update->setupflag == 0) && break_flag;
|
||||||
|
|
||||||
for (n = 0; n < nbondlist; n++) {
|
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,
|
breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
|
||||||
torque1on2, torque2on1);
|
torque1on2, torque2on1);
|
||||||
|
|
||||||
if ((breaking >= 1.0) && break_flag) {
|
if ((breaking >= 1.0) && allow_breaks) {
|
||||||
bondlist[n][2] = 0;
|
bondlist[n][2] = 0;
|
||||||
process_broken(i1, i2);
|
process_broken(i1, i2);
|
||||||
continue;
|
continue;
|
||||||
@ -763,6 +765,7 @@ void BondBPMRotational::read_restart(FILE *fp)
|
|||||||
void BondBPMRotational::write_restart_settings(FILE *fp)
|
void BondBPMRotational::write_restart_settings(FILE *fp)
|
||||||
{
|
{
|
||||||
fwrite(&smooth_flag, sizeof(int), 1, 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)
|
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(&smooth_flag, 1, MPI_INT, 0, world);
|
||||||
|
MPI_Bcast(&normalize_flag, 1, MPI_INT, 0, world);
|
||||||
}
|
}
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|||||||
@ -26,6 +26,7 @@
|
|||||||
#include "memory.h"
|
#include "memory.h"
|
||||||
#include "modify.h"
|
#include "modify.h"
|
||||||
#include "neighbor.h"
|
#include "neighbor.h"
|
||||||
|
#include "update.h"
|
||||||
|
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
@ -218,6 +219,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
|
|||||||
double invdim = 1.0 / dim;
|
double invdim = 1.0 / dim;
|
||||||
|
|
||||||
double **bondstore = fix_bond_history->bondstore;
|
double **bondstore = fix_bond_history->bondstore;
|
||||||
|
const bool allow_breaks = (update->setupflag == 0) && break_flag;
|
||||||
|
|
||||||
for (n = 0; n < nbondlist; n++) {
|
for (n = 0; n < nbondlist; n++) {
|
||||||
|
|
||||||
@ -249,7 +251,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
|
|||||||
r = sqrt(rsq);
|
r = sqrt(rsq);
|
||||||
e = (r - r0) / r0;
|
e = (r - r0) / r0;
|
||||||
|
|
||||||
if ((fabs(e) > ecrit[type]) && break_flag) {
|
if ((fabs(e) > ecrit[type]) && allow_breaks) {
|
||||||
bondlist[n][2] = 0;
|
bondlist[n][2] = 0;
|
||||||
process_broken(i1, i2);
|
process_broken(i1, i2);
|
||||||
|
|
||||||
|
|||||||
@ -26,6 +26,7 @@
|
|||||||
#include "force.h"
|
#include "force.h"
|
||||||
#include "memory.h"
|
#include "memory.h"
|
||||||
#include "neighbor.h"
|
#include "neighbor.h"
|
||||||
|
#include "update.h"
|
||||||
|
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
@ -185,6 +186,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag)
|
|||||||
int newton_bond = force->newton_bond;
|
int newton_bond = force->newton_bond;
|
||||||
|
|
||||||
double **bondstore = fix_bond_history->bondstore;
|
double **bondstore = fix_bond_history->bondstore;
|
||||||
|
const bool allow_breaks = (update->setupflag == 0) && break_flag;
|
||||||
|
|
||||||
for (n = 0; n < nbondlist; n++) {
|
for (n = 0; n < nbondlist; n++) {
|
||||||
|
|
||||||
@ -217,7 +219,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag)
|
|||||||
r = sqrt(rsq);
|
r = sqrt(rsq);
|
||||||
e = (r - r0) / r0;
|
e = (r - r0) / r0;
|
||||||
|
|
||||||
if ((fabs(e) > ecrit[type]) && break_flag) {
|
if ((fabs(e) > ecrit[type]) && allow_breaks) {
|
||||||
bondlist[n][2] = 0;
|
bondlist[n][2] = 0;
|
||||||
process_broken(i1, i2);
|
process_broken(i1, i2);
|
||||||
continue;
|
continue;
|
||||||
|
|||||||
@ -26,6 +26,7 @@
|
|||||||
#include "neigh_list.h"
|
#include "neigh_list.h"
|
||||||
#include "neighbor.h"
|
#include "neighbor.h"
|
||||||
|
|
||||||
|
#include <set>
|
||||||
#include <utility>
|
#include <utility>
|
||||||
|
|
||||||
using namespace LAMMPS_NS;
|
using namespace LAMMPS_NS;
|
||||||
@ -267,18 +268,40 @@ void FixUpdateSpecialBonds::post_run()
|
|||||||
|
|
||||||
void FixUpdateSpecialBonds::add_broken_bond(int i, int j)
|
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);
|
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)
|
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);
|
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);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/* ----------------------------------------------------------------------
|
/* ----------------------------------------------------------------------
|
||||||
|
|||||||
@ -22,6 +22,7 @@ FixStyle(UPDATE_SPECIAL_BONDS,FixUpdateSpecialBonds);
|
|||||||
|
|
||||||
#include "fix.h"
|
#include "fix.h"
|
||||||
|
|
||||||
|
#include <set>
|
||||||
#include <utility>
|
#include <utility>
|
||||||
|
|
||||||
namespace LAMMPS_NS {
|
namespace LAMMPS_NS {
|
||||||
@ -39,15 +40,15 @@ class FixUpdateSpecialBonds : public Fix {
|
|||||||
void write_restart(FILE *) override;
|
void write_restart(FILE *) override;
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
// Create two arrays to store bonds broken this timestep (new)
|
// Create array to store bonds broken this timestep (new)
|
||||||
// and since the last neighbor list build
|
// and a set for those broken since the last neighbor list build
|
||||||
std::vector<std::pair<tagint, tagint>> new_broken_pairs;
|
std::vector<std::pair<tagint, tagint>> new_broken_pairs;
|
||||||
std::vector<std::pair<tagint, tagint>> broken_pairs;
|
std::set<std::pair<tagint, tagint>> broken_pairs;
|
||||||
|
|
||||||
// Create two arrays to store newly created this timestep (new)
|
// Create arrays to store newly created this timestep (new)
|
||||||
// and since the last neighbor list build
|
// and a set for those created since the last neighbor list build
|
||||||
std::vector<std::pair<tagint, tagint>> new_created_pairs;
|
std::vector<std::pair<tagint, tagint>> new_created_pairs;
|
||||||
std::vector<std::pair<tagint, tagint>> created_pairs;
|
std::set<std::pair<tagint, tagint>> created_pairs;
|
||||||
};
|
};
|
||||||
|
|
||||||
} // namespace LAMMPS_NS
|
} // namespace LAMMPS_NS
|
||||||
|
|||||||
@ -44,7 +44,7 @@ using namespace RHEO_NS;
|
|||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
|
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)
|
id_fix(nullptr), compute_surface(nullptr)
|
||||||
{
|
{
|
||||||
partial_flag = 1;
|
partial_flag = 1;
|
||||||
@ -71,7 +71,6 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
|
|||||||
modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix));
|
modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix));
|
||||||
index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
|
index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
|
||||||
}
|
}
|
||||||
nbond = atom->ivector[index_nb];
|
|
||||||
|
|
||||||
//Store non-persistent per atom quantities, intermediate
|
//Store non-persistent per atom quantities, intermediate
|
||||||
|
|
||||||
@ -181,6 +180,7 @@ void BondRHEOShell::compute(int eflag, int vflag)
|
|||||||
double **v = atom->v;
|
double **v = atom->v;
|
||||||
double **f = atom->f;
|
double **f = atom->f;
|
||||||
tagint *tag = atom->tag;
|
tagint *tag = atom->tag;
|
||||||
|
int *nbond = atom->ivector[index_nb];
|
||||||
int *status = atom->rheo_status;
|
int *status = atom->rheo_status;
|
||||||
int **bondlist = neighbor->bondlist;
|
int **bondlist = neighbor->bondlist;
|
||||||
int nbondlist = neighbor->nbondlist;
|
int nbondlist = neighbor->nbondlist;
|
||||||
|
|||||||
@ -45,7 +45,7 @@ class BondRHEOShell : public BondBPM {
|
|||||||
double *k, *ecrit, *gamma;
|
double *k, *ecrit, *gamma;
|
||||||
double tform, rmax, rsurf;
|
double tform, rmax, rsurf;
|
||||||
|
|
||||||
int *dbond, *nbond;
|
int *dbond;
|
||||||
int index_nb, nmax_store;
|
int index_nb, nmax_store;
|
||||||
char *id_fix;
|
char *id_fix;
|
||||||
|
|
||||||
|
|||||||
@ -39,12 +39,10 @@ using namespace LAMMPS_NS;
|
|||||||
using namespace RHEO_NS;
|
using namespace RHEO_NS;
|
||||||
using namespace MathExtra;
|
using namespace MathExtra;
|
||||||
|
|
||||||
static constexpr double EPSILON = 1e-1;
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
|
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),
|
norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), compute_kernel(nullptr),
|
||||||
fix_pressure(nullptr)
|
fix_pressure(nullptr)
|
||||||
|
|
||||||
@ -64,13 +62,12 @@ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
|
|||||||
// between timesteps (fix property atom will handle callbacks)
|
// between timesteps (fix property atom will handle callbacks)
|
||||||
|
|
||||||
int tmp1, tmp2;
|
int tmp1, tmp2;
|
||||||
int index = atom->find_custom("fp_store", tmp1, tmp2);
|
index_fp_store = atom->find_custom("fp_store", tmp1, tmp2);
|
||||||
if (index == -1) {
|
if (index_fp_store == -1) {
|
||||||
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
|
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));
|
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()
|
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];
|
double xtmp, ytmp, ztmp, rsq, w, dot, dx[3];
|
||||||
|
|
||||||
int inum, *ilist, *jlist, *numneigh, **firstneigh;
|
int inum, *ilist, *jlist, *numneigh, **firstneigh;
|
||||||
@ -126,6 +123,7 @@ void ComputeRHEOInterface::compute_peratom()
|
|||||||
int newton = force->newton;
|
int newton = force->newton;
|
||||||
int *status = atom->rheo_status;
|
int *status = atom->rheo_status;
|
||||||
double *rho = atom->rho;
|
double *rho = atom->rho;
|
||||||
|
double **fp_store = atom->darray[index_fp_store];
|
||||||
|
|
||||||
inum = list->inum;
|
inum = list->inum;
|
||||||
ilist = list->ilist;
|
ilist = list->ilist;
|
||||||
@ -151,7 +149,6 @@ void ComputeRHEOInterface::compute_peratom()
|
|||||||
xtmp = x[i][0];
|
xtmp = x[i][0];
|
||||||
ytmp = x[i][1];
|
ytmp = x[i][1];
|
||||||
ztmp = x[i][2];
|
ztmp = x[i][2];
|
||||||
itype = type[i];
|
|
||||||
fluidi = !(status[i] & PHASECHECK);
|
fluidi = !(status[i] & PHASECHECK);
|
||||||
jlist = firstneigh[i];
|
jlist = firstneigh[i];
|
||||||
jnum = numneigh[i];
|
jnum = numneigh[i];
|
||||||
@ -212,9 +209,10 @@ void ComputeRHEOInterface::compute_peratom()
|
|||||||
if (status[i] & PHASECHECK) {
|
if (status[i] & PHASECHECK) {
|
||||||
if (normwf[i] != 0.0) {
|
if (normwf[i] != 0.0) {
|
||||||
// Stores rho for solid particles 1+Pw in Adami Adams 2012
|
// 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 {
|
} 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;
|
int m = 0;
|
||||||
double *rho = atom->rho;
|
double *rho = atom->rho;
|
||||||
|
double **fp_store = atom->darray[index_fp_store];
|
||||||
|
|
||||||
for (int i = 0; i < n; i++) {
|
for (int i = 0; i < n; i++) {
|
||||||
int j = list[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)
|
void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
|
||||||
{
|
{
|
||||||
double *rho = atom->rho;
|
double *rho = atom->rho;
|
||||||
|
double **fp_store = atom->darray[index_fp_store];
|
||||||
|
|
||||||
int m = 0;
|
int m = 0;
|
||||||
int last = first + n;
|
int last = first + n;
|
||||||
for (int i = first; i < last; i++) {
|
for (int i = first; i < last; i++) {
|
||||||
@ -335,6 +336,7 @@ void ComputeRHEOInterface::store_forces()
|
|||||||
double *mass = atom->mass;
|
double *mass = atom->mass;
|
||||||
double *rmass = atom->rmass;
|
double *rmass = atom->rmass;
|
||||||
double **f = atom->f;
|
double **f = atom->f;
|
||||||
|
double **fp_store = atom->darray[index_fp_store];
|
||||||
|
|
||||||
// When this is called, fp_store stores the pressure force
|
// When this is called, fp_store stores the pressure force
|
||||||
// After this method, fp_store instead stores non-pressure forces
|
// After this method, fp_store instead stores non-pressure forces
|
||||||
|
|||||||
@ -40,7 +40,8 @@ class ComputeRHEOInterface : public Compute {
|
|||||||
double correct_rho(int);
|
double correct_rho(int);
|
||||||
void store_forces();
|
void store_forces();
|
||||||
|
|
||||||
double *chi, **fp_store;
|
double *chi;
|
||||||
|
int index_fp_store;
|
||||||
class FixRHEO *fix_rheo;
|
class FixRHEO *fix_rheo;
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
|||||||
@ -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);
|
return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji);
|
||||||
if (kernel_style == QUINTIC)
|
if (kernel_style == QUINTIC)
|
||||||
return calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
|
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;
|
double wp;
|
||||||
int corrections_i = check_corrections(i);
|
int corrections_i = check_corrections(i);
|
||||||
|
|||||||
@ -502,7 +502,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n)
|
|||||||
|
|
||||||
void ComputeRHEOPropertyAtom::pack_nbond_shell(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 *mask = atom->mask;
|
||||||
int nlocal = atom->nlocal;
|
int nlocal = atom->nlocal;
|
||||||
|
|
||||||
|
|||||||
@ -39,6 +39,7 @@ using namespace RHEO_NS;
|
|||||||
using namespace FixConst;
|
using namespace FixConst;
|
||||||
|
|
||||||
static const char cite_rheo[] =
|
static const char cite_rheo[] =
|
||||||
|
"RHEO package: doi:10.1063/5.0228823\n\n"
|
||||||
"@article{Palermo2024,\n"
|
"@article{Palermo2024,\n"
|
||||||
" journal = {Physics of Fluids},\n"
|
" journal = {Physics of Fluids},\n"
|
||||||
" title = {Reproducing hydrodynamics and elastic objects: A hybrid mesh-free model framework for dynamic multi-phase flows},\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;
|
shift_flag = 1;
|
||||||
memory->create(shift_type, n + 1, "rheo:shift_type");
|
memory->create(shift_type, n + 1, "rheo:shift_type");
|
||||||
for (i = 1; i <= n; i++) shift_type[i] = 1;
|
for (i = 1; i <= n; i++) shift_type[i] = 1;
|
||||||
while (iarg < narg) { // optional sub-arguments
|
while (iarg + 1 < narg) { // optional sub-arguments
|
||||||
if (strcmp(arg[iarg], "scale/cross/type") == 0) {
|
if (strcmp(arg[iarg + 1], "scale/cross/type") == 0) {
|
||||||
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
|
if (iarg + 4 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
|
||||||
shift_cross_type_flag = 1;
|
shift_cross_type_flag = 1;
|
||||||
shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
|
shift_scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
|
||||||
shift_cmin = 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 + 3], false, lmp);
|
shift_wmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
|
||||||
iarg += 3;
|
iarg += 3;
|
||||||
} else if (strcmp(arg[iarg], "exclude/type") == 0) {
|
} else if (strcmp(arg[iarg + 1], "exclude/type") == 0) {
|
||||||
if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
|
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
|
||||||
utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error);
|
utils::bounds(FLERR, arg[iarg + 2], 1, n, nlo, nhi, error);
|
||||||
for (i = nlo; i <= nhi; i++) shift_type[i] = 0;
|
for (i = nlo; i <= nhi; i++) shift_type[i] = 0;
|
||||||
iarg += 1;
|
iarg += 1;
|
||||||
} else {
|
} else {
|
||||||
|
|||||||
@ -38,6 +38,7 @@ using namespace FixConst;
|
|||||||
enum { NONE, CONSTANT };
|
enum { NONE, CONSTANT };
|
||||||
|
|
||||||
static const char cite_rheo_oxide[] =
|
static const char cite_rheo_oxide[] =
|
||||||
|
"RHEO oxidation: doi:10.1016/j.apm.2024.02.027\n\n"
|
||||||
"@article{ApplMathModel.130.310,\n"
|
"@article{ApplMathModel.130.310,\n"
|
||||||
" title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
|
" title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
|
||||||
" journal = {Applied Mathematical Modelling},\n"
|
" journal = {Applied Mathematical Modelling},\n"
|
||||||
@ -109,7 +110,6 @@ void FixRHEOOxidation::init()
|
|||||||
int tmp1, tmp2;
|
int tmp1, tmp2;
|
||||||
index_nb = atom->find_custom("shell_nbond", 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");
|
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
|
// need a half neighbor list
|
||||||
auto req = neighbor->add_request(this, NeighConst::REQ_FULL);
|
auto req = neighbor->add_request(this, NeighConst::REQ_FULL);
|
||||||
|
|||||||
@ -39,11 +39,11 @@ class FixRHEOOxidation : public Fix {
|
|||||||
void post_force(int) override;
|
void post_force(int) override;
|
||||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||||
void unpack_forward_comm(int, int, double *) override;
|
void unpack_forward_comm(int, int, double *) override;
|
||||||
int *nbond;
|
|
||||||
double rsurf, cut;
|
double rsurf, cut;
|
||||||
|
int index_nb;
|
||||||
|
|
||||||
private:
|
private:
|
||||||
int btype, index_nb;
|
int btype;
|
||||||
double cutsq;
|
double cutsq;
|
||||||
|
|
||||||
class NeighList *list;
|
class NeighList *list;
|
||||||
|
|||||||
@ -221,11 +221,6 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
|
|||||||
|
|
||||||
FixRHEOThermal::~FixRHEOThermal()
|
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(cv_style);
|
||||||
memory->destroy(Tc_style);
|
memory->destroy(Tc_style);
|
||||||
memory->destroy(kappa_style);
|
memory->destroy(kappa_style);
|
||||||
|
|||||||
@ -115,7 +115,7 @@ void PairRHEO::compute(int eflag, int vflag)
|
|||||||
|
|
||||||
double **fp_store, *chi;
|
double **fp_store, *chi;
|
||||||
if (compute_interface) {
|
if (compute_interface) {
|
||||||
fp_store = compute_interface->fp_store;
|
fp_store = atom->darray[compute_interface->index_fp_store];
|
||||||
chi = compute_interface->chi;
|
chi = compute_interface->chi;
|
||||||
|
|
||||||
for (i = 0; i < atom->nmax; i++) {
|
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)
|
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 m = 0;
|
||||||
int last = first + n;
|
int last = first + n;
|
||||||
for (int i = first; i < last; i++) {
|
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)
|
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;
|
int m = 0;
|
||||||
for (int i = 0; i < n; i++) {
|
for (int i = 0; i < n; i++) {
|
||||||
int j = list[i];
|
int j = list[i];
|
||||||
|
|||||||
@ -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_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 |
|
| 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 |
|
| 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:
|
Common command line flags:
|
||||||
|
|
||||||
|
|||||||
278
tools/tabulate/plot_forces.py
Executable file
278
tools/tabulate/plot_forces.py
Executable file
@ -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.")
|
||||||
Reference in New Issue
Block a user