diff --git a/doc/src/Packages_details.txt b/doc/src/Packages_details.txt index 16f867b7c0..36071064c0 100644 --- a/doc/src/Packages_details.txt +++ b/doc/src/Packages_details.txt @@ -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 diff --git a/doc/src/Packages_user.txt b/doc/src/Packages_user.txt index 731cf7ff10..5b9b8734b3 100644 --- a/doc/src/Packages_user.txt +++ b/doc/src/Packages_user.txt @@ -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 diff --git a/doc/src/compute_mop.txt b/doc/src/compute_mop.txt deleted file mode 100644 index 732f55657a..0000000000 --- a/doc/src/compute_mop.txt +++ /dev/null @@ -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). diff --git a/doc/src/computes.txt b/doc/src/computes.txt index c431f1c184..158389bda9 100644 --- a/doc/src/computes.txt +++ b/doc/src/computes.txt @@ -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 diff --git a/doc/src/lammps.book b/doc/src/lammps.book index 0aa382ae39..c177757695 100644 --- a/doc/src/lammps.book +++ b/doc/src/lammps.book @@ -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 diff --git a/examples/USER/mop/in.compute_mop b/examples/USER/mop/in.compute_mop deleted file mode 100644 index ab9552796f..0000000000 --- a/examples/USER/mop/in.compute_mop +++ /dev/null @@ -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 - diff --git a/examples/USER/mop/in.compute_stress_mop b/examples/USER/mop/in.compute_stress_mop old mode 100755 new mode 100644 diff --git a/examples/USER/mop/log.31Aug18.compute_mop.g++.1 b/examples/USER/mop/log.31Aug18.compute_mop.g++.1 deleted file mode 100644 index 98b403b444..0000000000 --- a/examples/USER/mop/log.31Aug18.compute_mop.g++.1 +++ /dev/null @@ -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 diff --git a/examples/USER/mop/log.31Aug18.compute_mop.g++.4 b/examples/USER/mop/log.31Aug18.compute_mop.g++.4 deleted file mode 100644 index 16ef8a2ce1..0000000000 --- a/examples/USER/mop/log.31Aug18.compute_mop.g++.4 +++ /dev/null @@ -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 diff --git a/src/.gitignore b/src/.gitignore index c7aaa8ab64..3a6ba21426 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/USER-MOP/compute_mop.cpp b/src/USER-MOP/compute_mop.cpp deleted file mode 100644 index ba6b458172..0000000000 --- a/src/USER-MOP/compute_mop.cpp +++ /dev/null @@ -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 -#include -#include -#include - -#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 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; mmass; - 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 (mx[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]pos1) && (xj[dir]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)) || ((xi[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]pos1) && (xj[dir]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) -#include -#include -#include - -#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; ibinmass; - 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 (mx[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;ibinpos) && (xj[dir]pos1) && (xj[dir]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)) - || ((xi[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;ibinpos) && (xj[dir]pos1) && (xj[dir]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;ibinboxlo; - 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 ((origin - boxlo[dir]) * invdelta); - lo = origin - n*delta; - } - if (origin < boxhi[dir]) { - n = static_cast ((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 ((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]; - } - } -} diff --git a/src/USER-MOP/compute_mop_profile.h b/src/USER-MOP/compute_mop_profile.h deleted file mode 100644 index 3cb46087ae..0000000000 --- a/src/USER-MOP/compute_mop_profile.h +++ /dev/null @@ -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. - - - -*/ -