Examples - new local option

doc/src/compute_pair_entropy_atom.txt

@@ -42,7 +42,7 @@ 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
+"(Piaggi)"_#Piaggi and "(Nettleton)"_#Nettleton
 
 :c,image(Eqs/pair_entropy.jpg)
 
@@ -119,3 +119,6 @@ The default value for the optional keyword is avg = no.
 
 :link(Piaggi)
 [(Piaggi)] Piaggi and Parrinello, J Chem Phys, 147, 114112 (2017).
+
+:link(Nettleton}
+[(Nettleton)] Nettleton and Green, J Chem Phys, 29, 6 (1958).

examples/USER/misc/pair_entropy/start.lmp

@@ -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 pentropy/atom 0.25 7.75
+# Local density, no average
+compute         2 all pentropy/atom 0.25 7.75 local yes
+# Global density, average over neighbors
+compute         3 all pentropy/atom 0.25 7.75 avg yes 5.
+# Local density, average over neighbors
+compute         4 all pentropy/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             100000
+

src/compute_pair_entropy_atom.cpp

@@ -44,7 +44,7 @@ ComputePairEntropyAtom::ComputePairEntropyAtom(LAMMPS *lmp, int narg, char **arg
   Compute(lmp, narg, arg),
   pair_entropy(NULL), pair_entropy_avg(NULL)
 {
-  if (narg < 5 || narg > 8)
+  if (narg < 5 || narg > 10)
     error->all(FLERR,"Illegal compute pentropy/atom command");
 
   // Arguments are: sigma cutoff avg yes/no cutoff2
@@ -59,6 +59,7 @@ ComputePairEntropyAtom::ComputePairEntropyAtom(LAMMPS *lmp, int narg, char **arg
   if (cutoff < 0.0) error->all(FLERR,"Illegal compute pentropy/atom command; negative cutoff");
 
   avg_flag = 0;
+  local_flag = 0;
 
   // optional keywords
   int iarg = 5;
@@ -73,7 +74,12 @@ ComputePairEntropyAtom::ComputePairEntropyAtom(LAMMPS *lmp, int narg, char **arg
       if (cutoff2 < 0.0) error->all(FLERR,"Illegal compute pentropy/atom command; negative cutoff2");
       cutsq2 = cutoff2*cutoff2;
       iarg += 3;
-    } else error->all(FLERR,"Illegal compute pentropy/atom argument after sigma and cutoff should be avg");
+    } 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 pentropy/atom argument after local should be yes or no");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute pentropy/atom argument after sigma and cutoff should be avg or local");
   }
 
 
@@ -209,7 +215,13 @@ void ComputePairEntropyAtom::compute_peratom()
       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
       double normConstantBase = 4*MY_PI*density; // Normalization of g(r)
       normConstantBase *= sqrt(2.*MY_PI)*sigma; // Normalization of gaussian

src/compute_pair_entropy_atom.h

@@ -43,6 +43,7 @@ class ComputePairEntropyAtom : public Compute {
   int deltabin;
   double invNormConstantBase;
   int avg_flag;
+  int local_flag;
 };
 
 }