Merge pull request #1092 from RomainVermorel/master
New stress/mop and stress/mop/profile computes for USER-MISC
This commit is contained in:
Axel Kohlmeyer
111
doc/src/compute_stress_mop.txt
Normal file
111
doc/src/compute_stress_mop.txt
Normal file
@ -0,0 +1,111 @@
|
||||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
compute stress/mop command :h3
|
||||
compute stress/mop/profile command :h3
|
||||
|
||||
|
||||
[Syntax:]
|
||||
|
||||
compute ID group-ID style dir args keywords ... :pre
|
||||
|
||||
ID, group-ID are documented in "compute"_compute.html command
|
||||
style = {stress/mop} or {stress/mop/profile}
|
||||
dir = {x} or {y} or {z} is the direction normal to the plane
|
||||
args = argument specific to the compute style
|
||||
keywords = {kin} or {conf} or {total} (one of more can be specified) :ul
|
||||
{stress/mop} args = pos
|
||||
pos = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the plane
|
||||
{stress/mop/profile} args = origin delta
|
||||
origin = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the first plane
|
||||
delta = value (distance units) is the distance between planes :pre
|
||||
|
||||
compute 1 all stress/mop x lower total
|
||||
compute 1 liquid stress/mop z 0.0 kin conf
|
||||
fix 1 all ave/time 10 1000 10000 c_1\[*\] file mop.time
|
||||
fix 1 all ave/time 10 1000 10000 c_1\[2\] file mop.time :pre
|
||||
|
||||
compute 1 all stress/mop/profile x lower 0.1 total
|
||||
compute 1 liquid stress/mop/profile z 0.0 0.25 kin conf
|
||||
fix 1 all ave/time 500 20 10000 c_1\[*\] ave running overwrite file mopp.time mode vector :pre
|
||||
|
||||
|
||||
[Description:]
|
||||
|
||||
Compute {stress/mop} and compute {stress/mop/profile} define computations that
|
||||
calculate components of the local stress tensor using the method of
|
||||
planes "(Todd)"_#mop-todd. Specifically in compute {stress/mop} calculates 3
|
||||
components are computed in directions {dir},{x}; {dir},{y}; and
|
||||
{dir},{z}; where {dir} is the direction normal to the plane, while
|
||||
in compute {stress/mop/profile} the profile of the stress is computed.
|
||||
|
||||
Contrary to methods based on histograms of atomic stress (i.e. using
|
||||
"compute stress/atom"_compute_stress_atom.html), the method of planes is
|
||||
compatible with mechanical balance in heterogeneous systems and at
|
||||
interfaces "(Todd)"_#mop-todd.
|
||||
|
||||
The stress tensor is the sum of a kinetic term and a configurational
|
||||
term, which are given respectively by Eq. (21) and Eq. (16) in
|
||||
"(Todd)"_#mop-todd. For the kinetic part, the algorithm considers that
|
||||
atoms have crossed the plane if their positions at times t-dt and t are
|
||||
one on either side of the plane, and uses the velocity at time t-dt/2
|
||||
given by the velocity-Verlet algorithm.
|
||||
|
||||
Between one and three keywords can be used to indicate which
|
||||
contributions to the stress must be computed: kinetic stress (kin),
|
||||
configurational stress (conf), and/or total stress (total).
|
||||
|
||||
NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID.
|
||||
|
||||
NOTE 2: The local stress does not include any Lennard-Jones tail
|
||||
corrections to the pressure added by the "pair_modify tail
|
||||
yes"_pair_modify.html command, since those are contributions to the global system pressure.
|
||||
|
||||
[Output info:]
|
||||
|
||||
Compute {stress/mop} calculates a global vector (indices starting at 1), with 3
|
||||
values for each declared keyword (in the order the keywords have been
|
||||
declared). For each keyword, the stress tensor components are ordered as
|
||||
follows: stress_dir,x, stress_dir,y, and stress_dir,z.
|
||||
|
||||
Compute {stress/mop/profile} instead calculates a global array, with 1 column
|
||||
giving the position of the planes where the stress tensor was computed,
|
||||
and with 3 columns of values for each declared keyword (in the order the
|
||||
keywords have been declared). For each keyword, the profiles of stress
|
||||
tensor components are ordered as follows: stress_dir,x; stress_dir,y;
|
||||
and stress_dir,z.
|
||||
|
||||
The values are in pressure "units"_units.html.
|
||||
|
||||
The values produced by this compute can be accessed by various "output commands"_Howto_output.html. For instance, the results can be written to a file using the "fix ave/time"_fix_ave_time.html command. Please see the example in the examples/USER/mop folder.
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
These styles are part of the USER-MISC package. They are only enabled if
|
||||
LAMMPS is built with that package. See the "Build package"_Build_package.html
|
||||
doc page on for more info.
|
||||
|
||||
The method is only implemented for 3d orthogonal simulation boxes whose
|
||||
size does not change in time, and axis-aligned planes.
|
||||
|
||||
The method only works with two-body pair interactions, because it
|
||||
requires the class method pair->single() to be implemented. In
|
||||
particular, it does not work with more than two-body pair interactions,
|
||||
intra-molecular interactions, and long range (kspace) interactions.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"compute stress/atom"_compute_stress_atom.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(mop-todd)
|
||||
[(Todd)] B. D. Todd, Denis J. Evans, and Peter J. Daivis: "Pressure tensor for inhomogeneous fluids",
|
||||
Phys. Rev. E 52, 1627 (1995).
|
||||
@ -100,6 +100,7 @@ Computes :h1
|
||||
compute_sna_atom
|
||||
compute_spin
|
||||
compute_stress_atom
|
||||
compute_stress_mop
|
||||
compute_tally
|
||||
compute_tdpd_cc_atom
|
||||
compute_temp
|
||||
|
||||
@ -497,6 +497,7 @@ compute_smd_vol.html
|
||||
compute_sna_atom.html
|
||||
compute_spin.html
|
||||
compute_stress_atom.html
|
||||
compute_stress_mop.html
|
||||
compute_tally.html
|
||||
compute_tdpd_cc_atom.html
|
||||
compute_temp.html
|
||||
|
||||
2467
examples/USER/misc/mop/data.mop
Normal file
2467
examples/USER/misc/mop/data.mop
Normal file
File diff suppressed because it is too large
Load Diff
40
examples/USER/misc/mop/in.compute_stress_mop
Normal file
40
examples/USER/misc/mop/in.compute_stress_mop
Normal file
@ -0,0 +1,40 @@
|
||||
variable T equal 0.8
|
||||
variable p_solid equal 0.05
|
||||
|
||||
read_data data.mop
|
||||
|
||||
pair_style lj/cut 2.5
|
||||
pair_coeff * * 1.0 1.0
|
||||
pair_coeff 1 2 0.5 1.0
|
||||
pair_coeff 2 2 0.0 0.0
|
||||
neigh_modify delay 0
|
||||
|
||||
group liquid type 1
|
||||
group solid type 2
|
||||
region bottom block INF INF INF INF INF 7.0
|
||||
group bottom region bottom
|
||||
group solid_bottom intersect solid bottom
|
||||
group solid_up subtract solid solid_bottom
|
||||
|
||||
variable faSolid equal ${p_solid}*lx*ly/count(solid_up)
|
||||
fix piston_up solid_up aveforce NULL NULL -${faSolid}
|
||||
fix freeze_up solid_up setforce 0.0 0.0 NULL
|
||||
fix freeze_bottom solid_bottom setforce 0.0 0.0 0.0
|
||||
fix nvesol solid nve
|
||||
compute Tliq liquid temp
|
||||
fix nvtliq liquid nvt temp $T $T 0.5
|
||||
fix_modify nvtliq temp Tliq
|
||||
|
||||
thermo 1000
|
||||
thermo_modify flush yes temp Tliq
|
||||
|
||||
fix fxbal all balance 1000 1.05 shift z 10 1.05
|
||||
|
||||
compute mopz0 all stress/mop z center kin conf
|
||||
fix mopz0t all ave/time 1 1 1 c_mopz0[*] file mopz0.time
|
||||
|
||||
compute moppz liquid stress/mop/profile z 0.0 0.1 kin conf
|
||||
fix moppzt all ave/time 1 1 1 c_moppz[*] ave running overwrite file moppz.time mode vector
|
||||
|
||||
run 0
|
||||
|
||||
111
examples/USER/misc/mop/log.5Sep18.compute_stress_mop.g++.1
Normal file
111
examples/USER/misc/mop/log.5Sep18.compute_stress_mop.g++.1
Normal file
@ -0,0 +1,111 @@
|
||||
LAMMPS (5 Sep 2018)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:87)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
variable T equal 0.8
|
||||
variable p_solid equal 0.05
|
||||
|
||||
read_data data.mop
|
||||
orthogonal box = (0 0 -2) to (9.52441 9.52441 16)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
1224 atoms
|
||||
reading velocities ...
|
||||
1224 velocities
|
||||
|
||||
pair_style lj/cut 2.5
|
||||
pair_coeff * * 1.0 1.0
|
||||
pair_coeff 1 2 0.5 1.0
|
||||
pair_coeff 2 2 0.0 0.0
|
||||
neigh_modify delay 0
|
||||
|
||||
group liquid type 1
|
||||
792 atoms in group liquid
|
||||
group solid type 2
|
||||
432 atoms in group solid
|
||||
region bottom block INF INF INF INF INF 7.0
|
||||
group bottom region bottom
|
||||
630 atoms in group bottom
|
||||
group solid_bottom intersect solid bottom
|
||||
216 atoms in group solid_bottom
|
||||
group solid_up subtract solid solid_bottom
|
||||
216 atoms in group solid_up
|
||||
|
||||
variable faSolid equal ${p_solid}*lx*ly/count(solid_up)
|
||||
variable faSolid equal 0.05*lx*ly/count(solid_up)
|
||||
fix piston_up solid_up aveforce NULL NULL -${faSolid}
|
||||
fix piston_up solid_up aveforce NULL NULL -0.0209986841649146
|
||||
fix freeze_up solid_up setforce 0.0 0.0 NULL
|
||||
fix freeze_bottom solid_bottom setforce 0.0 0.0 0.0
|
||||
fix nvesol solid nve
|
||||
compute Tliq liquid temp
|
||||
fix nvtliq liquid nvt temp $T $T 0.5
|
||||
fix nvtliq liquid nvt temp 0.8 $T 0.5
|
||||
fix nvtliq liquid nvt temp 0.8 0.8 0.5
|
||||
fix_modify nvtliq temp Tliq
|
||||
WARNING: Temperature for fix modify is not for group all (src/fix_nh.cpp:1404)
|
||||
|
||||
thermo 1000
|
||||
thermo_modify flush yes temp Tliq
|
||||
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:488)
|
||||
|
||||
fix fxbal all balance 1000 1.05 shift z 10 1.05
|
||||
|
||||
compute mopz0 all stress/mop z center kin conf
|
||||
fix mopz0t all ave/time 1 1 1 c_mopz0[*] file mopz0.time
|
||||
|
||||
compute moppz liquid stress/mop/profile z 0.0 0.1 kin conf
|
||||
fix moppzt all ave/time 1 1 1 c_moppz[*] ave running overwrite file moppz.time mode vector
|
||||
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 0 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 2.8
|
||||
ghost atom cutoff = 2.8
|
||||
binsize = 1.4, bins = 7 7 13
|
||||
3 neighbor lists, perpetual/occasional/extra = 1 2 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
(2) compute stress/mop, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(3) compute stress/mop/profile, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.596 | 3.596 | 3.596 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 0.82011245 -3.0642111 0 -2.2692246 0.16906107 1632.8577
|
||||
Loop time of 1.19209e-06 on 1 procs for 0 steps with 1224 atoms
|
||||
|
||||
167.8% 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 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 1.192e-06 | | |100.00
|
||||
|
||||
Nlocal: 1224 ave 1224 max 1224 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 2975 ave 2975 max 2975 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 40241 ave 40241 max 40241 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 40241
|
||||
Ave neighs/atom = 32.8766
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:00:00
|
||||
111
examples/USER/misc/mop/log.5Sep18.compute_stress_mop.g++.4
Normal file
111
examples/USER/misc/mop/log.5Sep18.compute_stress_mop.g++.4
Normal file
@ -0,0 +1,111 @@
|
||||
LAMMPS (5 Sep 2018)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:87)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
variable T equal 0.8
|
||||
variable p_solid equal 0.05
|
||||
|
||||
read_data data.mop
|
||||
orthogonal box = (0 0 -2) to (9.52441 9.52441 16)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
1224 atoms
|
||||
reading velocities ...
|
||||
1224 velocities
|
||||
|
||||
pair_style lj/cut 2.5
|
||||
pair_coeff * * 1.0 1.0
|
||||
pair_coeff 1 2 0.5 1.0
|
||||
pair_coeff 2 2 0.0 0.0
|
||||
neigh_modify delay 0
|
||||
|
||||
group liquid type 1
|
||||
792 atoms in group liquid
|
||||
group solid type 2
|
||||
432 atoms in group solid
|
||||
region bottom block INF INF INF INF INF 7.0
|
||||
group bottom region bottom
|
||||
630 atoms in group bottom
|
||||
group solid_bottom intersect solid bottom
|
||||
216 atoms in group solid_bottom
|
||||
group solid_up subtract solid solid_bottom
|
||||
216 atoms in group solid_up
|
||||
|
||||
variable faSolid equal ${p_solid}*lx*ly/count(solid_up)
|
||||
variable faSolid equal 0.05*lx*ly/count(solid_up)
|
||||
fix piston_up solid_up aveforce NULL NULL -${faSolid}
|
||||
fix piston_up solid_up aveforce NULL NULL -0.0209986841649146
|
||||
fix freeze_up solid_up setforce 0.0 0.0 NULL
|
||||
fix freeze_bottom solid_bottom setforce 0.0 0.0 0.0
|
||||
fix nvesol solid nve
|
||||
compute Tliq liquid temp
|
||||
fix nvtliq liquid nvt temp $T $T 0.5
|
||||
fix nvtliq liquid nvt temp 0.8 $T 0.5
|
||||
fix nvtliq liquid nvt temp 0.8 0.8 0.5
|
||||
fix_modify nvtliq temp Tliq
|
||||
WARNING: Temperature for fix modify is not for group all (src/fix_nh.cpp:1404)
|
||||
|
||||
thermo 1000
|
||||
thermo_modify flush yes temp Tliq
|
||||
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:488)
|
||||
|
||||
fix fxbal all balance 1000 1.05 shift z 10 1.05
|
||||
|
||||
compute mopz0 all stress/mop z center kin conf
|
||||
fix mopz0t all ave/time 1 1 1 c_mopz0[*] file mopz0.time
|
||||
|
||||
compute moppz liquid stress/mop/profile z 0.0 0.1 kin conf
|
||||
fix moppzt all ave/time 1 1 1 c_moppz[*] ave running overwrite file moppz.time mode vector
|
||||
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 0 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 2.8
|
||||
ghost atom cutoff = 2.8
|
||||
binsize = 1.4, bins = 7 7 13
|
||||
3 neighbor lists, perpetual/occasional/extra = 1 2 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
(2) compute stress/mop, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(3) compute stress/mop/profile, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.509 | 3.51 | 3.511 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 0.82011245 -3.0642111 0 -2.2692246 0.16906107 1632.8577
|
||||
Loop time of 4.06504e-05 on 4 procs for 0 steps with 1224 atoms
|
||||
|
||||
65.2% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 4.065e-05 | | |100.00
|
||||
|
||||
Nlocal: 306 ave 320 max 295 min
|
||||
Histogram: 1 1 0 0 0 0 1 0 0 1
|
||||
Nghost: 1450.25 ave 1485 max 1422 min
|
||||
Histogram: 2 0 0 0 0 0 0 1 0 1
|
||||
Neighs: 10060.2 ave 10866 max 9507 min
|
||||
Histogram: 2 0 0 0 0 1 0 0 0 1
|
||||
|
||||
Total # of neighbors = 40241
|
||||
Ave neighs/atom = 32.8766
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:00:00
|
||||
185
examples/USER/misc/mop/moppz.time.reference
Normal file
185
examples/USER/misc/mop/moppz.time.reference
Normal file
@ -0,0 +1,185 @@
|
||||
# Time-averaged data for fix moppzt
|
||||
# TimeStep Number-of-rows
|
||||
# Row c_moppz[1] c_moppz[2] c_moppz[3] c_moppz[4] c_moppz[5] c_moppz[6] c_moppz[7]
|
||||
0 181
|
||||
1 -2 0 0 0 0 0 0
|
||||
2 -1.9 0 0 0 0 0 0
|
||||
3 -1.8 0 0 0 0 0 0
|
||||
4 -1.7 0 0 0 0 0 0
|
||||
5 -1.6 0 0 0 0 0 0
|
||||
6 -1.5 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
7 -1.4 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
8 -1.3 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
9 -1.2 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
10 -1.1 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
11 -1 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
12 -0.9 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
13 -0.8 0 0 0 -9.81273e-05 0.000228605 -0.00421138
|
||||
14 -0.7 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
15 -0.6 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
16 -0.5 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
17 -0.4 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
18 -0.3 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
19 -0.2 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
20 -0.1 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
21 0 0 0 0 -0.000370675 -0.00240125 -0.26848
|
||||
22 0.1 0 0 0 0.190761 -0.491728 0.287704
|
||||
23 0.2 0 0 0 0.190761 -0.491728 0.287704
|
||||
24 0.3 0 0 0 0.190761 -0.491728 0.287704
|
||||
25 0.4 0 0 0 0.190761 -0.491728 0.287704
|
||||
26 0.5 0 0 0 0.190761 -0.491728 0.287704
|
||||
27 0.6 0 0 0 0.190761 -0.491728 0.287704
|
||||
28 0.7 0 0 0 0.190761 -0.491728 0.287704
|
||||
29 0.8 0 0 0 -0.181602 -0.198457 -0.0964774
|
||||
30 0.9 0 0 0 -0.15138 0.183353 0.206848
|
||||
31 1 0 0 0 0.174362 1.27701 0.600545
|
||||
32 1.1 0 0 0 0.160987 0.563442 0.494994
|
||||
33 1.2 0 0 0 0.218876 0.59796 0.398527
|
||||
34 1.3 0 0 0 0.187614 0.558909 0.372353
|
||||
35 1.4 0 0 0 0.118586 0.410013 0.331945
|
||||
36 1.5 0 0 0 -0.0514208 0.40381 0.128097
|
||||
37 1.6 3.08628 0.241189 5.90817 -0.198262 0.324128 -0.0449302
|
||||
38 1.7 0 0 0 -0.104542 0.256677 -0.332854
|
||||
39 1.8 0.222123 2.43524 1.10089 -0.324638 -0.168682 -1.06238
|
||||
40 1.9 0 0 0 -0.175732 -0.186846 -0.163062
|
||||
41 2 0 0 0 -0.137995 0.0920401 -0.260106
|
||||
42 2.1 -0.179621 -2.59775 1.80077 -0.480624 -0.0439511 -0.0824913
|
||||
43 2.2 0 0 0 -0.499868 -0.0106185 -0.108924
|
||||
44 2.3 0 0 0 -0.703301 0.124555 -0.0880158
|
||||
45 2.4 0 0 0 -0.581211 -0.244281 -0.250071
|
||||
46 2.5 1.05274 -2.86043 3.36339 -0.575104 -0.148715 -0.249092
|
||||
47 2.6 0 0 0 0.66061 -0.157649 -0.357141
|
||||
48 2.7 0 0 0 0.299971 -0.302298 -0.572714
|
||||
49 2.8 0 0 0 0.33107 -0.201699 -0.470466
|
||||
50 2.9 0 0 0 0.822686 1.08427 -0.390511
|
||||
51 3 0 0 0 0.716428 0.750998 -0.698174
|
||||
52 3.1 0.805189 0.571878 4.31938 0.121891 0.922727 -0.932582
|
||||
53 3.2 0 0 0 0.0442642 1.02537 -1.03066
|
||||
54 3.3 2.54289 -1.93701 4.88355 0.0731321 1.09091 -0.83075
|
||||
55 3.4 0 0 0 0.426589 0.821174 -0.765855
|
||||
56 3.5 0 0 0 0.445135 0.299996 -1.48972
|
||||
57 3.6 0 0 0 0.362916 -1.28673 -0.853897
|
||||
58 3.7 0.952867 -1.07044 1.04141 0.12517 -1.00353 -0.785272
|
||||
59 3.8 0.617661 0.991499 1.80973 -0.182369 -1.04057 -1.00435
|
||||
60 3.9 0.60295 -2.41888 3.98011 0.0347345 -1.01302 -0.88314
|
||||
61 4 -2.97421 -2.01531 2.98586 0.43463 -0.465643 -0.801128
|
||||
62 4.1 -3.23318 -3.31281 0.956525 0.732752 0.140718 -1.10583
|
||||
63 4.2 0 0 0 0.969872 0.298566 -0.823464
|
||||
64 4.3 0 0 0 0.7707 0.557002 -0.836549
|
||||
65 4.4 0 0 0 0.395828 0.66755 -1.53454
|
||||
66 4.5 0 0 0 0.104451 0.46777 -1.32358
|
||||
67 4.6 0 0 0 0.402084 0.464983 -1.22051
|
||||
68 4.7 0 0 0 0.352808 0.0794986 -1.31292
|
||||
69 4.8 0 0 0 0.0215512 0.284343 -0.975326
|
||||
70 4.9 0 0 0 -0.133637 0.250925 -1.33918
|
||||
71 5 0 0 0 -0.066208 0.104514 -1.27412
|
||||
72 5.1 0 0 0 -0.184391 0.479805 -1.15139
|
||||
73 5.2 0 0 0 -0.200251 0.527142 -1.34307
|
||||
74 5.3 0 0 0 0.043532 -0.0788824 -0.998406
|
||||
75 5.4 0 0 0 -0.531846 0.126289 -1.05818
|
||||
76 5.5 0 0 0 -0.259593 0.0818463 -1.58939
|
||||
77 5.6 0 0 0 -0.373828 -0.343977 -1.50908
|
||||
78 5.7 -0.294161 -1.07567 3.46536 -0.0644873 -0.424333 -1.28548
|
||||
79 5.8 0 0 0 -0.293233 -0.201133 -1.19085
|
||||
80 5.9 0.961568 -1.44949 2.42101 -0.632816 -0.0669315 -0.85119
|
||||
81 6 0 0 0 -0.0559892 -0.0194478 -1.04541
|
||||
82 6.1 0 0 0 -0.339753 0.286693 -1.24366
|
||||
83 6.2 0 0 0 -0.376208 0.444053 -1.7662
|
||||
84 6.3 0 0 0 -0.718923 0.555398 -1.93862
|
||||
85 6.4 0 0 0 -1.10631 0.263525 -1.79723
|
||||
86 6.5 0 0 0 -0.217948 -0.0489491 -2.07833
|
||||
87 6.6 0 0 0 -0.376248 -0.0588682 -2.45322
|
||||
88 6.7 -2.12742 4.22609 2.36568 -0.236703 -0.279582 -1.56434
|
||||
89 6.8 0.869072 -0.141389 3.92123 0.0540986 -0.00271606 -0.930143
|
||||
90 6.9 0 0 0 1.08829 -1.11737 -0.808187
|
||||
91 7 1.62633 1.08234 0.844097 1.18575 -0.408792 -0.752394
|
||||
92 7.1 0 0 0 1.03324 -0.470631 -0.486767
|
||||
93 7.2 0 0 0 0.950164 -0.112451 -0.479409
|
||||
94 7.3 -2.66121 -0.326607 7.83093 0.359 -0.482493 0.154384
|
||||
95 7.4 0 0 0 0.359089 -1.12337 0.409711
|
||||
96 7.5 -1.88971 1.34806 3.56893 0.394677 -1.0109 0.548348
|
||||
97 7.6 -1.34494 -0.896214 2.06959 0.231398 -0.728529 0.313513
|
||||
98 7.7 0 0 0 0.415681 -0.45268 0.507181
|
||||
99 7.8 0 0 0 0.259423 -0.11638 0.464208
|
||||
100 7.9 -1.97572 -1.20836 3.95731 0.252257 -0.0845701 -0.249345
|
||||
101 8 0 0 0 0.0688154 0.290386 -0.462467
|
||||
102 8.1 0.25925 -0.458269 3.33086 0.360399 -0.0409494 -0.656911
|
||||
103 8.2 0 0 0 -0.0587033 0.347698 -0.340604
|
||||
104 8.3 0 0 0 -0.377192 0.153096 -0.914654
|
||||
105 8.4 0 0 0 -0.431553 0.274996 -0.946252
|
||||
106 8.5 0 0 0 -0.898366 0.146653 -1.36383
|
||||
107 8.6 0 0 0 -0.889593 0.385951 0.125116
|
||||
108 8.7 0 0 0 -0.0139171 -0.162302 -0.0287854
|
||||
109 8.8 0 0 0 -0.266284 -0.148945 0.393533
|
||||
110 8.9 0 0 0 -0.00920376 -0.0770818 0.334642
|
||||
111 9 0 0 0 -0.0949156 0.0113352 -0.0761263
|
||||
112 9.1 0 0 0 0.0688045 0.104558 -0.101891
|
||||
113 9.2 3.79773 0.0255401 3.75032 0.419832 0.295402 0.652533
|
||||
114 9.3 0 0 0 0.594267 0.70396 0.836434
|
||||
115 9.4 0 0 0 0.174722 1.00483 1.42787
|
||||
116 9.5 0 0 0 0.0626835 0.518952 0.269158
|
||||
117 9.6 0 0 0 -0.302859 -0.265212 -0.0145578
|
||||
118 9.7 0 0 0 -0.114026 -0.201336 -0.539522
|
||||
119 9.8 0 0 0 0.104008 -0.30236 -0.0789062
|
||||
120 9.9 0 0 0 -0.0482778 -0.553118 0.45214
|
||||
121 10 0 0 0 -0.0554938 -0.402692 0.141112
|
||||
122 10.1 0 0 0 0.174338 0.556958 -0.0922154
|
||||
123 10.2 0 0 0 -1.06045 0.541565 -0.0409312
|
||||
124 10.3 0 0 0 -1.20782 0.464574 -0.413871
|
||||
125 10.4 0 0 0 -0.891701 0.327653 -0.286438
|
||||
126 10.5 0 0 0 0.231227 -0.064277 -0.89684
|
||||
127 10.6 -1.27989 -4.87365 9.40433 0.211278 0.230826 -1.23536
|
||||
128 10.7 -2.1001 -0.417817 1.17745 0.425856 0.078728 -1.44229
|
||||
129 10.8 0 0 0 0.30965 0.450884 -1.74985
|
||||
130 10.9 0 0 0 0.36735 0.990032 -1.19971
|
||||
131 11 0.253834 -1.84303 3.91828 1.01826 0.0660896 -0.481086
|
||||
132 11.1 0 0 0 0.744006 0.0906555 -0.897417
|
||||
133 11.2 0 0 0 0.339073 0.361038 -0.545084
|
||||
134 11.3 -1.9974 -0.431998 3.46296 0.611295 0.17282 0.0341483
|
||||
135 11.4 0 0 0 -0.491432 -0.958871 1.28001
|
||||
136 11.5 0 0 0 0.0431048 -1.50924 1.24037
|
||||
137 11.6 0 0 0 -0.684419 -0.0163951 1.06179
|
||||
138 11.7 0 0 0 -0.425278 -0.127741 0.757298
|
||||
139 11.8 -2.09164 0.00894897 2.22812 -0.0955178 -0.310572 0.661289
|
||||
140 11.9 0 0 0 0.156959 -0.233409 0.802568
|
||||
141 12 0 0 0 -0.05541 -0.346448 0.541571
|
||||
142 12.1 0 0 0 0.706767 0.182767 0.25767
|
||||
143 12.2 0 0 0 0.4791 0.464612 -0.212887
|
||||
144 12.3 0 0 0 0.81454 0.440323 -0.461359
|
||||
145 12.4 0 0 0 -0.110025 0.200698 -0.996706
|
||||
146 12.5 0 0 0 -0.149791 0.165599 -1.02233
|
||||
147 12.6 0 0 0 -0.170933 0.0644682 -0.866174
|
||||
148 12.7 0 0 0 -0.122869 -0.0196287 -0.801348
|
||||
149 12.8 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
150 12.9 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
151 13 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
152 13.1 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
153 13.2 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
154 13.3 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
155 13.4 0 0 0 -0.0693832 -0.0673091 -0.382802
|
||||
156 13.5 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
157 13.6 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
158 13.7 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
159 13.8 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
160 13.9 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
161 14 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
162 14.1 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
163 14.2 0 0 0 -0.000502433 0.000137492 -0.227425
|
||||
164 14.3 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
165 14.4 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
166 14.5 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
167 14.6 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
168 14.7 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
169 14.8 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
170 14.9 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
171 15 0 0 0 5.79042e-05 4.68687e-05 -0.00286094
|
||||
172 15.1 0 0 0 0 0 0
|
||||
173 15.2 0 0 0 0 0 0
|
||||
174 15.3 0 0 0 0 0 0
|
||||
175 15.4 0 0 0 0 0 0
|
||||
176 15.5 0 0 0 0 0 0
|
||||
177 15.6 0 0 0 0 0 0
|
||||
178 15.7 0 0 0 0 0 0
|
||||
179 15.8 0 0 0 0 0 0
|
||||
180 15.9 0 0 0 0 0 0
|
||||
181 16 0 0 0 0 0 0
|
||||
3
examples/USER/misc/mop/mopz0.time.reference
Normal file
3
examples/USER/misc/mop/mopz0.time.reference
Normal file
@ -0,0 +1,3 @@
|
||||
# Time-averaged data for fix mopz0t
|
||||
# TimeStep c_mopz0[1] c_mopz0[2] c_mopz0[3] c_mopz0[4] c_mopz0[5] c_mopz0[6]
|
||||
0 1.62633 1.08234 0.844097 1.18575 -0.408792 -0.752394
|
||||
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -302,6 +302,10 @@
|
||||
/compute_rigid_local.h
|
||||
/compute_spec_atom.cpp
|
||||
/compute_spec_atom.h
|
||||
/compute_stress_mop.cpp
|
||||
/compute_stress_mop.h
|
||||
/compute_stress_mop.profile.cpp
|
||||
/compute_stress_mop.profile.h
|
||||
/compute_stress_tally.cpp
|
||||
/compute_stress_tally.h
|
||||
/compute_temp_asphere.cpp
|
||||
|
||||
@ -30,6 +30,8 @@ compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
|
||||
compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
|
||||
compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017
|
||||
compute entropy/atom, Pablo Piaggi (EPFL), pablo.piaggi at epfl.ch, 15 June 2018
|
||||
compute stress/mop, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
|
||||
compute stress/mop/profile, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
|
||||
compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
||||
compute PRESSURE/GREM, David Stelter, dstelter@bu.edu, 22 Nov 16
|
||||
dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||
|
||||
424
src/USER-MISC/compute_stress_mop.cpp
Normal file
424
src/USER-MISC/compute_stress_mop.cpp
Normal file
@ -0,0 +1,424 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
|
||||
--------------------------------------------------------------------------*/
|
||||
|
||||
#include <mpi.h>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
|
||||
#include "compute_stress_mop.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "domain.h"
|
||||
#include "group.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "neighbor.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
#include "error.h"
|
||||
#include "memory.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{X,Y,Z};
|
||||
enum{TOTAL,CONF,KIN};
|
||||
|
||||
#define BIG 1000000000
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 6) error->all(FLERR,"Illegal compute stress/mop command");
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
// set compute mode and direction of plane(s) for pressure calculation
|
||||
|
||||
if (strcmp(arg[3],"x")==0) {
|
||||
dir = X;
|
||||
} else if (strcmp(arg[3],"y")==0) {
|
||||
dir = Y;
|
||||
} else if (strcmp(arg[3],"z")==0) {
|
||||
dir = Z;
|
||||
} else error->all(FLERR,"Illegal compute stress/mop command");
|
||||
|
||||
// Position of the plane
|
||||
|
||||
if (strcmp(arg[4],"lower")==0) {
|
||||
pos = domain->boxlo[dir];
|
||||
} else if (strcmp(arg[4],"upper")==0) {
|
||||
pos = domain->boxhi[dir];
|
||||
} else if (strcmp(arg[4],"center")==0) {
|
||||
pos = 0.5*(domain->boxlo[dir]+domain->boxhi[dir]);
|
||||
} else pos = force->numeric(FLERR,arg[4]);
|
||||
|
||||
if ( pos < (domain->boxlo[dir]+domain->prd_half[dir]) ) {
|
||||
pos1 = pos + domain->prd[dir];
|
||||
} else {
|
||||
pos1 = pos - domain->prd[dir];
|
||||
}
|
||||
|
||||
// parse values until one isn't recognized
|
||||
|
||||
which = new int[3*(narg-5)];
|
||||
nvalues = 0;
|
||||
int i;
|
||||
|
||||
int iarg=5;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"conf") == 0) {
|
||||
for (i=0; i<3; i++) {
|
||||
which[nvalues] = CONF;
|
||||
nvalues++;
|
||||
}
|
||||
} else if (strcmp(arg[iarg],"kin") == 0) {
|
||||
for (i=0; i<3; i++) {
|
||||
which[nvalues] = KIN;
|
||||
nvalues++;
|
||||
}
|
||||
} else if (strcmp(arg[iarg],"total") == 0) {
|
||||
for (i=0; i<3; i++) {
|
||||
which[nvalues] = TOTAL;
|
||||
nvalues++;
|
||||
}
|
||||
} else error->all(FLERR, "Illegal compute stress/mop command"); //break;
|
||||
|
||||
iarg++;
|
||||
}
|
||||
|
||||
// Error check
|
||||
// 3D only
|
||||
|
||||
if (domain->dimension < 3)
|
||||
error->all(FLERR, "Compute stress/mop incompatible with simulation dimension");
|
||||
|
||||
// orthogonal simulation box
|
||||
if (domain->triclinic != 0)
|
||||
error->all(FLERR, "Compute stress/mop incompatible with triclinic simulation box");
|
||||
// plane inside the box
|
||||
if (pos >domain->boxhi[dir] || pos <domain->boxlo[dir])
|
||||
error->all(FLERR, "Plane for compute stress/mop is out of bounds");
|
||||
|
||||
|
||||
// Initialize some variables
|
||||
|
||||
values_local = values_global = vector = NULL;
|
||||
|
||||
// this fix produces a global vector
|
||||
|
||||
memory->create(vector,nvalues,"stress/mop:vector");
|
||||
memory->create(values_local,nvalues,"stress/mop/spatial:values_local");
|
||||
memory->create(values_global,nvalues,"stress/mop/spatial:values_global");
|
||||
size_vector = nvalues;
|
||||
|
||||
vector_flag = 1;
|
||||
extvector = 0;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeStressMop::~ComputeStressMop()
|
||||
{
|
||||
|
||||
delete [] which;
|
||||
|
||||
memory->destroy(values_local);
|
||||
memory->destroy(values_global);
|
||||
memory->destroy(vector);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMop::init()
|
||||
{
|
||||
|
||||
// Conversion constants
|
||||
|
||||
nktv2p = force->nktv2p;
|
||||
ftm2v = force->ftm2v;
|
||||
|
||||
// Plane area
|
||||
|
||||
area = 1;
|
||||
int i;
|
||||
for (i=0; i<3; i++){
|
||||
if (i!=dir) area = area*domain->prd[i];
|
||||
}
|
||||
|
||||
// Timestep Value
|
||||
|
||||
dt = update->dt;
|
||||
|
||||
// Error check
|
||||
|
||||
// Compute stress/mop requires fixed simulation box
|
||||
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
|
||||
error->all(FLERR, "Compute stress/mop requires a fixed simulation box");
|
||||
|
||||
// This compute requires a pair style with pair_single method implemented
|
||||
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"No pair style is defined for compute stress/mop");
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all(FLERR,"Pair style does not support compute stress/mop");
|
||||
|
||||
// Warnings
|
||||
|
||||
if (me==0){
|
||||
|
||||
//Compute stress/mop only accounts for pair interactions.
|
||||
// issue a warning if any intramolecular potential or Kspace is defined.
|
||||
|
||||
if (force->bond!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop does not account for bond potentials");
|
||||
if (force->angle!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop does not account for angle potentials");
|
||||
if (force->dihedral!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop does not account for dihedral potentials");
|
||||
if (force->improper!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop does not account for improper potentials");
|
||||
if (force->kspace!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop does not account for kspace contributions");
|
||||
}
|
||||
|
||||
// need an occasional half neighbor list
|
||||
int irequest = neighbor->request((void *) this);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMop::init_list(int id, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute output vector
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMop::compute_vector()
|
||||
{
|
||||
invoked_array = update->ntimestep;
|
||||
|
||||
//Compute pressures on separate procs
|
||||
compute_pairs();
|
||||
|
||||
// sum pressure contributions over all procs
|
||||
MPI_Allreduce(values_local,values_global,nvalues,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
int m;
|
||||
for (m=0; m<nvalues; m++) {
|
||||
vector[m] = values_global[m];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
compute pressure contribution of local proc
|
||||
-------------------------------------------------------------------------*/
|
||||
|
||||
void ComputeStressMop::compute_pairs()
|
||||
|
||||
{
|
||||
int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
|
||||
double delx,dely,delz;
|
||||
double rsq,eng,fpair,factor_coul,factor_lj;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
|
||||
// zero out arrays for one sample
|
||||
|
||||
for (i = 0; i < nvalues; i++) values_local[i] = 0.0;
|
||||
|
||||
// invoke half neighbor list (will copy or build if necessary)
|
||||
|
||||
neighbor->build_one(list);
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
Pair *pair = force->pair;
|
||||
double **cutsq = force->pair->cutsq;
|
||||
|
||||
// Parse values
|
||||
|
||||
double xi[3];
|
||||
double vi[3];
|
||||
double fi[3];
|
||||
double xj[3];
|
||||
|
||||
m = 0;
|
||||
while (m<nvalues) {
|
||||
if (which[m] == CONF || which[m] == TOTAL) {
|
||||
|
||||
//Compute configurational contribution to pressure
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
|
||||
xi[0] = atom->x[i][0];
|
||||
xi[1] = atom->x[i][1];
|
||||
xi[2] = atom->x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
// skip if neither I nor J are in group
|
||||
if (!(mask[i] & groupbit || mask[j] & groupbit)) continue;
|
||||
|
||||
xj[0] = atom->x[j][0];
|
||||
xj[1] = atom->x[j][1];
|
||||
xj[2] = atom->x[j][2];
|
||||
delx = xi[0] - xj[0];
|
||||
dely = xi[1] - xj[1];
|
||||
delz = xi[2] - xj[2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
if (rsq >= cutsq[itype][jtype]) continue;
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
|
||||
//check if ij pair is accross plane, add contribution to pressure
|
||||
if ( ((xi[dir]>pos) && (xj[dir]<pos)) || ((xi[dir]>pos1) && (xj[dir]<pos1)) ) {
|
||||
|
||||
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
|
||||
values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p;
|
||||
values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
|
||||
values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
|
||||
}
|
||||
else if ( ((xi[dir]<pos) && (xj[dir]>pos)) || ((xi[dir]<pos1) && (xj[dir]>pos1)) ){
|
||||
|
||||
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
|
||||
values_local[m] -= fpair*(xi[0]-xj[0])/area*nktv2p;
|
||||
values_local[m+1] -= fpair*(xi[1]-xj[1])/area*nktv2p;
|
||||
values_local[m+2] -= fpair*(xi[2]-xj[2])/area*nktv2p;
|
||||
}
|
||||
|
||||
} else {
|
||||
|
||||
if ( ((xi[dir]>pos) && (xj[dir]<pos)) || ((xi[dir]>pos1) && (xj[dir]<pos1)) ) {
|
||||
|
||||
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
|
||||
values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p;
|
||||
values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
|
||||
values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
// Compute kinetic contribution to pressure
|
||||
// counts local particles transfers across the plane
|
||||
|
||||
if (which[m] == KIN || which[m] == TOTAL){
|
||||
double vcm[3];
|
||||
double masstotal,sgn;
|
||||
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
|
||||
// skip if I is not in group
|
||||
if (mask[i] & groupbit){
|
||||
|
||||
itype = type[i];
|
||||
|
||||
//coordinates at t
|
||||
xi[0] = atom->x[i][0];
|
||||
xi[1] = atom->x[i][1];
|
||||
xi[2] = atom->x[i][2];
|
||||
|
||||
//velocities at t
|
||||
vi[0] = atom->v[i][0];
|
||||
vi[1] = atom->v[i][1];
|
||||
vi[2] = atom->v[i][2];
|
||||
|
||||
//forces at t
|
||||
fi[0] = atom->f[i][0];
|
||||
fi[1] = atom->f[i][1];
|
||||
fi[2] = atom->f[i][2];
|
||||
|
||||
//coordinates at t-dt (based on Velocity-Verlet alg.)
|
||||
xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;
|
||||
xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;
|
||||
xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
|
||||
|
||||
// because LAMMPS does not put atoms back in the box
|
||||
// at each timestep, must check atoms going through the
|
||||
// image of the plane that is closest to the box
|
||||
|
||||
double pos_temp = pos+copysign(1,domain->prd_half[dir]-pos)*domain->prd[dir];
|
||||
if (fabs(xi[dir]-pos)<fabs(xi[dir]-pos_temp)) pos_temp = pos;
|
||||
|
||||
if (((xi[dir]-pos_temp)*(xj[dir]-pos_temp)<0)){
|
||||
|
||||
//sgn = copysign(1,vi[dir]-vcm[dir]);
|
||||
sgn = copysign(1,vi[dir]);
|
||||
|
||||
//approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)
|
||||
double vcross[3];
|
||||
vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
|
||||
vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
|
||||
vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;
|
||||
|
||||
values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
|
||||
values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
|
||||
values_local[m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
m+=3;
|
||||
}
|
||||
}
|
||||
|
||||
102
src/USER-MISC/compute_stress_mop.h
Normal file
102
src/USER-MISC/compute_stress_mop.h
Normal file
@ -0,0 +1,102 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
|
||||
--------------------------------------------------------------------------*/
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(stress/mop,ComputeStressMop)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_STRESS_MOP_H
|
||||
#define LMP_COMPUTE_STRESS_MOP_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeStressMop : public Compute {
|
||||
public:
|
||||
ComputeStressMop(class LAMMPS *, int, char **);
|
||||
virtual ~ComputeStressMop();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void compute_vector();
|
||||
|
||||
private:
|
||||
|
||||
void compute_pairs();
|
||||
|
||||
int me,nvalues,dir;
|
||||
int *which;
|
||||
|
||||
double *values_local,*values_global;
|
||||
double pos,pos1,dt,nktv2p,ftm2v;
|
||||
double area;
|
||||
class NeighList *list;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Compute stress/mop incompatible with simulation dimension
|
||||
|
||||
Compute stress/mop only works with 3D simulations.
|
||||
|
||||
E: Compute stress/mop incompatible with triclinic simulation box
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute stress/mop requires a fixed simulation box
|
||||
|
||||
Compute stress/mop is not compatible with any change of volume or shape
|
||||
or boundary conditions of the simulation box.
|
||||
|
||||
E: No pair style is defined for compute stress/mop
|
||||
|
||||
Self-explanatory. Compute stress/mop requires the definition of a pair style.
|
||||
|
||||
E: Pair style does not support compute stress/mop
|
||||
|
||||
The pair style does not have a single() function, so it can
|
||||
not be invoked by compute stress/mop.
|
||||
|
||||
W: compute stress/mop does not account for bond potentials
|
||||
|
||||
W: compute stress/mop does not account for angle potentials
|
||||
|
||||
W: compute stress/mop does not account for dihedral potentials
|
||||
|
||||
W: compute stress/mop does not account for improper potentials
|
||||
|
||||
W: compute stress/mop does not account for kspace contributions
|
||||
|
||||
Compute stress/mop only accounts for pairwise additive interactions for
|
||||
the computation of local stress tensor components.
|
||||
|
||||
*/
|
||||
|
||||
496
src/USER-MISC/compute_stress_mop_profile.cpp
Normal file
496
src/USER-MISC/compute_stress_mop_profile.cpp
Normal file
@ -0,0 +1,496 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
|
||||
--------------------------------------------------------------------------*/
|
||||
|
||||
#include <mpi.h>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
|
||||
#include "compute_stress_mop_profile.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "domain.h"
|
||||
#include "group.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "neighbor.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
#include "error.h"
|
||||
#include "memory.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{X,Y,Z};
|
||||
enum{LOWER,CENTER,UPPER,COORD};
|
||||
enum{TOTAL,CONF,KIN};
|
||||
|
||||
#define BIG 1000000000
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 7) error->all(FLERR,"Illegal compute stress/mop/profile command");
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
// set compute mode and direction of plane(s) for pressure calculation
|
||||
|
||||
if (strcmp(arg[3],"x")==0) {
|
||||
dir = X;
|
||||
} else if (strcmp(arg[3],"y")==0) {
|
||||
dir = Y;
|
||||
} else if (strcmp(arg[3],"z")==0) {
|
||||
dir = Z;
|
||||
} else error->all(FLERR,"Illegal compute stress/mop/profile command");
|
||||
|
||||
// bin parameters
|
||||
|
||||
if (strcmp(arg[4],"lower") == 0) originflag = LOWER;
|
||||
else if (strcmp(arg[4],"center") == 0) originflag = CENTER;
|
||||
else if (strcmp(arg[4],"upper") == 0) originflag = UPPER;
|
||||
else originflag = COORD;
|
||||
if (originflag == COORD)
|
||||
origin = force->numeric(FLERR,arg[4]);
|
||||
delta = force->numeric(FLERR,arg[5]);
|
||||
invdelta = 1.0/delta;
|
||||
|
||||
// parse values until one isn't recognized
|
||||
|
||||
which = new int[3*(narg-6)];
|
||||
nvalues = 0;
|
||||
int i;
|
||||
|
||||
int iarg=6;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"conf") == 0) {
|
||||
for (i=0; i<3; i++) {
|
||||
which[nvalues] = CONF;
|
||||
nvalues++;
|
||||
}
|
||||
} else if (strcmp(arg[iarg],"kin") == 0) {
|
||||
for (i=0; i<3; i++) {
|
||||
which[nvalues] = KIN;
|
||||
nvalues++;
|
||||
}
|
||||
} else if (strcmp(arg[iarg],"total") == 0) {
|
||||
for (i=0; i<3; i++) {
|
||||
which[nvalues] = TOTAL;
|
||||
nvalues++;
|
||||
}
|
||||
} else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break;
|
||||
|
||||
iarg++;
|
||||
}
|
||||
|
||||
// check domain related errors
|
||||
|
||||
// 3D only
|
||||
|
||||
if (domain->dimension < 3)
|
||||
error->all(FLERR, "Compute stress/mop/profile incompatible with simulation dimension");
|
||||
|
||||
// orthogonal simulation box
|
||||
|
||||
if (domain->triclinic != 0)
|
||||
error->all(FLERR, "Compute stress/mop/profile incompatible with triclinic simulation box");
|
||||
|
||||
// initialize some variables
|
||||
|
||||
nbins = 0;
|
||||
coord = coordp = NULL;
|
||||
values_local = values_global = array = NULL;
|
||||
|
||||
// bin setup
|
||||
|
||||
setup_bins();
|
||||
|
||||
// this fix produces a global array
|
||||
|
||||
memory->create(array,nbins,1+nvalues,"stress/mop/profile:array");
|
||||
size_array_rows = nbins;
|
||||
size_array_cols = 1 + nvalues;
|
||||
|
||||
array_flag = 1;
|
||||
extarray = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeStressMopProfile::~ComputeStressMopProfile()
|
||||
{
|
||||
|
||||
delete [] which;
|
||||
|
||||
memory->destroy(coord);
|
||||
memory->destroy(coordp);
|
||||
memory->destroy(values_local);
|
||||
memory->destroy(values_global);
|
||||
memory->destroy(array);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMopProfile::init()
|
||||
{
|
||||
|
||||
// conversion constants
|
||||
|
||||
nktv2p = force->nktv2p;
|
||||
ftm2v = force->ftm2v;
|
||||
|
||||
// plane area
|
||||
|
||||
area = 1;
|
||||
int i;
|
||||
for (i=0; i<3; i++) {
|
||||
if (i!=dir) area = area*domain->prd[i];
|
||||
}
|
||||
|
||||
// timestep Value
|
||||
|
||||
dt = update->dt;
|
||||
|
||||
// Error check
|
||||
// Compute stress/mop/profile requires fixed simulation box
|
||||
|
||||
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
|
||||
error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box");
|
||||
|
||||
//This compute requires a pair style with pair_single method implemented
|
||||
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"No pair style is defined for compute stress/mop/profile");
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all(FLERR,"Pair style does not support compute stress/mop/profile");
|
||||
|
||||
// Warnings
|
||||
|
||||
if (me==0){
|
||||
|
||||
//Compute stress/mop/profile only accounts for pair interactions.
|
||||
// issue a warning if any intramolecular potential or Kspace is defined.
|
||||
|
||||
if (force->bond!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop/profile does not account for bond potentials");
|
||||
if (force->angle!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop/profile does not account for angle potentials");
|
||||
if (force->dihedral!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop/profile does not account for dihedral potentials");
|
||||
if (force->improper!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop/profile does not account for improper potentials");
|
||||
if (force->kspace!=NULL)
|
||||
error->warning(FLERR,"compute stress/mop/profile does not account for kspace contributions");
|
||||
}
|
||||
|
||||
// need an occasional half neighbor list
|
||||
|
||||
int irequest = neighbor->request((void *) this);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMopProfile::init_list(int id, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute output array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMopProfile::compute_array()
|
||||
{
|
||||
invoked_array = update->ntimestep;
|
||||
|
||||
//Compute pressures on separate procs
|
||||
compute_pairs();
|
||||
|
||||
// sum pressure contributions over all procs
|
||||
MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
int ibin,m,mo;
|
||||
for (ibin=0; ibin<nbins; ibin++) {
|
||||
array[ibin][0] = coord[ibin][0];
|
||||
mo=1;
|
||||
|
||||
m = 0;
|
||||
while (m<nvalues) {
|
||||
array[ibin][m+mo] = values_global[ibin][m];
|
||||
m++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
compute pressure contribution of local proc
|
||||
-------------------------------------------------------------------------*/
|
||||
|
||||
void ComputeStressMopProfile::compute_pairs()
|
||||
|
||||
{
|
||||
int i,j,m,n,ii,jj,inum,jnum,itype,jtype,ibin;
|
||||
double delx,dely,delz;
|
||||
double rsq,eng,fpair,factor_coul,factor_lj;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double pos,pos1,pos_temp;
|
||||
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
|
||||
// zero out arrays for one sample
|
||||
for (m = 0; m < nbins; m++) {
|
||||
for (i = 0; i < nvalues; i++) values_local[m][i] = 0.0;
|
||||
}
|
||||
|
||||
// invoke half neighbor list (will copy or build if necessary)
|
||||
|
||||
neighbor->build_one(list);
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
Pair *pair = force->pair;
|
||||
double **cutsq = force->pair->cutsq;
|
||||
|
||||
// parse values
|
||||
|
||||
double xi[3];
|
||||
double vi[3];
|
||||
double fi[3];
|
||||
double xj[3];
|
||||
|
||||
m = 0;
|
||||
while (m<nvalues) {
|
||||
if (which[m] == CONF || which[m] == TOTAL) {
|
||||
|
||||
// Compute configurational contribution to pressure
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
|
||||
xi[0] = atom->x[i][0];
|
||||
xi[1] = atom->x[i][1];
|
||||
xi[2] = atom->x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
// skip if neither I nor J are in group
|
||||
|
||||
if (!(mask[i] & groupbit || mask[j] & groupbit)) continue;
|
||||
|
||||
xj[0] = atom->x[j][0];
|
||||
xj[1] = atom->x[j][1];
|
||||
xj[2] = atom->x[j][2];
|
||||
delx = xi[0] - xj[0];
|
||||
dely = xi[1] - xj[1];
|
||||
delz = xi[2] - xj[2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
if (rsq >= cutsq[itype][jtype]) continue;
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
|
||||
for (ibin=0;ibin<nbins;ibin++) {
|
||||
pos = coord[ibin][0];
|
||||
pos1 = coordp[ibin][0];
|
||||
|
||||
//check if ij pair is accross plane, add contribution to pressure
|
||||
|
||||
if ( ((xi[dir]>pos) && (xj[dir]<pos))
|
||||
|| ((xi[dir]>pos1) && (xj[dir]<pos1)) ) {
|
||||
|
||||
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
|
||||
values_local[ibin][m] += fpair*(xi[0]-xj[0])/area*nktv2p;
|
||||
values_local[ibin][m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
|
||||
values_local[ibin][m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
|
||||
|
||||
} else if ( ((xi[dir]<pos) && (xj[dir]>pos))
|
||||
|| ((xi[dir]<pos1) && (xj[dir]>pos1)) ) {
|
||||
|
||||
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
|
||||
values_local[ibin][m] -= fpair*(xi[0]-xj[0])/area*nktv2p;
|
||||
values_local[ibin][m+1] -= fpair*(xi[1]-xj[1])/area*nktv2p;
|
||||
values_local[ibin][m+2] -= fpair*(xi[2]-xj[2])/area*nktv2p;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
|
||||
for (ibin=0;ibin<nbins;ibin++) {
|
||||
pos = coord[ibin][0];
|
||||
pos1 = coordp[ibin][0];
|
||||
|
||||
//check if ij pair is accross plane, add contribution to pressure
|
||||
|
||||
if ( ((xi[dir]>pos) && (xj[dir]<pos))
|
||||
|| ((xi[dir]>pos1) && (xj[dir]<pos1)) ) {
|
||||
|
||||
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
|
||||
values_local[ibin][m] += fpair*(xi[0]-xj[0])/area*nktv2p;
|
||||
values_local[ibin][m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
|
||||
values_local[ibin][m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// compute kinetic contribution to pressure
|
||||
// counts local particles transfers across the plane
|
||||
|
||||
if (which[m] == KIN || which[m] == TOTAL){
|
||||
|
||||
double vcm[3];
|
||||
double masstotal,sgn;
|
||||
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
|
||||
// skip if I is not in group
|
||||
|
||||
if (mask[i] & groupbit){
|
||||
|
||||
itype = type[i];
|
||||
|
||||
//coordinates at t
|
||||
xi[0] = atom->x[i][0];
|
||||
xi[1] = atom->x[i][1];
|
||||
xi[2] = atom->x[i][2];
|
||||
|
||||
//velocities at t
|
||||
vi[0] = atom->v[i][0];
|
||||
vi[1] = atom->v[i][1];
|
||||
vi[2] = atom->v[i][2];
|
||||
|
||||
//forces at t
|
||||
fi[0] = atom->f[i][0];
|
||||
fi[1] = atom->f[i][1];
|
||||
fi[2] = atom->f[i][2];
|
||||
|
||||
//coordinates at t-dt (based on Velocity-Verlet alg.)
|
||||
xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;
|
||||
xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;
|
||||
xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
|
||||
|
||||
for (ibin=0;ibin<nbins;ibin++) {
|
||||
pos = coord[ibin][0];
|
||||
pos1 = coordp[ibin][0];
|
||||
|
||||
if (((xi[dir]-pos)*(xj[dir]-pos)*(xi[dir]-pos1)*(xj[dir]-pos1)<0)){
|
||||
|
||||
sgn = copysign(1,vi[dir]);
|
||||
|
||||
//approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)
|
||||
double vcross[3];
|
||||
vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
|
||||
vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
|
||||
vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;
|
||||
|
||||
values_local[ibin][m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
|
||||
values_local[ibin][m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
|
||||
values_local[ibin][m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
m+=3;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
setup 1d bins and their extent and coordinates
|
||||
called at init()
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeStressMopProfile::setup_bins()
|
||||
{
|
||||
int i,j,k,m,n;
|
||||
double lo,hi,coord1,coord2;
|
||||
|
||||
double *boxlo,*boxhi,*prd;
|
||||
boxlo = domain->boxlo;
|
||||
boxhi = domain->boxhi;
|
||||
prd = domain->prd;
|
||||
|
||||
if (originflag == LOWER) origin = boxlo[dir];
|
||||
else if (originflag == UPPER) origin = boxhi[dir];
|
||||
else if (originflag == CENTER)
|
||||
origin = 0.5 * (boxlo[dir] + boxhi[dir]);
|
||||
|
||||
if (origin < boxlo[dir]) {
|
||||
error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" );
|
||||
} else {
|
||||
n = static_cast<int> ((origin - boxlo[dir]) * invdelta);
|
||||
lo = origin - n*delta;
|
||||
}
|
||||
if (origin < boxhi[dir]) {
|
||||
n = static_cast<int> ((boxhi[dir] - origin) * invdelta);
|
||||
hi = origin + n*delta;
|
||||
} else {
|
||||
error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" );
|
||||
}
|
||||
|
||||
offset = lo;
|
||||
nbins = static_cast<int> ((hi-lo) * invdelta + 1.5);
|
||||
|
||||
//allocate bin arrays
|
||||
memory->create(coord,nbins,1,"stress/mop/profile:coord");
|
||||
memory->create(coordp,nbins,1,"stress/mop/profile:coordp");
|
||||
memory->create(values_local,nbins,nvalues,"stress/mop/profile:values_local");
|
||||
memory->create(values_global,nbins,nvalues,"stress/mop/profile:values_global");
|
||||
|
||||
// set bin coordinates
|
||||
for (i = 0; i < nbins; i++) {
|
||||
coord[i][0] = offset + i*delta;
|
||||
if ( coord[i][0] < (domain->boxlo[dir]+domain->prd_half[dir]) ) {
|
||||
coordp[i][0] = coord[i][0] + domain->prd[dir];
|
||||
} else {
|
||||
coordp[i][0] = coord[i][0] - domain->prd[dir];
|
||||
}
|
||||
}
|
||||
}
|
||||
113
src/USER-MISC/compute_stress_mop_profile.h
Normal file
113
src/USER-MISC/compute_stress_mop_profile.h
Normal file
@ -0,0 +1,113 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
|
||||
--------------------------------------------------------------------------*/
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(stress/mop/profile,ComputeStressMopProfile)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_STRESS_MOP_PROFILE_H
|
||||
#define LMP_COMPUTE_STRESS_MOP_PROFILE_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeStressMopProfile : public Compute {
|
||||
public:
|
||||
ComputeStressMopProfile(class LAMMPS *, int, char **);
|
||||
virtual ~ComputeStressMopProfile();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void compute_array();
|
||||
|
||||
private:
|
||||
|
||||
void compute_pairs();
|
||||
void setup_bins();
|
||||
|
||||
int me,nvalues,dir;
|
||||
int *which;
|
||||
|
||||
int originflag;
|
||||
double origin,delta,offset,invdelta;
|
||||
int nbins;
|
||||
double **coord,**coordp;
|
||||
double **values_local,**values_global;
|
||||
|
||||
int ndim;
|
||||
double dt,nktv2p,ftm2v;
|
||||
double area;
|
||||
class NeighList *list;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Compute stress/mop/profile incompatible with simulation dimension
|
||||
|
||||
Compute stress/mop/profile only works with 3D simulations.
|
||||
|
||||
E: Compute stress/mop/profile incompatible with triclinic simulation box
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute stress/mop/profile requires a fixed simulation box
|
||||
|
||||
Compute stress/mop/profile is not compatible with any change of volume or shape
|
||||
or boundary conditions of the simulation box.
|
||||
|
||||
E: No pair style is defined for compute stress/mop/profile
|
||||
|
||||
Self-explanatory. Compute stress/mop/profile requires the definition of a pair style.
|
||||
|
||||
E: Pair style does not support compute stress/mop/profile
|
||||
|
||||
The pair style does not have a single() function, so it can
|
||||
not be invoked by compute stress/mop/profile.
|
||||
|
||||
E: Origin of bins for compute stress/mop/profile is out of bounds
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
W: compute stress/mop/profile does not account for bond potentials
|
||||
|
||||
W: compute stress/mop/profile does not account for angle potentials
|
||||
|
||||
W: compute stress/mop/profile does not account for dihedral potentials
|
||||
|
||||
W: compute stress/mop/profile does not account for improper potentials
|
||||
|
||||
W: compute stress/mop/profile does not account for kspace contributions
|
||||
|
||||
Compute stress/mop/profile only accounts for pairwise additive interactions for
|
||||
the computation of local stress tensor components.
|
||||
|
||||
*/
|
||||
|
||||
Reference in New Issue
Block a user