From d4385ade159f0a2d163ace5fb38328e2d7d8da8e Mon Sep 17 00:00:00 2001
From: "Steven J. Plimpton" <sjplimp@singsing.sandia.gov>
Date: Mon, 16 Jul 2018 11:22:22 -0600
Subject: [PATCH] more changes to ATM source and doc file
---
doc/src/Section_commands.txt | 3 +-
doc/src/pair_atm.txt | 132 +++++++++++++-----------------
doc/src/pair_style.txt | 1 +
examples/atm/log.9Jul18.atm.g++.1 | 30 +++----
examples/atm/log.9Jul18.atm.g++.4 | 30 +++----
src/MANYBODY/pair_atm.cpp | 101 +++++++++++++----------
6 files changed, 151 insertions(+), 146 deletions(-)
diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index e44595455f..49ade2ca1b 100644
--- a/doc/src/Section_commands.txt
+++ b/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,
diff --git a/doc/src/pair_atm.txt b/doc/src/pair_atm.txt
index 96e002f2f8..e442c32814 100644
--- a/doc/src/pair_atm.txt
+++ b/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,
-and the distances between pairs of atoms r12, r23 and r31
-and the angles gamma1, gamma2 and gamma3
-are shown at the diagram:
+where nu is the three-body interaction strength. The distances
+between pairs of atoms r12, r23, r31 and the angles gamma1, gamma2,
+gamma3 are as shown in this diagram:
:c,image(JPG/pair_atm_dia.jpg)
-There is no \"central\" atom, the interaction is symmetric with respect
-to permutation of atom types.
+Note that for the interaction between a triplet of atoms I,J,K, there
+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
-"neighbor list"_neighbor.html
-and the distances between atoms satisfy r12 r23 r31 > cutoff^3.
+The potential for a triplet of atom is calculated only if all 3
+distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff.
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
-nu (energy/distance^9 units) :ul
+K = atom type of the third atom (1 to Ntypes)
+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.
-For a two-atom type simulation, more pair_coeff commands can be used.
-ATM interaction is symmetric with respect to permutation of atoms,
-it is necessary to provide pair_coeff command for one permutation.
-For example, the command
+Note that a pair_coeff command can override a previous setting for the
+same I,J,K triplet. For example, these commands set nu for all I,J.K
+triplets, then overwrite nu for just the I,J,K = 2,3,4 triplet:
-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 2 1 1 nu :pre
+pair_coeff 1 1 1 0.25 :per
-Thus, to specify all ATM interactions between two atom types (eg 1 and 2)
-it is sufficient to provide four pair_coeff commands, eg:
+For a simulation with two atom types, four pair_coeff commands will
+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
-(eg 111, 112, 113, 122, 123, 133, 222, 223, 233, 333)
-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.
+For a simulation with three atom types, ten pair_coeff commands will
+specify all possible nu values:
---------------------
-
-[Steve:]
-
-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.
-
-[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.
+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
+
+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
+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
+combination and all its permutations. However, as with all pair
+styles, it is required to specify a pair_coeff command for all I,J
+combinations, else an error will result.
:line
diff --git a/doc/src/pair_style.txt b/doc/src/pair_style.txt
index 475761add7..7e29544a70 100644
--- a/doc/src/pair_style.txt
+++ b/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
diff --git a/examples/atm/log.9Jul18.atm.g++.1 b/examples/atm/log.9Jul18.atm.g++.1
index 27ba23e3a9..0613533f33 100644
--- a/examples/atm/log.9Jul18.atm.g++.1
+++ b/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
- 5 1.0335109 -4.5924034 0 -3.0425247 -3.6376817
- 10 1.0347484 -4.5941952 0 -3.0424606 -3.6032204
- 15 1.0357954 -4.5962409 0 -3.0429363 -3.5421887
- 20 1.0350606 -4.595891 0 -3.0436883 -3.4478779
- 25 1.0301813 -4.5896962 0 -3.0448107 -3.3111695
-Loop time of 34.5602 on 1 procs for 25 steps with 4000 atoms
+ 0 1.033 -4.733482 0 -3.1843694 -3.924644
+ 5 1.0335743 -4.7330528 0 -3.1830789 -3.9119065
+ 10 1.0349987 -4.7346788 0 -3.1825689 -3.8769962
+ 15 1.0362024 -4.7401425 0 -3.1862275 -3.8225106
+ 20 1.0352365 -4.7459627 0 -3.1934962 -3.7403572
+ 25 1.0295963 -4.7444397 0 -3.2004313 -3.6132651
+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
-Output | 0.00018287 | 0.00018287 | 0.00018287 | 0.0 | 0.00
-Modify | 0.0016184 | 0.0016184 | 0.0016184 | 0.0 | 0.00
-Other | | 0.0006311 | | | 0.00
+Comm | 0.0015321 | 0.0015321 | 0.0015321 | 0.0 | 0.01
+Output | 0.00016594 | 0.00016594 | 0.00016594 | 0.0 | 0.00
+Modify | 0.001785 | 0.001785 | 0.001785 | 0.0 | 0.01
+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
diff --git a/examples/atm/log.9Jul18.atm.g++.4 b/examples/atm/log.9Jul18.atm.g++.4
index 300f82d581..5627348fab 100644
--- a/examples/atm/log.9Jul18.atm.g++.4
+++ b/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
- 5 1.0335109 -4.5924034 0 -3.0425247 -3.6376817
- 10 1.0347484 -4.5941952 0 -3.0424606 -3.6032204
- 15 1.0357954 -4.5962409 0 -3.0429363 -3.5421887
- 20 1.0350606 -4.595891 0 -3.0436883 -3.4478779
- 25 1.0301813 -4.5896962 0 -3.0448107 -3.3111695
-Loop time of 10.061 on 4 procs for 25 steps with 4000 atoms
+ 0 1.033 -4.733482 0 -3.1843694 -3.924644
+ 5 1.0335743 -4.7330528 0 -3.1830789 -3.9119065
+ 10 1.0349987 -4.7346788 0 -3.1825689 -3.8769962
+ 15 1.0362024 -4.7401425 0 -3.1862275 -3.8225106
+ 20 1.0352365 -4.7459627 0 -3.1934962 -3.7403572
+ 25 1.0295963 -4.7444397 0 -3.2004313 -3.6132651
+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
-Output | 7.2241e-05 | 8.2552e-05 | 0.00011158 | 0.0 | 0.00
-Modify | 0.00053763 | 0.00055265 | 0.00056624 | 0.0 | 0.01
-Other | | 0.0002971 | | | 0.00
+Comm | 0.032945 | 0.0509 | 0.084664 | 9.1 | 0.70
+Output | 0.00010133 | 0.00011051 | 0.00013804 | 0.0 | 0.00
+Modify | 0.00059557 | 0.00060683 | 0.00061274 | 0.0 | 0.01
+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
diff --git a/src/MANYBODY/pair_atm.cpp b/src/MANYBODY/pair_atm.cpp
index 0db076acd0..2ba8912899 100644
--- a/src/MANYBODY/pair_atm.cpp
+++ b/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)
- error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg != 4) 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[j][i][k] = nu[j][k][i] =
- nu[k][i][j] = nu[k][j][i] = nu_one;
+ nu[i][j][k] = 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,25 +316,14 @@ 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
------------------------------------------------------------------------- */
void PairATM::interaction_ddd(double nu,
- double r6, double rij2, double rik2, double rjk2,
- double *rij, double *rik, double *rjk,
- double *fj, double *fk, int eflag, double &eng)
+ double r6, double rij2, double rik2, double rjk2,
+ double *rij, double *rik, double *rjk,
+ double *fj, double *fk, int eflag, double &eng)
{
double r5inv,rri,rrj,rrk,rrr;
r5inv = nu / (r6*r6*sqrt(r6));
@@ -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++) {
- fj[i] = rrj*(rrk - rri)*rik[i]
- - (rrk*rri - rjk2*rik2 + rrr/rij2)*rij[i]
- + (rrk*rri - rik2*rij2 + rrr/rjk2)*rjk[i];
+ for (int i = 0; i < 3; i++) {
+ fj[i] = rrj*(rrk - rri)*rik[i] -
+ (rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] +
+ (rrk*rri - rik2*rij2 + rrr/rjk2) * rjk[i];
fj[i] *= 3.0*r5inv;
- fk[i] = rrk*(rri + rrj)*rij[i]
- + (rri*rrj + rik2*rij2 - rrr/rjk2)*rjk[i]
- + (rri*rrj + rij2*rjk2 - rrr/rik2)*rik[i];
+ fk[i] = rrk*(rri + rrj)*rij[i] +
+ (rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] +
+ (rri*rrj + rij2*rjk2 - rrr/rik2) * rik[i];
fk[i] *= 3.0*r5inv;
}
if (eflag) eng = (r6 - 0.6*rrr)*r5inv;
--
GitLab