mirror of https://github.com/lammps/lammps.git
update docs and references for name changes in USER-MOP package, remove obsoleted files
This commit is contained in:
Axel Kohlmeyer
parent
a797a0d193
commit
cb4ffaf95c
|
|
@ -1577,11 +1577,11 @@ USER-MOP package :link(PKG-USER-MOP),h4
|
|||
|
||||
[Contents:]
|
||||
|
||||
This package contains two compute styles: compute mop and compute mop/profile.
|
||||
This package contains two compute styles: compute stress/mop and compute stress/mop/profile.
|
||||
|
||||
compute mop calculates components of the local stress tensor using the method of planes, applied on one plane.
|
||||
compute stress/mop calculates components of the local stress tensor using the method of planes, applied on one plane.
|
||||
|
||||
compute mop/profile calculates profiles of components of the local stress tensor using the method of planes, applied on an array of regularly spaced planes.
|
||||
compute stress/mop/profile calculates profiles of components of the local stress tensor using the method of planes, applied on an array of regularly spaced planes.
|
||||
|
||||
[Author:] Laurent Joly (University of Lyon 1)
|
||||
Romain Vermorel at (University of Pau and Pays de l'Adour)
|
||||
|
|
@ -1590,8 +1590,8 @@ Romain Vermorel at (University of Pau and Pays de l'Adour)
|
|||
|
||||
src/USER-MOP: filenames -> commands
|
||||
src/USER-MOP/README
|
||||
"compute mop"_compute_mop.html
|
||||
"compute mop/profile"_compute_mop.html
|
||||
"compute stress/mop"_compute_stress_mop.html
|
||||
"compute stress/mop/profile"_compute_stress_mop.html
|
||||
examples/USER/mop :ul
|
||||
|
||||
:line
|
||||
|
|
|
|||
|
|
@ -57,7 +57,7 @@ Package, Description, Doc page, Example, Library
|
|||
"USER-MESO"_Packages_details.html#PKG-USER-MESO, mesoscale DPD models, "pair_style edpd"_pair_meso.html, USER/meso, no
|
||||
"USER-MGPT"_Packages_details.html#PKG-USER-MGPT, fast MGPT multi-ion potentials, "pair_style mgpt"_pair_mgpt.html, USER/mgpt, no
|
||||
"USER-MISC"_Packages_details.html#PKG-USER-MISC, single-file contributions, USER-MISC/README, USER/misc, no
|
||||
"USER-MOP"_Packages_details.html#PKG-USER-MOP, compute stress via method-of-planes, "compute_style mop"_compute_mop.html, USER/mop, no
|
||||
"USER-MOP"_Packages_details.html#PKG-USER-MOP, compute stress via method-of-planes, "compute_style stress/mop"_compute_stress_mop.html, USER/mop, no
|
||||
"USER-MOFFF"_Packages_details.html#PKG-USER-MOFFF, styles for "MOF-FF"_MOFplus force field, "pair_style buck6d/coul/gauss"_pair_buck6d_coul_gauss.html, USER/mofff, no
|
||||
"USER-MOLFILE"_Packages_details.html#PKG-USER-MOLFILE, "VMD"_https://www.ks.uiuc.edu/Research/vmd/ molfile plug-ins,"dump molfile"_dump_molfile.html, n/a, ext
|
||||
"USER-NETCDF"_Packages_details.html#PKG-USER-NETCDF, dump output via NetCDF,"dump netcdf"_dump_netcdf.html, n/a, ext
|
||||
|
|
|
|||
|
|
@ -1,111 +0,0 @@
|
|||
"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 mop command :h3
|
||||
compute 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 = {mop} or {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
|
||||
{mop} args = pos
|
||||
pos = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the plane
|
||||
{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 mop x lower total
|
||||
compute 1 liquid 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 mop/profile x lower 0.1 total
|
||||
compute 1 liquid 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 {mop} and compute {mop/profile} define computations that
|
||||
calculate components of the local stress tensor using the method of
|
||||
planes "(Todd)"_#mop-todd. Specifically in compute {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 {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 {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 {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:]
|
||||
|
||||
This style is part of the USER-MOP package. It is 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).
|
||||
|
|
@ -55,7 +55,6 @@ Computes :h1
|
|||
compute_meso_e_atom
|
||||
compute_meso_rho_atom
|
||||
compute_meso_t_atom
|
||||
compute_mop
|
||||
compute_msd
|
||||
compute_msd_chunk
|
||||
compute_msd_nongauss
|
||||
|
|
@ -99,6 +98,7 @@ Computes :h1
|
|||
compute_sna_atom
|
||||
compute_spin
|
||||
compute_stress_atom
|
||||
compute_stress_mop
|
||||
compute_tally
|
||||
compute_tdpd_cc_atom
|
||||
compute_temp
|
||||
|
|
|
|||
|
|
@ -446,7 +446,6 @@ compute_ke_rigid.html
|
|||
compute_meso_e_atom.html
|
||||
compute_meso_rho_atom.html
|
||||
compute_meso_t_atom.html
|
||||
compute_mop.html
|
||||
compute_msd.html
|
||||
compute_msd_chunk.html
|
||||
compute_msd_nongauss.html
|
||||
|
|
@ -490,6 +489,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
|
||||
|
|
|
|||
|
|
@ -1,62 +0,0 @@
|
|||
variable T equal 0.8
|
||||
variable p_solid equal 0.05
|
||||
|
||||
lattice fcc 1.0
|
||||
region box block 0.0 6.0 0.0 6.0 -2.0 12.0
|
||||
create_box 2 box
|
||||
|
||||
mass * 1.0
|
||||
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
|
||||
|
||||
region solid_bottom block INF INF INF INF -1.1 0.1
|
||||
region liquid block INF INF INF INF 1.1 8.9
|
||||
region solid_up block INF INF INF INF 9.9 11.1
|
||||
|
||||
create_atoms 1 region liquid
|
||||
delete_atoms porosity liquid 0.26 88765
|
||||
group liquid region liquid
|
||||
|
||||
create_atoms 2 region solid_bottom
|
||||
group solid_bottom region solid_bottom
|
||||
create_atoms 2 region solid_up
|
||||
group solid_up region solid_up
|
||||
group solid union solid_bottom solid_up
|
||||
|
||||
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 10000
|
||||
thermo_modify flush yes temp Tliq
|
||||
|
||||
# dump 1 all atom 10000 dump.lammpstrj
|
||||
|
||||
fix fxbal all balance 1000 1.05 shift z 10 1.05
|
||||
velocity liquid create $T 47298 dist gaussian rot yes
|
||||
run 50000
|
||||
# undump 1
|
||||
reset_timestep 0
|
||||
|
||||
compute bin_z liquid chunk/atom bin/1d z 0.0 0.1 units box
|
||||
|
||||
compute liquidStress_ke liquid stress/atom NULL ke
|
||||
compute liquidStress_vir liquid stress/atom NULL virial
|
||||
fix profile_z liquid ave/chunk 10 1000 10000 bin_z density/number temp c_liquidStress_ke[1] c_liquidStress_ke[2] c_liquidStress_ke[3] c_liquidStress_ke[4] c_liquidStress_ke[5] c_liquidStress_ke[6] c_liquidStress_vir[1] c_liquidStress_vir[2] c_liquidStress_vir[3] c_liquidStress_vir[4] c_liquidStress_vir[5] c_liquidStress_vir[6] ave running overwrite file profile.z
|
||||
|
||||
compute mopz0 all mop z center kin conf
|
||||
fix mopz0t all ave/time 10 1000 10000 c_mopz0[*] file mopz0.time
|
||||
|
||||
compute moppz liquid mop/profile z 0.0 0.1 kin conf
|
||||
fix moppzt all ave/time 100 100 10000 c_moppz[*] ave running overwrite file moppz.time mode vector
|
||||
|
||||
run 40000
|
||||
|
||||
|
|
@ -1,186 +0,0 @@
|
|||
LAMMPS (31 Aug 2018)
|
||||
variable T equal 0.8
|
||||
variable p_solid equal 0.05
|
||||
|
||||
lattice fcc 1.0
|
||||
Lattice spacing in x,y,z = 1.5874 1.5874 1.5874
|
||||
region box block 0.0 6.0 0.0 6.0 -2.0 12.0
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 -3.1748) to (9.52441 9.52441 19.0488)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
|
||||
mass * 1.0
|
||||
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
|
||||
|
||||
region solid_bottom block INF INF INF INF -1.1 0.1
|
||||
region liquid block INF INF INF INF 1.1 8.9
|
||||
region solid_up block INF INF INF INF 9.9 11.1
|
||||
|
||||
create_atoms 1 region liquid
|
||||
Created 1080 atoms
|
||||
Time spent = 0.000631094 secs
|
||||
delete_atoms porosity liquid 0.26 88765
|
||||
Deleted 257 atoms, new total = 823
|
||||
group liquid region liquid
|
||||
823 atoms in group liquid
|
||||
|
||||
create_atoms 2 region solid_bottom
|
||||
Created 216 atoms
|
||||
Time spent = 0.000221014 secs
|
||||
group solid_bottom region solid_bottom
|
||||
216 atoms in group solid_bottom
|
||||
create_atoms 2 region solid_up
|
||||
Created 216 atoms
|
||||
Time spent = 0.000169039 secs
|
||||
group solid_up region solid_up
|
||||
216 atoms in group solid_up
|
||||
group solid union solid_bottom solid_up
|
||||
432 atoms in group solid
|
||||
|
||||
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 (../fix_nh.cpp:1404)
|
||||
|
||||
thermo 10000
|
||||
thermo_modify flush yes temp Tliq
|
||||
WARNING: Temperature for thermo pressure is not for group all (../thermo.cpp:488)
|
||||
|
||||
# dump 1 all atom 10000 dump.lammpstrj
|
||||
|
||||
fix fxbal all balance 1000 1.05 shift z 10 1.05
|
||||
velocity liquid create $T 47298 dist gaussian rot yes
|
||||
velocity liquid create 0.8 47298 dist gaussian rot yes
|
||||
run 50000
|
||||
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 16
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.187 | 3.187 | 3.187 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 0.8 -3.7188573 0 -2.9328812 -0.65673631 2016
|
||||
10000 0.80645968 -3.2160592 0 -2.4237366 -0.17266146 2016
|
||||
20000 0.85787011 -3.2087113 0 -2.3658796 0.032152534 2016
|
||||
30000 0.76357074 -3.1761807 0 -2.4259953 -0.034009443 2016
|
||||
40000 0.8116632 -3.1965608 0 -2.399126 -0.030033303 2016
|
||||
50000 0.83600463 -3.1878245 0 -2.3664749 0.054158845 2016
|
||||
Loop time of 40.6472 on 1 procs for 50000 steps with 1255 atoms
|
||||
|
||||
Performance: 531401.825 tau/day, 1230.097 timesteps/s
|
||||
99.7% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 28.102 | 28.102 | 28.102 | 0.0 | 69.14
|
||||
Neigh | 9.8377 | 9.8377 | 9.8377 | 0.0 | 24.20
|
||||
Comm | 0.91642 | 0.91642 | 0.91642 | 0.0 | 2.25
|
||||
Output | 0.00024581 | 0.00024581 | 0.00024581 | 0.0 | 0.00
|
||||
Modify | 1.4027 | 1.4027 | 1.4027 | 0.0 | 3.45
|
||||
Other | | 0.3879 | | | 0.95
|
||||
|
||||
Nlocal: 1255 ave 1255 max 1255 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 2268 ave 2268 max 2268 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 40149 ave 40149 max 40149 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 40149
|
||||
Ave neighs/atom = 31.9912
|
||||
Neighbor list builds = 5714
|
||||
Dangerous builds = 0
|
||||
# undump 1
|
||||
reset_timestep 0
|
||||
|
||||
compute bin_z liquid chunk/atom bin/1d z 0.0 0.1 units box
|
||||
|
||||
compute liquidStress_ke liquid stress/atom NULL ke
|
||||
compute liquidStress_vir liquid stress/atom NULL virial
|
||||
fix profile_z liquid ave/chunk 10 1000 10000 bin_z density/number temp c_liquidStress_ke[1] c_liquidStress_ke[2] c_liquidStress_ke[3] c_liquidStress_ke[4] c_liquidStress_ke[5] c_liquidStress_ke[6] c_liquidStress_vir[1] c_liquidStress_vir[2] c_liquidStress_vir[3] c_liquidStress_vir[4] c_liquidStress_vir[5] c_liquidStress_vir[6] ave running overwrite file profile.z
|
||||
|
||||
compute mopz0 all mop z center kin conf
|
||||
fix mopz0t all ave/time 10 1000 10000 c_mopz0[*] file mopz0.time
|
||||
|
||||
compute moppz liquid mop/profile z 0.0 0.1 kin conf
|
||||
fix moppzt all ave/time 100 100 10000 c_moppz[*] ave running overwrite file moppz.time mode vector
|
||||
|
||||
run 40000
|
||||
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 16
|
||||
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 mop, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(3) compute 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) = 4.186 | 4.186 | 4.186 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 0.83600463 -3.1878245 0 -2.3664749 0.054158845 2016
|
||||
10000 0.80600425 -3.1900933 0 -2.3982182 -0.042750667 2016
|
||||
20000 0.78117978 -3.2121798 0 -2.444694 0.0092285451 2016
|
||||
30000 0.81008494 -3.1728363 0 -2.376952 0.053347807 2016
|
||||
40000 0.80440016 -3.1768915 0 -2.3865924 -0.0075214347 2016
|
||||
Loop time of 45.4499 on 1 procs for 40000 steps with 1255 atoms
|
||||
|
||||
Performance: 380198.462 tau/day, 880.089 timesteps/s
|
||||
99.8% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 25.456 | 25.456 | 25.456 | 0.0 | 56.01
|
||||
Neigh | 7.9422 | 7.9422 | 7.9422 | 0.0 | 17.47
|
||||
Comm | 0.73477 | 0.73477 | 0.73477 | 0.0 | 1.62
|
||||
Output | 0.00010562 | 0.00010562 | 0.00010562 | 0.0 | 0.00
|
||||
Modify | 11.004 | 11.004 | 11.004 | 0.0 | 24.21
|
||||
Other | | 0.3131 | | | 0.69
|
||||
|
||||
Nlocal: 1255 ave 1255 max 1255 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 2269 ave 2269 max 2269 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 40325 ave 40325 max 40325 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 40325
|
||||
Ave neighs/atom = 32.1315
|
||||
Neighbor list builds = 4587
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:01:26
|
||||
|
|
@ -1,188 +0,0 @@
|
|||
LAMMPS (31 Aug 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
|
||||
|
||||
lattice fcc 1.0
|
||||
Lattice spacing in x,y,z = 1.5874 1.5874 1.5874
|
||||
region box block 0.0 6.0 0.0 6.0 -2.0 12.0
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 -3.1748) to (9.52441 9.52441 19.0488)
|
||||
1 by 1 by 4 MPI processor grid
|
||||
|
||||
mass * 1.0
|
||||
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
|
||||
|
||||
region solid_bottom block INF INF INF INF -1.1 0.1
|
||||
region liquid block INF INF INF INF 1.1 8.9
|
||||
region solid_up block INF INF INF INF 9.9 11.1
|
||||
|
||||
create_atoms 1 region liquid
|
||||
Created 1080 atoms
|
||||
Time spent = 0.000378132 secs
|
||||
delete_atoms porosity liquid 0.26 88765
|
||||
Deleted 288 atoms, new total = 792
|
||||
group liquid region liquid
|
||||
792 atoms in group liquid
|
||||
|
||||
create_atoms 2 region solid_bottom
|
||||
Created 216 atoms
|
||||
Time spent = 0.000345945 secs
|
||||
group solid_bottom region solid_bottom
|
||||
216 atoms in group solid_bottom
|
||||
create_atoms 2 region solid_up
|
||||
Created 216 atoms
|
||||
Time spent = 0.000124454 secs
|
||||
group solid_up region solid_up
|
||||
216 atoms in group solid_up
|
||||
group solid union solid_bottom solid_up
|
||||
432 atoms in group solid
|
||||
|
||||
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 10000
|
||||
thermo_modify flush yes temp Tliq
|
||||
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:488)
|
||||
|
||||
# dump 1 all atom 10000 dump.lammpstrj
|
||||
|
||||
fix fxbal all balance 1000 1.05 shift z 10 1.05
|
||||
velocity liquid create $T 47298 dist gaussian rot yes
|
||||
velocity liquid create 0.8 47298 dist gaussian rot yes
|
||||
run 50000
|
||||
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 16
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 3.122 | 3.135 | 3.147 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 0.8 -3.4905808 0 -2.7150906 -0.59565852 2016
|
||||
10000 0.82075861 -3.1822235 0 -2.3866107 0.013840263 2016
|
||||
20000 0.76467575 -3.0955084 0 -2.3542602 -0.076868925 2016
|
||||
30000 0.75803557 -3.1011543 0 -2.3663428 -0.052887049 2016
|
||||
40000 0.81732724 -3.064259 0 -2.2719724 0.070708808 2016
|
||||
50000 0.75874551 -3.1070261 0 -2.3715265 -0.074970431 2016
|
||||
Loop time of 22.1566 on 4 procs for 50000 steps with 1224 atoms
|
||||
|
||||
Performance: 974879.887 tau/day, 2256.666 timesteps/s
|
||||
98.5% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 10.732 | 13.012 | 14.39 | 38.0 | 58.73
|
||||
Neigh | 2.47 | 3.7351 | 4.3661 | 38.4 | 16.86
|
||||
Comm | 1.881 | 3.4383 | 5.8722 | 79.7 | 15.52
|
||||
Output | 0.00014567 | 0.0003581 | 0.0009892 | 0.0 | 0.00
|
||||
Modify | 1.1006 | 1.5188 | 2.6121 | 51.3 | 6.85
|
||||
Other | | 0.4521 | | | 2.04
|
||||
|
||||
Nlocal: 306 ave 312 max 295 min
|
||||
Histogram: 1 0 0 0 0 0 0 1 1 1
|
||||
Nghost: 1242.75 ave 1373 max 944 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 1 2
|
||||
Neighs: 9770.25 ave 10807 max 8736 min
|
||||
Histogram: 1 0 0 1 0 0 1 0 0 1
|
||||
|
||||
Total # of neighbors = 39081
|
||||
Ave neighs/atom = 31.9289
|
||||
Neighbor list builds = 5704
|
||||
Dangerous builds = 0
|
||||
# undump 1
|
||||
reset_timestep 0
|
||||
|
||||
compute bin_z liquid chunk/atom bin/1d z 0.0 0.1 units box
|
||||
|
||||
compute liquidStress_ke liquid stress/atom NULL ke
|
||||
compute liquidStress_vir liquid stress/atom NULL virial
|
||||
fix profile_z liquid ave/chunk 10 1000 10000 bin_z density/number temp c_liquidStress_ke[1] c_liquidStress_ke[2] c_liquidStress_ke[3] c_liquidStress_ke[4] c_liquidStress_ke[5] c_liquidStress_ke[6] c_liquidStress_vir[1] c_liquidStress_vir[2] c_liquidStress_vir[3] c_liquidStress_vir[4] c_liquidStress_vir[5] c_liquidStress_vir[6] ave running overwrite file profile.z
|
||||
|
||||
compute mopz0 all mop z center kin conf
|
||||
fix mopz0t all ave/time 10 1000 10000 c_mopz0[*] file mopz0.time
|
||||
|
||||
compute moppz liquid mop/profile z 0.0 0.1 kin conf
|
||||
fix moppzt all ave/time 100 100 10000 c_moppz[*] ave running overwrite file moppz.time mode vector
|
||||
|
||||
run 40000
|
||||
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 16
|
||||
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 mop, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(3) compute 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) = 4.139 | 4.147 | 4.15 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 0.75874551 -3.1070261 0 -2.3715265 -0.074970431 2016
|
||||
10000 0.82372476 -3.1299329 0 -2.3314448 -0.14706101 2016
|
||||
20000 0.80692892 -3.1278896 0 -2.3456828 -0.085123604 2016
|
||||
30000 0.78458951 -3.0966006 0 -2.3360488 0.13637007 2016
|
||||
40000 0.80106495 -3.1135836 0 -2.3370611 -0.14404185 2016
|
||||
Loop time of 31.4145 on 4 procs for 40000 steps with 1224 atoms
|
||||
|
||||
Performance: 550065.249 tau/day, 1273.299 timesteps/s
|
||||
92.6% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 10.199 | 12.307 | 13.428 | 35.4 | 39.18
|
||||
Neigh | 2.1261 | 3.1416 | 3.6373 | 33.5 | 10.00
|
||||
Comm | 3.5381 | 4.476 | 6.229 | 48.9 | 14.25
|
||||
Output | 0.00062943 | 0.0031546 | 0.0040004 | 2.6 | 0.01
|
||||
Modify | 10.186 | 10.862 | 12.26 | 24.8 | 34.58
|
||||
Other | | 0.6247 | | | 1.99
|
||||
|
||||
Nlocal: 306 ave 315 max 299 min
|
||||
Histogram: 1 0 0 1 1 0 0 0 0 1
|
||||
Nghost: 1221.5 ave 1347 max 912 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 1 2
|
||||
Neighs: 9710.25 ave 10301 max 8980 min
|
||||
Histogram: 1 0 0 0 0 1 1 0 0 1
|
||||
|
||||
Total # of neighbors = 38841
|
||||
Ave neighs/atom = 31.7328
|
||||
Neighbor list builds = 4573
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:00:53
|
||||
|
|
@ -288,10 +288,6 @@
|
|||
/compute_meso_rho_atom.h
|
||||
/compute_meso_t_atom.cpp
|
||||
/compute_meso_t_atom.h
|
||||
/compute_mop.cpp
|
||||
/compute_mop.h
|
||||
/compute_mop_profile.cpp
|
||||
/compute_mop_profile.h
|
||||
/compute_msd_nongauss.cpp
|
||||
/compute_msd_nongauss.h
|
||||
/compute_pe_tally.cpp
|
||||
|
|
@ -306,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
|
||||
|
|
|
|||
|
|
@ -1,424 +0,0 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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_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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeMop::ComputeMop(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 6) error->all(FLERR,"Illegal compute 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 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 mop command"); //break;
|
||||
|
||||
iarg++;
|
||||
}
|
||||
|
||||
// Error check
|
||||
// 3D only
|
||||
|
||||
if (domain->dimension < 3)
|
||||
error->all(FLERR, "Compute mop incompatible with simulation dimension");
|
||||
|
||||
// orthogonal simulation box
|
||||
if (domain->triclinic != 0)
|
||||
error->all(FLERR, "Compute 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 mop is out of bounds");
|
||||
|
||||
|
||||
// Initialize some variables
|
||||
|
||||
values_local = values_global = vector = NULL;
|
||||
|
||||
// this fix produces a global vector
|
||||
|
||||
memory->create(vector,nvalues,"mop:vector");
|
||||
memory->create(values_local,nvalues,"mop/spatial:values_local");
|
||||
memory->create(values_global,nvalues,"mop/spatial:values_global");
|
||||
size_vector = nvalues;
|
||||
|
||||
vector_flag = 1;
|
||||
extvector = 0;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeMop::~ComputeMop()
|
||||
{
|
||||
|
||||
delete [] which;
|
||||
|
||||
memory->destroy(values_local);
|
||||
memory->destroy(values_global);
|
||||
memory->destroy(vector);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMop::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 mop requires fixed simulation box
|
||||
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
|
||||
error->all(FLERR, "Compute 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 mop");
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all(FLERR,"Pair style does not support compute mop");
|
||||
|
||||
// Warnings
|
||||
|
||||
if (me==0){
|
||||
|
||||
//Compute 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 mop does not account for bond potentials");
|
||||
if (force->angle!=NULL)
|
||||
error->warning(FLERR,"compute mop does not account for angle potentials");
|
||||
if (force->dihedral!=NULL)
|
||||
error->warning(FLERR,"compute mop does not account for dihedral potentials");
|
||||
if (force->improper!=NULL)
|
||||
error->warning(FLERR,"compute mop does not account for improper potentials");
|
||||
if (force->kspace!=NULL)
|
||||
error->warning(FLERR,"compute 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 ComputeMop::init_list(int id, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute output vector
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMop::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 ComputeMop::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;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1,74 +0,0 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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(mop,ComputeMop)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_MOP_H
|
||||
#define LMP_COMPUTE_MOP_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeMop : public Compute {
|
||||
public:
|
||||
ComputeMop(class LAMMPS *, int, char **);
|
||||
virtual ~ComputeMop();
|
||||
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: Pair style does not support compute mop
|
||||
|
||||
The pair style does not have a single() function, so it can
|
||||
not be invoked by compute mop/spatial.
|
||||
|
||||
|
||||
|
||||
*/
|
||||
|
||||
|
|
@ -1,496 +0,0 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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_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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeMopProfile::ComputeMopProfile(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 7) error->all(FLERR,"Illegal compute 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 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 mop/profile command"); //break;
|
||||
|
||||
iarg++;
|
||||
}
|
||||
|
||||
// check domain related errors
|
||||
|
||||
// 3D only
|
||||
|
||||
if (domain->dimension < 3)
|
||||
error->all(FLERR, "Compute mop/profile incompatible with simulation dimension");
|
||||
|
||||
// orthogonal simulation box
|
||||
|
||||
if (domain->triclinic != 0)
|
||||
error->all(FLERR, "Compute 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,"mop/profile:array");
|
||||
size_array_rows = nbins;
|
||||
size_array_cols = 1 + nvalues;
|
||||
|
||||
array_flag = 1;
|
||||
extarray = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeMopProfile::~ComputeMopProfile()
|
||||
{
|
||||
|
||||
delete [] which;
|
||||
|
||||
memory->destroy(coord);
|
||||
memory->destroy(coordp);
|
||||
memory->destroy(values_local);
|
||||
memory->destroy(values_global);
|
||||
memory->destroy(array);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMopProfile::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 mop/profile requires fixed simulation box
|
||||
|
||||
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
|
||||
error->all(FLERR, "Compute 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 mop/profile");
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all(FLERR,"Pair style does not support compute mop/profile");
|
||||
|
||||
// Warnings
|
||||
|
||||
if (me==0){
|
||||
|
||||
//Compute 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 mop/profile does not account for bond potentials");
|
||||
if (force->angle!=NULL)
|
||||
error->warning(FLERR,"compute mop/profile does not account for angle potentials");
|
||||
if (force->dihedral!=NULL)
|
||||
error->warning(FLERR,"compute mop/profile does not account for dihedral potentials");
|
||||
if (force->improper!=NULL)
|
||||
error->warning(FLERR,"compute mop/profile does not account for improper potentials");
|
||||
if (force->kspace!=NULL)
|
||||
error->warning(FLERR,"compute 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 ComputeMopProfile::init_list(int id, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute output array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMopProfile::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 ComputeMopProfile::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 ComputeMopProfile::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 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 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,"mop/profile:coord");
|
||||
memory->create(coordp,nbins,1,"mop/profile:coordp");
|
||||
memory->create(values_local,nbins,nvalues,"mop/profile:values_local");
|
||||
memory->create(values_global,nbins,nvalues,"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];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,81 +0,0 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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(mop/profile,ComputeMopProfile)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_MOP_PROFILE_H
|
||||
#define LMP_COMPUTE_MOP_PROFILE_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeMopProfile : public Compute {
|
||||
public:
|
||||
ComputeMopProfile(class LAMMPS *, int, char **);
|
||||
virtual ~ComputeMopProfile();
|
||||
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: Pair style does not support compute mop/profile
|
||||
|
||||
The pair style does not have a single() function, so it can
|
||||
not be invoked by compute mop/profile.
|
||||
|
||||
|
||||
|
||||
*/
|
||||
|
||||
Loading…
Reference in New Issue