forked from lijiext/lammps
Merge pull request #925 from PabloPiaggi/pair_entropy
Pair entropy fingerprint
This commit is contained in:
Steve Plimpton
GitHub
commit
8e14143908
Binary file not shown.
|
After Width: | Height: | Size: 8.7 KiB |
|
|
@ -0,0 +1,10 @@
|
||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
\thispagestyle{empty}
|
||||||
|
|
||||||
|
$$
|
||||||
|
s_S^i=-2\pi\rho k_B \int\limits_0^{r_m} \left [ g(r) \ln g(r) - g(r) + 1 \right ] r^2 dr ,
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
Binary file not shown.
|
After Width: | Height: | Size: 8.3 KiB |
|
|
@ -0,0 +1,10 @@
|
||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
\thispagestyle{empty}
|
||||||
|
|
||||||
|
$$
|
||||||
|
g_m^i(r) = \frac{1}{4 \pi \rho r^2} \sum\limits_{j} \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(r-r_{ij})^2/(2\sigma^2)} ,
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
Binary file not shown.
|
After Width: | Height: | Size: 5.6 KiB |
|
|
@ -0,0 +1,10 @@
|
||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
\thispagestyle{empty}
|
||||||
|
|
||||||
|
$$
|
||||||
|
\bar{s}_S^i = \frac{\sum_j s_S^j + s_S^i}{N + 1} ,
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
|
|
@ -863,6 +863,7 @@ package"_Section_start.html#start_3.
|
||||||
"dpd"_compute_dpd.html,
|
"dpd"_compute_dpd.html,
|
||||||
"dpd/atom"_compute_dpd_atom.html,
|
"dpd/atom"_compute_dpd_atom.html,
|
||||||
"edpd/temp/atom"_compute_edpd_temp_atom.html,
|
"edpd/temp/atom"_compute_edpd_temp_atom.html,
|
||||||
|
"entropy/atom"_compute_entropy_atom.html,
|
||||||
"fep"_compute_fep.html,
|
"fep"_compute_fep.html,
|
||||||
"force/tally"_compute_tally.html,
|
"force/tally"_compute_tally.html,
|
||||||
"heat/flux/tally"_compute_tally.html,
|
"heat/flux/tally"_compute_tally.html,
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,130 @@
|
||||||
|
"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 entropy/atom command :h3
|
||||||
|
|
||||||
|
[Syntax:]
|
||||||
|
|
||||||
|
compute ID group-ID entropy/atom sigma cutoff keyword value ... :pre
|
||||||
|
|
||||||
|
ID, group-ID are documented in "compute"_compute.html command :l
|
||||||
|
entropy/atom = style name of this compute command :l
|
||||||
|
sigma = width of gaussians used in the g(r) smoothening :l
|
||||||
|
cutoff = cutoff for the g(r) calculation :l
|
||||||
|
one or more keyword/value pairs may be appended :l
|
||||||
|
keyword = {avg} or {local}
|
||||||
|
{avg} values = {yes} or {no} cutoff2
|
||||||
|
{yes} = average the pair entropy over neighbors
|
||||||
|
{no} = do not average the pair entropy over neighbors
|
||||||
|
cutoff2 = cutoff for the averaging over neighbors
|
||||||
|
{local} values = {yes} or {no} = use the local density around each atom to normalize the g(r) :pre
|
||||||
|
:ule
|
||||||
|
|
||||||
|
[Examples:]
|
||||||
|
|
||||||
|
compute 1 all entropy/atom 0.25 5.
|
||||||
|
compute 1 all entropy/atom 0.25 5. avg yes 5.
|
||||||
|
compute 1 all entropy/atom 0.125 7.3 avg yes 5.1 local yes :pre
|
||||||
|
|
||||||
|
[Description:]
|
||||||
|
|
||||||
|
Define a computation that calculates the pair entropy fingerprint for
|
||||||
|
each atom in the group. The fingerprint is useful to distinguish between
|
||||||
|
ordered and disordered environments, for instance liquid and solid-like
|
||||||
|
environments, or glassy and crystalline-like environments. Some
|
||||||
|
applications could be the identification of grain boundaries, a
|
||||||
|
melt-solid interface, or a solid cluster emerging from the melt.
|
||||||
|
The advantage of this parameter over others is that no a priori
|
||||||
|
information about the solid structure is required.
|
||||||
|
|
||||||
|
This parameter for atom i is computed using the following formula from
|
||||||
|
"(Piaggi)"_#Piaggi and "(Nettleton)"_#Nettleton ,
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_entropy.jpg)
|
||||||
|
|
||||||
|
where r is a distance, g(r) is the radial distribution function of atom
|
||||||
|
i and rho is the density of the system. The g(r) computed for each
|
||||||
|
atom i can be noisy and therefore it is smoothened using:
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_entropy2.jpg)
|
||||||
|
|
||||||
|
where the sum in j goes through the neighbors of atom i, and sigma is a
|
||||||
|
parameter to control the smoothening.
|
||||||
|
|
||||||
|
The input parameters are {sigma} the smoothening parameter, and the
|
||||||
|
{cutoff} for the calculation of g(r).
|
||||||
|
|
||||||
|
If the keyword {avg} has the setting {yes}, then this compute also
|
||||||
|
averages the parameter over the neighbors of atom i according to:
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_entropy3.jpg)
|
||||||
|
|
||||||
|
where the sum j goes over the neighbors of atom i and N is the number
|
||||||
|
of neighbors. This procedure provides a sharper distinction between
|
||||||
|
order and disorder environments. In this case the input parameter
|
||||||
|
{cutoff2} is the cutoff for the averaging over the neighbors and
|
||||||
|
must also be specified.
|
||||||
|
|
||||||
|
If the {avg yes} option is used, the effective cutoff of the neighbor
|
||||||
|
list should be {cutoff}+{cutoff2} and therefore it might be necessary
|
||||||
|
to increase the skin of the neighbor list with:
|
||||||
|
|
||||||
|
neighbor skin bin :pre
|
||||||
|
|
||||||
|
See "neighbor"_neighbor.html for details.
|
||||||
|
|
||||||
|
If the {local yes} option is used, the g(r) is normalized by the
|
||||||
|
local density around each atom, that is to say the density around each
|
||||||
|
atom is the number of neighbors within the neighbor list cutoff divided
|
||||||
|
by the corresponding volume. This option can be useful when dealing with
|
||||||
|
inhomogeneus systems such as those that have surfaces.
|
||||||
|
|
||||||
|
Here are typical input parameters for fcc aluminum (lattice
|
||||||
|
constant 4.05 Angstroms),
|
||||||
|
|
||||||
|
compute 1 all entropy/atom 0.25 5.7 avg yes 3.7 :pre
|
||||||
|
|
||||||
|
and for bcc sodium (lattice constant 4.23 Angstroms),
|
||||||
|
|
||||||
|
compute 1 all entropy/atom 0.25 7.3 avg yes 5.1 :pre
|
||||||
|
|
||||||
|
|
||||||
|
[Output info:]
|
||||||
|
|
||||||
|
By default, this compute calculates the pair entropy value for each
|
||||||
|
atom as a per-atom vector, which can be accessed by any command that
|
||||||
|
uses per-atom values from a compute as input. See "Section
|
||||||
|
6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output
|
||||||
|
options.
|
||||||
|
|
||||||
|
The pair entropy values have units of the Boltzmann constant. They are
|
||||||
|
always negative, and lower values (lower entropy) correspond to more
|
||||||
|
ordered environments.
|
||||||
|
|
||||||
|
[Restrictions:]
|
||||||
|
|
||||||
|
This compute is part of the USER-MISC package. It is only enabled if
|
||||||
|
LAMMPS was built with that package. See the "Making
|
||||||
|
LAMMPS"_Section_start.html#start_3 section for more info.
|
||||||
|
|
||||||
|
[Related commands:]
|
||||||
|
|
||||||
|
"compute cna/atom"_compute_cna_atom.html
|
||||||
|
"compute centro/atom"_compute_centro_atom.html
|
||||||
|
|
||||||
|
[Default:]
|
||||||
|
|
||||||
|
The default values for the optional keywords are avg = no and local = no.
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
:link(Piaggi)
|
||||||
|
[(Piaggi)] Piaggi and Parrinello, J Chem Phys, 147, 114112 (2017).
|
||||||
|
|
||||||
|
:link(Nettleton)
|
||||||
|
[(Nettleton)] Nettleton and Green, J Chem Phys, 29, 6 (1958).
|
||||||
|
|
@ -31,6 +31,7 @@ Computes :h1
|
||||||
compute_dpd
|
compute_dpd
|
||||||
compute_dpd_atom
|
compute_dpd_atom
|
||||||
compute_edpd_temp_atom
|
compute_edpd_temp_atom
|
||||||
|
compute_entropy_atom
|
||||||
compute_erotate_asphere
|
compute_erotate_asphere
|
||||||
compute_erotate_rigid
|
compute_erotate_rigid
|
||||||
compute_erotate_sphere
|
compute_erotate_sphere
|
||||||
|
|
|
||||||
|
|
@ -317,6 +317,7 @@ compute_displace_atom.html
|
||||||
compute_dpd.html
|
compute_dpd.html
|
||||||
compute_dpd_atom.html
|
compute_dpd_atom.html
|
||||||
compute_edpd_temp_atom.html
|
compute_edpd_temp_atom.html
|
||||||
|
compute_entropy_atom.html
|
||||||
compute_erotate_asphere.html
|
compute_erotate_asphere.html
|
||||||
compute_erotate_rigid.html
|
compute_erotate_rigid.html
|
||||||
compute_erotate_sphere.html
|
compute_erotate_sphere.html
|
||||||
|
|
|
||||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,34 @@
|
||||||
|
echo both
|
||||||
|
|
||||||
|
units metal
|
||||||
|
atom_style full
|
||||||
|
|
||||||
|
read_data data.interface
|
||||||
|
mass 1 22.98977
|
||||||
|
|
||||||
|
neigh_modify delay 10 every 1
|
||||||
|
pair_style eam/fs
|
||||||
|
pair_coeff * * Na_MendelevM_2014.eam.fs Na
|
||||||
|
timestep 0.002
|
||||||
|
thermo 500
|
||||||
|
|
||||||
|
neighbor 4. bin
|
||||||
|
|
||||||
|
# Define computes
|
||||||
|
# Global density, no average
|
||||||
|
compute 1 all entropy/atom 0.25 7.75
|
||||||
|
# Local density, no average
|
||||||
|
compute 2 all entropy/atom 0.25 7.75 local yes
|
||||||
|
# Global density, average over neighbors
|
||||||
|
compute 3 all entropy/atom 0.25 7.75 avg yes 5.
|
||||||
|
# Local density, average over neighbors
|
||||||
|
compute 4 all entropy/atom 0.25 7.75 avg yes 5. local yes
|
||||||
|
|
||||||
|
dump myDump all custom 500 dump.interface id type x y z c_1 c_2 c_3 c_4
|
||||||
|
|
||||||
|
|
||||||
|
fix 1 all nph x 1. 1. 10.
|
||||||
|
fix 2 all temp/csvr 350. 350. 0.1 64582
|
||||||
|
|
||||||
|
run 1000
|
||||||
|
|
||||||
|
|
@ -0,0 +1,124 @@
|
||||||
|
LAMMPS (30 Mar 2018)
|
||||||
|
|
||||||
|
units metal
|
||||||
|
atom_style full
|
||||||
|
|
||||||
|
read_data data.interface
|
||||||
|
Reading data file ...
|
||||||
|
triclinic box = (0 0 0) to (138.4 34.57 34.57) with tilt (0 0 0)
|
||||||
|
4 by 1 by 1 MPI processor grid
|
||||||
|
reading atoms ...
|
||||||
|
4096 atoms
|
||||||
|
reading velocities ...
|
||||||
|
4096 velocities
|
||||||
|
Finding 1-2 1-3 1-4 neighbors ...
|
||||||
|
special bond factors lj: 0 0 0
|
||||||
|
special bond factors coul: 0 0 0
|
||||||
|
0 = max # of 1-2 neighbors
|
||||||
|
0 = max # of 1-3 neighbors
|
||||||
|
0 = max # of 1-4 neighbors
|
||||||
|
1 = max # of special neighbors
|
||||||
|
mass 1 22.98977
|
||||||
|
|
||||||
|
neigh_modify delay 10 every 1
|
||||||
|
pair_style eam/fs
|
||||||
|
pair_coeff * * Na_MendelevM_2014.eam.fs Na
|
||||||
|
timestep 0.002
|
||||||
|
thermo 500
|
||||||
|
|
||||||
|
neighbor 4. bin
|
||||||
|
|
||||||
|
# Define computes
|
||||||
|
# Global density, no average
|
||||||
|
compute 1 all entropy/atom 0.25 7.75
|
||||||
|
# Local density, no average
|
||||||
|
compute 2 all entropy/atom 0.25 7.75 local yes
|
||||||
|
# Global density, average over neighbors
|
||||||
|
compute 3 all entropy/atom 0.25 7.75 avg yes 5.
|
||||||
|
# Local density, average over neighbors
|
||||||
|
compute 4 all entropy/atom 0.25 7.75 avg yes 5. local yes
|
||||||
|
|
||||||
|
dump myDump all custom 500 dump.interface id type x y z c_1 c_2 c_3 c_4
|
||||||
|
|
||||||
|
|
||||||
|
fix 1 all nph x 1. 1. 10.
|
||||||
|
fix 2 all temp/csvr 350. 350. 0.1 64582
|
||||||
|
|
||||||
|
run 1000
|
||||||
|
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||||
|
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||||
|
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||||
|
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||||
|
Neighbor list info ...
|
||||||
|
update every 1 steps, delay 10 steps, check yes
|
||||||
|
max neighbors/atom: 2000, page size: 100000
|
||||||
|
master list distance cutoff = 13.2
|
||||||
|
ghost atom cutoff = 13.2
|
||||||
|
binsize = 6.6, bins = 21 6 6
|
||||||
|
5 neighbor lists, perpetual/occasional/extra = 5 0 0
|
||||||
|
(1) pair eam/fs, perpetual
|
||||||
|
attributes: half, newton on
|
||||||
|
pair build: half/bin/newton/tri
|
||||||
|
stencil: half/bin/3d/newton/tri
|
||||||
|
bin: standard
|
||||||
|
(2) compute entropy/atom, perpetual
|
||||||
|
attributes: full, newton on, ghost
|
||||||
|
pair build: full/bin/ghost
|
||||||
|
stencil: full/ghost/bin/3d
|
||||||
|
bin: standard
|
||||||
|
(3) compute entropy/atom, perpetual, copy from (2)
|
||||||
|
attributes: full, newton on, ghost
|
||||||
|
pair build: copy
|
||||||
|
stencil: none
|
||||||
|
bin: none
|
||||||
|
(4) compute entropy/atom, perpetual, copy from (2)
|
||||||
|
attributes: full, newton on, ghost
|
||||||
|
pair build: copy
|
||||||
|
stencil: none
|
||||||
|
bin: none
|
||||||
|
(5) compute entropy/atom, perpetual, copy from (2)
|
||||||
|
attributes: full, newton on, ghost
|
||||||
|
pair build: copy
|
||||||
|
stencil: none
|
||||||
|
bin: none
|
||||||
|
Setting up Verlet run ...
|
||||||
|
Unit style : metal
|
||||||
|
Current step : 0
|
||||||
|
Time step : 0.002
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 25.68 | 25.69 | 25.69 Mbytes
|
||||||
|
Step Temp E_pair E_mol TotEng Press Volume
|
||||||
|
0 346.29871 -4285.222 0 -4101.9191 594.65353 165399.75
|
||||||
|
500 359.33758 -4285.247 0 -4095.0423 471.98587 165847.18
|
||||||
|
1000 348.99659 -4276.2274 0 -4091.4964 149.27188 166966.18
|
||||||
|
Loop time of 5.3437 on 4 procs for 1000 steps with 4096 atoms
|
||||||
|
|
||||||
|
Performance: 32.337 ns/day, 0.742 hours/ns, 187.136 timesteps/s
|
||||||
|
99.8% CPU use with 4 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 4.2832 | 4.3257 | 4.3839 | 1.8 | 80.95
|
||||||
|
Bond | 0.00018309 | 0.00019825 | 0.00021418 | 0.0 | 0.00
|
||||||
|
Neigh | 0.42195 | 0.42512 | 0.42739 | 0.3 | 7.96
|
||||||
|
Comm | 0.051679 | 0.1101 | 0.14916 | 10.8 | 2.06
|
||||||
|
Output | 0.40909 | 0.4091 | 0.40911 | 0.0 | 7.66
|
||||||
|
Modify | 0.060869 | 0.061921 | 0.06327 | 0.4 | 1.16
|
||||||
|
Other | | 0.01161 | | | 0.22
|
||||||
|
|
||||||
|
Nlocal: 1024 ave 1040 max 1001 min
|
||||||
|
Histogram: 1 0 0 0 0 0 2 0 0 1
|
||||||
|
Nghost: 4614.25 ave 4700 max 4540 min
|
||||||
|
Histogram: 1 1 0 0 0 0 0 1 0 1
|
||||||
|
Neighs: 121747 ave 126398 max 116931 min
|
||||||
|
Histogram: 1 0 0 1 0 0 1 0 0 1
|
||||||
|
FullNghs: 243494 ave 252523 max 233842 min
|
||||||
|
Histogram: 1 0 0 1 0 0 1 0 0 1
|
||||||
|
|
||||||
|
Total # of neighbors = 973974
|
||||||
|
Ave neighs/atom = 237.787
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 13
|
||||||
|
Dangerous builds = 0
|
||||||
|
|
||||||
|
Total wall time: 0:00:05
|
||||||
|
|
@ -29,6 +29,7 @@ bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||||
compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
|
compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
|
||||||
compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
|
compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
|
||||||
compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017
|
compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017
|
||||||
|
compute entropy/atom, Pablo Piaggi (EPFL), pablo.piaggi at epfl.ch, 15 June 2018
|
||||||
compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
||||||
compute PRESSURE/GREM, David Stelter, dstelter@bu.edu, 22 Nov 16
|
compute PRESSURE/GREM, David Stelter, dstelter@bu.edu, 22 Nov 16
|
||||||
dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,340 @@
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
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 author: Pablo Piaggi (EPFL Lausanne)
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include <cmath>
|
||||||
|
#include <cstring>
|
||||||
|
#include <cstdlib>
|
||||||
|
#include "compute_entropy_atom.h"
|
||||||
|
#include "atom.h"
|
||||||
|
#include "update.h"
|
||||||
|
#include "modify.h"
|
||||||
|
#include "neighbor.h"
|
||||||
|
#include "neigh_list.h"
|
||||||
|
#include "neigh_request.h"
|
||||||
|
#include "force.h"
|
||||||
|
#include "pair.h"
|
||||||
|
#include "comm.h"
|
||||||
|
#include "math_extra.h"
|
||||||
|
#include "math_const.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include "error.h"
|
||||||
|
#include "domain.h"
|
||||||
|
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
using namespace MathConst;
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
ComputeEntropyAtom::
|
||||||
|
ComputeEntropyAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||||
|
Compute(lmp, narg, arg),
|
||||||
|
pair_entropy(NULL), pair_entropy_avg(NULL)
|
||||||
|
{
|
||||||
|
if (narg < 5 || narg > 10)
|
||||||
|
error->all(FLERR,"Illegal compute entropy/atom command; wrong number"
|
||||||
|
" of arguments");
|
||||||
|
|
||||||
|
// Arguments are: sigma cutoff avg yes/no cutoff2 local yes/no
|
||||||
|
// sigma is the gaussian width
|
||||||
|
// cutoff is the cutoff for the calculation of g(r)
|
||||||
|
// avg is optional and allows averaging the pair entropy over neighbors
|
||||||
|
// the next argument should be yes or no
|
||||||
|
// cutoff2 is the cutoff for the averaging
|
||||||
|
// local is optional and allows using the local density to normalize
|
||||||
|
// the g(r)
|
||||||
|
|
||||||
|
sigma = force->numeric(FLERR,arg[3]);
|
||||||
|
if (sigma < 0.0) error->all(FLERR,"Illegal compute entropy/atom"
|
||||||
|
" command; negative sigma");
|
||||||
|
cutoff = force->numeric(FLERR,arg[4]);
|
||||||
|
if (cutoff < 0.0) error->all(FLERR,"Illegal compute entropy/atom"
|
||||||
|
" command; negative cutoff");
|
||||||
|
|
||||||
|
avg_flag = 0;
|
||||||
|
local_flag = 0;
|
||||||
|
|
||||||
|
// optional keywords
|
||||||
|
int iarg = 5;
|
||||||
|
while (iarg < narg) {
|
||||||
|
if (strcmp(arg[iarg],"avg") == 0) {
|
||||||
|
if (iarg+2 > narg)
|
||||||
|
error->all(FLERR,"Illegal compute entropy/atom;"
|
||||||
|
" missing arguments after avg");
|
||||||
|
if (strcmp(arg[iarg+1],"yes") == 0) avg_flag = 1;
|
||||||
|
else if (strcmp(arg[iarg+1],"no") == 0) avg_flag = 0;
|
||||||
|
else error->all(FLERR,"Illegal compute entropy/atom;"
|
||||||
|
" argument after avg should be yes or no");
|
||||||
|
cutoff2 = force->numeric(FLERR,arg[iarg+2]);
|
||||||
|
if (cutoff2 < 0.0) error->all(FLERR,"Illegal compute entropy/atom"
|
||||||
|
" command; negative cutoff2");
|
||||||
|
cutsq2 = cutoff2*cutoff2;
|
||||||
|
iarg += 3;
|
||||||
|
} else if (strcmp(arg[iarg],"local") == 0) {
|
||||||
|
if (strcmp(arg[iarg+1],"yes") == 0) local_flag = 1;
|
||||||
|
else if (strcmp(arg[iarg+1],"no") == 0) local_flag = 0;
|
||||||
|
else error->all(FLERR,"Illegal compute entropy/atom;"
|
||||||
|
" argument after local should be yes or no");
|
||||||
|
iarg += 2;
|
||||||
|
} else error->all(FLERR,"Illegal compute entropy/atom; argument after"
|
||||||
|
" sigma and cutoff should be avg or local");
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
cutsq = cutoff*cutoff;
|
||||||
|
nbin = static_cast<int>(cutoff / sigma) + 1;
|
||||||
|
nmax = 0;
|
||||||
|
maxneigh = 0;
|
||||||
|
// Number of bins above and below the central one that will be
|
||||||
|
// considered as affected by the gaussian kernel
|
||||||
|
// 3 seems a good compromise between speed and good mollification
|
||||||
|
deltabin = 3;
|
||||||
|
deltar = sigma;
|
||||||
|
peratom_flag = 1;
|
||||||
|
size_peratom_cols = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
ComputeEntropyAtom::~ComputeEntropyAtom()
|
||||||
|
{
|
||||||
|
memory->destroy(pair_entropy);
|
||||||
|
if (avg_flag) memory->destroy(pair_entropy_avg);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void ComputeEntropyAtom::init()
|
||||||
|
{
|
||||||
|
if (force->pair == NULL)
|
||||||
|
error->all(FLERR,"Compute centro/atom requires a pair style be"
|
||||||
|
" defined");
|
||||||
|
|
||||||
|
if ( (cutoff+cutoff2) > (force->pair->cutforce + neighbor->skin) )
|
||||||
|
{
|
||||||
|
error->all(FLERR,"Compute entropy/atom cutoff is longer than the"
|
||||||
|
" pairwise cutoff. Increase the neighbor list skin"
|
||||||
|
" distance.");
|
||||||
|
}
|
||||||
|
|
||||||
|
int count = 0;
|
||||||
|
for (int i = 0; i < modify->ncompute; i++)
|
||||||
|
if (strcmp(modify->compute[i]->style,"entropy/atom") == 0) count++;
|
||||||
|
if (count > 1 && comm->me == 0)
|
||||||
|
error->warning(FLERR,"More than one compute entropy/atom");
|
||||||
|
|
||||||
|
// need a full neighbor list with neighbors of the ghost atoms
|
||||||
|
|
||||||
|
int irequest = neighbor->request(this,instance_me);
|
||||||
|
neighbor->requests[irequest]->pair = 0;
|
||||||
|
neighbor->requests[irequest]->compute = 1;
|
||||||
|
neighbor->requests[irequest]->half = 0;
|
||||||
|
neighbor->requests[irequest]->full = 1;
|
||||||
|
neighbor->requests[irequest]->occasional = 0;
|
||||||
|
neighbor->requests[irequest]->ghost = 1;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void ComputeEntropyAtom::init_list(int id, NeighList *ptr)
|
||||||
|
{
|
||||||
|
list = ptr;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void ComputeEntropyAtom::compute_peratom()
|
||||||
|
{
|
||||||
|
int i,j,ii,jj,inum,jnum;
|
||||||
|
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
|
||||||
|
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||||
|
double rbin[nbin], rbinsq[nbin];
|
||||||
|
|
||||||
|
invoked_peratom = update->ntimestep;
|
||||||
|
|
||||||
|
// Initialize distance vectors
|
||||||
|
for (int i = 0; i < nbin; i++) {
|
||||||
|
rbin[i] = i*deltar;
|
||||||
|
rbinsq[i] = rbin[i]*rbin[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
// grow pair_entropy and pair_entropy_avg array if necessary
|
||||||
|
|
||||||
|
if (atom->nmax > nmax) {
|
||||||
|
if (!avg_flag) {
|
||||||
|
memory->destroy(pair_entropy);
|
||||||
|
nmax = atom->nmax;
|
||||||
|
memory->create(pair_entropy,nmax,"entropy/atom:pair_entropy");
|
||||||
|
vector_atom = pair_entropy;
|
||||||
|
} else {
|
||||||
|
memory->destroy(pair_entropy);
|
||||||
|
memory->destroy(pair_entropy_avg);
|
||||||
|
nmax = atom->nmax;
|
||||||
|
memory->create(pair_entropy,nmax,"entropy/atom:pair_entropy");
|
||||||
|
memory->create(pair_entropy_avg,nmax,
|
||||||
|
"entropy/atom:pair_entropy_avg");
|
||||||
|
vector_atom = pair_entropy_avg;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
inum = list->inum + list->gnum;
|
||||||
|
ilist = list->ilist;
|
||||||
|
numneigh = list->numneigh;
|
||||||
|
firstneigh = list->firstneigh;
|
||||||
|
|
||||||
|
// Compute some constants
|
||||||
|
double sigmasq2=2*sigma*sigma;
|
||||||
|
double volume = domain->xprd * domain->yprd * domain->zprd;
|
||||||
|
double density = atom->natoms / volume;
|
||||||
|
|
||||||
|
// compute pair entropy for each atom in group
|
||||||
|
// use full neighbor list
|
||||||
|
|
||||||
|
double **x = atom->x;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
|
||||||
|
for (ii = 0; ii < inum; ii++) {
|
||||||
|
i = ilist[ii];
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
xtmp = x[i][0];
|
||||||
|
ytmp = x[i][1];
|
||||||
|
ztmp = x[i][2];
|
||||||
|
jlist = firstneigh[i];
|
||||||
|
jnum = numneigh[i];
|
||||||
|
|
||||||
|
// If local density is used, calculate it
|
||||||
|
if (local_flag) {
|
||||||
|
double neigh_cutoff = force->pair->cutforce + neighbor->skin;
|
||||||
|
double volume =
|
||||||
|
(4./3.)*MY_PI*neigh_cutoff*neigh_cutoff*neigh_cutoff;
|
||||||
|
density = jnum / volume;
|
||||||
|
}
|
||||||
|
|
||||||
|
// calculate kernel normalization
|
||||||
|
// Normalization of g(r)
|
||||||
|
double normConstantBase = 4*MY_PI*density;
|
||||||
|
// Normalization of gaussian
|
||||||
|
normConstantBase *= sqrt(2.*MY_PI)*sigma;
|
||||||
|
double invNormConstantBase = 1./normConstantBase;
|
||||||
|
|
||||||
|
// loop over list of all neighbors within force cutoff
|
||||||
|
|
||||||
|
// initialize gofr
|
||||||
|
double gofr[nbin];
|
||||||
|
for(int k=0;k<nbin;++k) gofr[k]=0.;
|
||||||
|
|
||||||
|
for (jj = 0; jj < jnum; jj++) {
|
||||||
|
j = jlist[jj];
|
||||||
|
j &= NEIGHMASK;
|
||||||
|
|
||||||
|
delx = xtmp - x[j][0];
|
||||||
|
dely = ytmp - x[j][1];
|
||||||
|
delz = ztmp - x[j][2];
|
||||||
|
rsq = delx*delx + dely*dely + delz*delz;
|
||||||
|
if (rsq < cutsq) {
|
||||||
|
// contribute to gofr
|
||||||
|
double r=sqrt(rsq);
|
||||||
|
int bin=floor(r/deltar);
|
||||||
|
int minbin, maxbin;
|
||||||
|
minbin=bin - deltabin;
|
||||||
|
if (minbin < 0) minbin=0;
|
||||||
|
if (minbin > (nbin-1)) minbin=nbin-1;
|
||||||
|
maxbin=bin + deltabin;
|
||||||
|
if (maxbin > (nbin-1)) maxbin=nbin-1;
|
||||||
|
for(int k=minbin;k<maxbin+1;k++) {
|
||||||
|
double invNormKernel=invNormConstantBase/rbinsq[k];
|
||||||
|
double distance = r - rbin[k];
|
||||||
|
gofr[k] += invNormKernel*exp(-distance*distance/sigmasq2);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Calculate integrand
|
||||||
|
double integrand[nbin];
|
||||||
|
for(int k=0;k<nbin;++k){
|
||||||
|
if (gofr[k]<1.e-10) {
|
||||||
|
integrand[k] = rbinsq[k];
|
||||||
|
} else {
|
||||||
|
integrand[k] = (gofr[k]*log(gofr[k])-gofr[k]+1)*rbinsq[k];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Integrate with trapezoid rule
|
||||||
|
double value = 0.;
|
||||||
|
for(int k=1;k<nbin-1;++k){
|
||||||
|
value += integrand[k];
|
||||||
|
}
|
||||||
|
value += 0.5*integrand[0];
|
||||||
|
value += 0.5*integrand[nbin-1];
|
||||||
|
value *= deltar;
|
||||||
|
|
||||||
|
pair_entropy[i] = -2*MY_PI*density*value;
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
if (avg_flag) {
|
||||||
|
for (ii = 0; ii < inum; ii++) {
|
||||||
|
i = ilist[ii];
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
xtmp = x[i][0];
|
||||||
|
ytmp = x[i][1];
|
||||||
|
ztmp = x[i][2];
|
||||||
|
jlist = firstneigh[i];
|
||||||
|
jnum = numneigh[i];
|
||||||
|
|
||||||
|
pair_entropy_avg[i] = pair_entropy[i];
|
||||||
|
double counter = 1;
|
||||||
|
|
||||||
|
// loop over list of all neighbors within force cutoff
|
||||||
|
for (jj = 0; jj < jnum; jj++) {
|
||||||
|
j = jlist[jj];
|
||||||
|
j &= NEIGHMASK;
|
||||||
|
|
||||||
|
delx = xtmp - x[j][0];
|
||||||
|
dely = ytmp - x[j][1];
|
||||||
|
delz = ztmp - x[j][2];
|
||||||
|
rsq = delx*delx + dely*dely + delz*delz;
|
||||||
|
if (rsq < cutsq2) {
|
||||||
|
pair_entropy_avg[i] += pair_entropy[j];
|
||||||
|
counter += 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
pair_entropy_avg[i] /= counter;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
memory usage of local atom-based array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
double ComputeEntropyAtom::memory_usage()
|
||||||
|
{
|
||||||
|
double bytes;
|
||||||
|
if (!avg_flag) {
|
||||||
|
bytes = nmax * sizeof(double);
|
||||||
|
} else {
|
||||||
|
bytes = 2 * nmax * sizeof(double);
|
||||||
|
}
|
||||||
|
return bytes;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,71 @@
|
||||||
|
/* -*- c++ -*- ----------------------------------------------------------
|
||||||
|
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.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#ifdef COMPUTE_CLASS
|
||||||
|
|
||||||
|
ComputeStyle(entropy/atom,ComputeEntropyAtom)
|
||||||
|
|
||||||
|
#else
|
||||||
|
|
||||||
|
#ifndef COMPUTE_ENTROPY_ATOM_H
|
||||||
|
#define COMPUTE_ENTROPY_ATOM_H
|
||||||
|
|
||||||
|
#include "compute.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class ComputeEntropyAtom : public Compute {
|
||||||
|
public:
|
||||||
|
ComputeEntropyAtom(class LAMMPS *, int, char **);
|
||||||
|
~ComputeEntropyAtom();
|
||||||
|
void init();
|
||||||
|
void init_list(int, class NeighList *);
|
||||||
|
void compute_peratom();
|
||||||
|
double memory_usage();
|
||||||
|
|
||||||
|
private:
|
||||||
|
int nmax,maxneigh, nbin;
|
||||||
|
class NeighList *list;
|
||||||
|
double *pair_entropy, *pair_entropy_avg;
|
||||||
|
double sigma, cutoff, cutoff2;
|
||||||
|
double cutsq, cutsq2;
|
||||||
|
double deltar;
|
||||||
|
int deltabin;
|
||||||
|
double invNormConstantBase;
|
||||||
|
int avg_flag;
|
||||||
|
int local_flag;
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#endif
|
||||||
|
|
||||||
|
/* ERROR/WARNING messages:
|
||||||
|
|
||||||
|
E: Illegal ... command
|
||||||
|
|
||||||
|
Self-explanatory. Check the input script syntax and compare to the
|
||||||
|
documentation for the command. You can use -echo screen as a
|
||||||
|
command-line option when running LAMMPS to see the offending line.
|
||||||
|
|
||||||
|
E: Compute entropy/atom requires a pair style be defined
|
||||||
|
|
||||||
|
This is because the computation of the pair entropy values
|
||||||
|
uses a pairwise neighbor list.
|
||||||
|
|
||||||
|
W: More than one compute entropy/atom
|
||||||
|
|
||||||
|
It is not efficient to use compute entropy/atom more than once.
|
||||||
|
|
||||||
|
*/
|
||||||
Loading…
Reference in New Issue