more changes to ATM source and doc file

doc/src/Section_commands.txt

@@ -582,7 +582,7 @@ USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT.
 "dt/reset"_fix_dt_reset.html,
 "efield"_fix_efield.html,
 "ehex"_fix_ehex.html,
-"enforce2d"_fix_enforce2d.html,
+"enforce2d (k)"_fix_enforce2d.html,
 "evaporate"_fix_evaporate.html,
 "external"_fix_external.html,
 "freeze"_fix_freeze.html,
@@ -929,6 +929,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "adp (o)"_pair_adp.html,
 "airebo (oi)"_pair_airebo.html,
 "airebo/morse (oi)"_pair_airebo.html,
+"atm"_pair_atm.html,
 "beck (go)"_pair_beck.html,
 "body"_pair_body.html,
 "bop"_pair_bop.html,

doc/src/pair_atm.txt

@@ -16,7 +16,7 @@ cutoff = global cutoff for 3-body interactions (distance units) :ul
 
 [Examples:]
 
-pair_style 2.5
+pair_style atm 2.5
 pair_coeff * * * 0.072 :pre
 
 pair_style hybrid/overlay lj/cut 6.5 atm 2.5
@@ -33,106 +33,92 @@ potential for the energy E of a system of atoms as
 
 :c,image(Eqs/pair_atm.jpg)
 
-where nu is the three-body interaction strength,
+where nu is the three-body interaction strength.  The distances
-and the distances between pairs of atoms r12, r23 and r31
+between pairs of atoms r12, r23, r31 and the angles gamma1, gamma2,
-and the angles gamma1, gamma2 and gamma3
+gamma3 are as shown in this diagram:
-are shown at the diagram:
 
 :c,image(JPG/pair_atm_dia.jpg)
 
-There is no \"central\" atom, the interaction is symmetric with respect
+Note that for the interaction between a triplet of atoms I,J,K, there
-to permutation of atom types.
+is no "central" atom.  The interaction is symmetric with respect to
+permutation of the three atoms.  In fact, it is computed 3 times by
+the code with each of I,J,K being atom 1.  Thus the nu value must be
+the same for all those permutations of the atom types of I,J,K as
+discussed below.
 
 The {atm} potential is typically used in combination with a two-body
 potential using the "pair_style hybrid/overlay"_pair_hybrid.html
 command as in the example above.
 
-The potential is calculated if all three atoms are in the
+The potential for a triplet of atom is calculated only if all 3
-"neighbor list"_neighbor.html
+distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff.
-and the distances between atoms satisfy r12 r23 r31 > cutoff^3.
 
 The following coefficients must be defined for each pair of atoms
 types via the "pair_coeff"_pair_coeff.html command as in the examples
 above, or in the restart files read by the
 "read_restart"_read_restart.html commands:
 
-K - the type of the third atom
+K = atom type of the third atom (1 to Ntypes)
-nu (energy/distance^9 units) :ul
+nu = prefactor (energy/distance^9 units) :ul
 
-For a single-atom type simulation, only a single entry is required, eg
+K can be specified in one of two ways.  An explicit numeric value can
+be used, as in the 2nd example above.  J <= K is required.  LAMMPS
+sets the coefficients for the other 5 symmetric interactions to the
+same values.  E.g. if I = 1, J = 2, K = 3, then these 6 values are set
+to the specified nu: nu123, nu132, nu213, nu231, nu312, nu321.  This
+enforces the symmetry discussed above.
 
-pair_coeff 1 1 1 nu :pre
+A wildcard asterisk can be used for K to set the coefficients for
+multiple triplets of atom types.  This takes the form "*" or "*n" or
+"n*" or "m*n".  If N = the number of atom types, then an asterisk with
+no numeric values means all types from 1 to N.  A leading asterisk
+means all types from 1 to n (inclusive).  A trailing asterisk means
+all types from n to N (inclusive).  A middle asterisk means all types
+from m to n (inclusive).  Note that only type triplets with J <= K are
+considered; if asterisks imply type triplets where K < J, they are
+ignored.
 
-where all three atoms are of type 1.
+Note that a pair_coeff command can override a previous setting for the
-For a two-atom type simulation, more pair_coeff commands can be used.
+same I,J,K triplet.  For example, these commands set nu for all I,J.K
-ATM interaction is symmetric with respect to permutation of atoms,
+triplets, then overwrite nu for just the I,J,K = 2,3,4 triplet:
-it is necessary to provide pair_coeff command for one permutation.
-For example, the command
 
-pair_coeff 1 1 2 nu :pre
+pair_coeff * * * 0.25
+pair_coeff 2 3 4 0.1 :pre
 
-also implies
+Note that for a simulation with a single atom type, only a single
+entry is required, e.g.
 
-pair_coeff 1 2 1 nu
+pair_coeff 1 1 1 0.25 :per
-pair_coeff 2 1 1 nu :pre
 
-Thus, to specify all ATM interactions between two atom types (eg 1 and 2)
+For a simulation with two atom types, four pair_coeff commands will
-it is sufficient to provide four pair_coeff commands, eg:
+specify all possible nu values:
 
 pair_coeff 1 1 1 nu1
 pair_coeff 1 1 2 nu2
 pair_coeff 1 2 2 nu3
 pair_coeff 2 2 2 nu4 :pre
 
-For 3 atom types, 10 pair_coeff commands
+For a simulation with three atom types, ten pair_coeff commands will
-(eg 111, 112, 113, 122, 123, 133, 222, 223, 233, 333)
+specify all possible nu values:
-will describe all possible ATM interactions, etc.
-It is not necessary to provide the pair_coeff commands for interactions
-of all combinations of atom types: if some combination is not defined then
-there is no ATM interaction for this combination and all its permutations.
 
---------------------
+pair_coeff 1 1 1 nu1
+pair_coeff 1 1 2 nu2
+pair_coeff 1 1 3 nu3
+pair_coeff 1 2 2 nu4
+pair_coeff 1 2 3 nu5
+pair_coeff 1 3 3 nu6
+pair_coeff 2 2 2 nu7
+pair_coeff 2 2 3 nu8
+pair_coeff 2 3 3 nu9
+pair_coeff 3 3 3 nu10 :pre
 
-[Steve:]
+By default the nu value for all triplets is set to 0.0.  Thus it is
-
+not required to provide pair_coeff commands that enumerate triplet
-I think a better syntax for the pair coeff command might be this:
+interactions for all K types.  If some I,J,K combination is not
-
+speficied, then there will be no 3-body ATM interactions for that
-pair_coeff I J v1 v2 ... vN
+combination and all its permutations.  However, as with all pair
-
+styles, it is required to specify a pair_coeff command for all I,J
-when 1,2,...N are the number of atom types defined.
+combinations, else an error will result.
-Then there be one pair_coeff command for each type pair,
-the same syntax as all other potentials in LAMMPS use.
-
-[Sergey:]
-
-The reason for my original pair_coeff command syntax is that
-I would like to reserve further arguments for possible extension of
-ATM potential in LAMMPS to further terms in the multipole expansion of
-many-body dispersion interaction.
-For example, the next term would be dipole-dipole-quadrupole, it may be
-activated if the next argument of pair_coeff is present
-without breaking backward compatibility.
-
-Also, the command you propose
-(i) will not account that the value of nu for different permutations
-of atom types is the same, and
-(ii) will make the numbering of atom types messy since there is
-no requirement to supply the values of nu for all triplets.
-
---------------------
-
-I think a better syntax for the pair coeff command might be this:
-
-pair_coeff I J v1 v2 ... vN
-
-when 1,2,...N are the number of atom types defined.
-Then there be one pair_coeff command for each type pair,
-the same syntax as all other potentials in LAMMPS use.
-
-Note that you refer to "elements", but the pair coeff command
-knows nothing about elements.  Only atom types.  There
-could be 10 atom types that all map to the same chemical
-element.
 
 :line
 

doc/src/pair_style.txt

@@ -103,6 +103,7 @@ in the pair section of "this page"_Section_commands.html#cmd_5.
 "pair_style adp"_pair_adp.html - angular dependent potential (ADP) of Mishin
 "pair_style airebo"_pair_airebo.html - AIREBO potential of Stuart
 "pair_style airebo/morse"_pair_airebo.html - AIREBO with Morse instead of LJ
+"pair_style atm"_pair_atm.html - Axilrod-Teller-Muto potential
 "pair_style beck"_pair_beck.html - Beck potential
 "pair_style body"_pair_body.html - interactions between body particles
 "pair_style bop"_pair_bop.html - BOP potential of Pettifor

examples/atm/log.9Jul18.atm.g++.1

@@ -26,7 +26,7 @@ Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252)
   1 by 1 by 1 MPI processor grid
 create_atoms	1 box
 Created 4000 atoms
-  Time spent = 0.00120211 secs
+  Time spent = 0.00125098 secs
 
 pair_style	hybrid/overlay lj/cut 4.5 atm 2.5
 pair_coeff	* * lj/cut 1.0 1.0
@@ -60,26 +60,26 @@ Neighbor list info ...
       bin: standard
 Per MPI rank memory allocation (min/avg/max) = 11.47 | 11.47 | 11.47 Mbytes
 Step Temp E_pair E_mol TotEng Press 
-       0        1.033   -4.5929584            0   -3.0438458   -3.6506231 
+       0        1.033    -4.733482            0   -3.1843694    -3.924644 
-       5    1.0335109   -4.5924034            0   -3.0425247   -3.6376817 
+       5    1.0335743   -4.7330528            0   -3.1830789   -3.9119065 
-      10    1.0347484   -4.5941952            0   -3.0424606   -3.6032204 
+      10    1.0349987   -4.7346788            0   -3.1825689   -3.8769962 
-      15    1.0357954   -4.5962409            0   -3.0429363   -3.5421887 
+      15    1.0362024   -4.7401425            0   -3.1862275   -3.8225106 
-      20    1.0350606    -4.595891            0   -3.0436883   -3.4478779 
+      20    1.0352365   -4.7459627            0   -3.1934962   -3.7403572 
-      25    1.0301813   -4.5896962            0   -3.0448107   -3.3111695 
+      25    1.0295963   -4.7444397            0   -3.2004313   -3.6132651 
-Loop time of 34.5602 on 1 procs for 25 steps with 4000 atoms
+Loop time of 27.841 on 1 procs for 25 steps with 4000 atoms
 
-Performance: 124.999 tau/day, 0.723 timesteps/s
+Performance: 155.167 tau/day, 0.898 timesteps/s
 100.0% CPU use with 1 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 34.556     | 34.556     | 34.556     |   0.0 | 99.99
+Pair    | 27.837     | 27.837     | 27.837     |   0.0 | 99.99
 Neigh   | 0          | 0          | 0          |   0.0 |  0.00
-Comm    | 0.001323   | 0.001323   | 0.001323   |   0.0 |  0.00
+Comm    | 0.0015321  | 0.0015321  | 0.0015321  |   0.0 |  0.01
-Output  | 0.00018287 | 0.00018287 | 0.00018287 |   0.0 |  0.00
+Output  | 0.00016594 | 0.00016594 | 0.00016594 |   0.0 |  0.00
-Modify  | 0.0016184  | 0.0016184  | 0.0016184  |   0.0 |  0.00
+Modify  | 0.001785   | 0.001785   | 0.001785   |   0.0 |  0.01
-Other   |            | 0.0006311  |            |       |  0.00
+Other   |            | 0.0006731  |            |       |  0.00
 
 Nlocal:    4000 ave 4000 max 4000 min
 Histogram: 1 0 0 0 0 0 0 0 0 0
@@ -97,4 +97,4 @@ Dangerous builds = 0
 
 Please see the log.cite file for references relevant to this simulation
 
-Total wall time: 0:00:36
+Total wall time: 0:00:29

examples/atm/log.9Jul18.atm.g++.4

@@ -26,7 +26,7 @@ Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252)
   1 by 2 by 2 MPI processor grid
 create_atoms	1 box
 Created 4000 atoms
-  Time spent = 0.000735998 secs
+  Time spent = 0.000769138 secs
 
 pair_style	hybrid/overlay lj/cut 4.5 atm 2.5
 pair_coeff	* * lj/cut 1.0 1.0
@@ -60,26 +60,26 @@ Neighbor list info ...
       bin: standard
 Per MPI rank memory allocation (min/avg/max) = 5.532 | 5.532 | 5.532 Mbytes
 Step Temp E_pair E_mol TotEng Press 
-       0        1.033   -4.5929584            0   -3.0438458   -3.6506231 
+       0        1.033    -4.733482            0   -3.1843694    -3.924644 
-       5    1.0335109   -4.5924034            0   -3.0425247   -3.6376817 
+       5    1.0335743   -4.7330528            0   -3.1830789   -3.9119065 
-      10    1.0347484   -4.5941952            0   -3.0424606   -3.6032204 
+      10    1.0349987   -4.7346788            0   -3.1825689   -3.8769962 
-      15    1.0357954   -4.5962409            0   -3.0429363   -3.5421887 
+      15    1.0362024   -4.7401425            0   -3.1862275   -3.8225106 
-      20    1.0350606    -4.595891            0   -3.0436883   -3.4478779 
+      20    1.0352365   -4.7459627            0   -3.1934962   -3.7403572 
-      25    1.0301813   -4.5896962            0   -3.0448107   -3.3111695 
+      25    1.0295963   -4.7444397            0   -3.2004313   -3.6132651 
-Loop time of 10.061 on 4 procs for 25 steps with 4000 atoms
+Loop time of 7.22029 on 4 procs for 25 steps with 4000 atoms
 
-Performance: 429.382 tau/day, 2.485 timesteps/s
+Performance: 598.314 tau/day, 3.462 timesteps/s
 100.0% CPU use with 4 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 9.8393     | 9.9439     | 10.008     |   2.0 | 98.84
+Pair    | 7.1346     | 7.1684     | 7.1863     |   0.8 | 99.28
 Neigh   | 0          | 0          | 0          |   0.0 |  0.00
-Comm    | 0.052154   | 0.11613    | 0.22077    |  18.7 |  1.15
+Comm    | 0.032945   | 0.0509     | 0.084664   |   9.1 |  0.70
-Output  | 7.2241e-05 | 8.2552e-05 | 0.00011158 |   0.0 |  0.00
+Output  | 0.00010133 | 0.00011051 | 0.00013804 |   0.0 |  0.00
-Modify  | 0.00053763 | 0.00055265 | 0.00056624 |   0.0 |  0.01
+Modify  | 0.00059557 | 0.00060683 | 0.00061274 |   0.0 |  0.01
-Other   |            | 0.0002971  |            |       |  0.00
+Other   |            | 0.000318   |            |       |  0.00
 
 Nlocal:    1000 ave 1000 max 1000 min
 Histogram: 4 0 0 0 0 0 0 0 0 0
@@ -97,4 +97,4 @@ Dangerous builds = 0
 
 Please see the log.cite file for references relevant to this simulation
 
-Total wall time: 0:00:10
+Total wall time: 0:00:07

src/MANYBODY/pair_atm.cpp

@@ -61,10 +61,12 @@ PairATM::PairATM(LAMMPS *lmp) : Pair(lmp)
 PairATM::~PairATM()
 {
   if (copymode) return;
+
   if (allocated) {
-    memory->destroy(nu);
     memory->destroy(setflag);
     memory->destroy(cutsq);
+
+    memory->destroy(nu);
   }
 }
 
@@ -89,6 +91,8 @@ void PairATM::compute(int eflag, int vflag)
   double **f = atom->f;
   int *type = atom->type;
 
+  double cutoffsq = cut_global*cut_global;
+
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -127,10 +131,11 @@ void PairATM::compute(int eflag, int vflag)
         rjk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2];
 
         r6 = rij2*rik2*rjk2;
-        if (r6 > cut_sixth) continue;
+        if (r6 > cutoffsq) continue;
 
         nu_local = nu[type[i]][type[j]][type[k]];
         if (nu_local == 0.0) continue;
+
         interaction_ddd(nu_local,
                         r6,rij2,rik2,rjk2,rij,rik,rjk,fj,fk,eflag,evdwl);
 
@@ -152,32 +157,50 @@ void PairATM::compute(int eflag, int vflag)
   if (vflag_fdotr) virial_fdotr_compute();
 }
 
+/* ---------------------------------------------------------------------- */
+
+void PairATM::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+  memory->create(nu,n+1,n+1,n+1,"pair:nu");
+
+  // initialize all nu values to 0.0
+
+  for (int i = 1; i <= n; i++)
+    for (int j = 1; j <= n; j++)
+      for (int k = 1; k <= n; k++)
+        nu[i][j][k] = 0.0;
+}
+
 /* ----------------------------------------------------------------------
-   reads the input script line with arguments you define
+   global settings
 ------------------------------------------------------------------------- */
 
 void PairATM::settings(int narg, char **arg)
 {
   if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+
   cut_global = force->numeric(FLERR,arg[0]);
 }
 
 /* ----------------------------------------------------------------------
-   set coefficients for one i,j,k type triplet
+   set coefficients for one I,J,K type triplet
 ------------------------------------------------------------------------- */
 
 void PairATM::coeff(int narg, char **arg)
 {
-  if (narg != 4)
+  if (narg != 4) error->all(FLERR,"Incorrect args for pair coefficients");
-    error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
-  int n = atom->ntypes;
-  for (int i = 0; i <= n; i++)
-    for (int j = 0; j <= n; j++)
-      for (int k = 0; k <= n; k++)
-        nu[i][j][k] = 0.0;
-
   int ilo,ihi,jlo,jhi,klo,khi;
   force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
   force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
@@ -185,16 +208,11 @@ void PairATM::coeff(int narg, char **arg)
 
   double nu_one = force->numeric(FLERR,arg[3]);
 
-  cut_sixth = cut_global*cut_global;
-  cut_sixth = cut_sixth*cut_sixth*cut_sixth;
-
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j<=jhi; j++) {
       for (int k = MAX(klo,j); k<=khi; k++) {
-        nu[i][j][k] = nu[i][k][j] = 
+        nu[i][j][k] = nu_one;
-        nu[j][i][k] = nu[j][k][i] = 
-        nu[k][i][j] = nu[k][j][i] = nu_one;
         count++;
       }
       setflag[i][j] = 1;
@@ -211,18 +229,28 @@ void PairATM::coeff(int narg, char **arg)
 void PairATM::init_style()
 {
   // need a full neighbor list
+
   int irequest = neighbor->request(this,instance_me);
   neighbor->requests[irequest]->half = 0;
   neighbor->requests[irequest]->full = 1;
 }
 
 /* ----------------------------------------------------------------------
-   perform initialization for one i,j type pair
+   init for one i,j type pair and corresponding j,i
+   also for all k type permutations
 ------------------------------------------------------------------------- */
 
 double PairATM::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+  // set all 6 symmetric permutations of I,J,K types to same nu value
+
+  int ntypes = atom->ntypes;
+  for (int k = j; k <= ntypes; k++)
+    nu[i][k][j] = nu[j][i][k] = nu[j][k][i] = nu[k][i][j] = nu[k][j][i] = 
+      nu[i][j][k];
+
   return cut_global;
 }
 
@@ -239,7 +267,7 @@ void PairATM::write_restart(FILE *fp)
     for (j = i; j <= atom->ntypes; j++) {
       fwrite(&setflag[i][j],sizeof(int),1,fp);
       if (setflag[i][j]) 
-        for (k = i; k <= atom->ntypes; k++) 
+        for (k = j; k <= atom->ntypes; k++) 
           fwrite(&nu[i][j][k],sizeof(double),1,fp);
     }
   }
@@ -260,7 +288,7 @@ void PairATM::read_restart(FILE *fp)
     for (j = i; j <= atom->ntypes; j++) {
       if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
-      if (setflag[i][j]) for (k = i; k <= atom->ntypes; k++) {
+      if (setflag[i][j]) for (k = j; k <= atom->ntypes; k++) {
         if (me == 0) fread(&nu[i][j][k],sizeof(double),1,fp);
         MPI_Bcast(&nu[i][j][k],1,MPI_DOUBLE,0,world);
       }
@@ -288,17 +316,6 @@ void PairATM::read_restart_settings(FILE *fp)
   MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
 }
 
-/* ---------------------------------------------------------------------- */
-
-void PairATM::allocate()
-{
-  allocated = 1;
-  int n = atom->ntypes;
-  memory->create(nu,n+1,n+1,n+1,"pair:a");
-  memory->create(setflag,n+1,n+1,"pair:setflag");
-  memory->create(cutsq,n+1,n+1,"pair:cutsq");
-}
-
 /* ----------------------------------------------------------------------
    Axilrod-Teller-Muto (dipole-dipole-dipole) potential
 ------------------------------------------------------------------------- */
@@ -314,14 +331,14 @@ void PairATM::interaction_ddd(double nu,
   rrj = rij[0]*rjk[0] + rij[1]*rjk[1] + rij[2]*rjk[2];
   rrk = rjk[0]*rik[0] + rjk[1]*rik[1] + rjk[2]*rik[2];
   rrr = 5.0*rri*rrj*rrk;
-  for (int i=0; i<3; i++) {
+  for (int i = 0; i < 3; i++) {
-    fj[i] = rrj*(rrk - rri)*rik[i]
+    fj[i] = rrj*(rrk - rri)*rik[i] - 
-      - (rrk*rri - rjk2*rik2 + rrr/rij2)*rij[i]
+      (rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] + 
-      + (rrk*rri - rik2*rij2 + rrr/rjk2)*rjk[i];
+      (rrk*rri - rik2*rij2 + rrr/rjk2) * rjk[i];
     fj[i] *= 3.0*r5inv;
-    fk[i] = rrk*(rri + rrj)*rij[i]
+    fk[i] = rrk*(rri + rrj)*rij[i] + 
-      + (rri*rrj + rik2*rij2 - rrr/rjk2)*rjk[i]
+      (rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] + 
-      + (rri*rrj + rij2*rjk2 - rrr/rik2)*rik[i];
+      (rri*rrj + rij2*rjk2 - rrr/rik2) * rik[i];
     fk[i] *= 3.0*r5inv;
   }
   if (eflag) eng = (r6 - 0.6*rrr)*r5inv;