Merge pull request #2556 from akohlmey/eam-he-pair-style

Add eam/he pair style for modeling He bubbles in metals

doc/src/Commands_pair.rst

@@ -96,6 +96,7 @@ OPT.
    * :doc:`eam/cd <pair_eam>`
    * :doc:`eam/cd/old <pair_eam>`
    * :doc:`eam/fs (gikot) <pair_eam>`
+   * :doc:`eam/he <pair_eam>`
    * :doc:`edip (o) <pair_edip>`
    * :doc:`edip/multi <pair_edip>`
    * :doc:`edpd <pair_mesodpd>`

doc/src/pair_eam.rst

@@ -18,6 +18,7 @@
 .. index:: pair_style eam/fs/kk
 .. index:: pair_style eam/fs/omp
 .. index:: pair_style eam/fs/opt
+.. index:: pair_style eam/he
 
 pair_style eam command
 ======================
@@ -38,6 +39,9 @@ pair_style eam/cd/old command
 pair_style eam/fs command
 =========================
 
+pair_style eam/he command
+=========================
+
 Accelerator Variants: *eam/fs/gpu*, *eam/fs/intel*, *eam/fs/kk*, *eam/fs/omp*, *eam/fs/opt*
 
 Syntax
@@ -47,7 +51,7 @@ Syntax
 
    pair_style style
 
-* style = *eam* or *eam/alloy* or *eam/cd* or *eam/cd/old* or *eam/fs*
+* style = *eam* or *eam/alloy* or *eam/cd* or *eam/cd/old* or *eam/fs* or *eam/he*
 
 Examples
 """"""""
@@ -67,6 +71,9 @@ Examples
    pair_style eam/fs
    pair_coeff * * NiAlH_jea.eam.fs Ni Al Ni Ni
 
+   pair_style eam/he
+   pair_coeff * * PdHHe.eam.he Pd H He
+
 Description
 """""""""""
 
@@ -104,8 +111,8 @@ are parameterized in terms of LAMMPS :doc:`metal units <units>`.
    potentials, the same way that DYNAMO does.  Alternatively, a single
    DYNAMO *setfl* file or Finnis/Sinclair EAM file can be used by LAMMPS
    to model alloy systems by invoking the *eam/alloy* or *eam/cd* or
-   *eam/fs* styles as described below.  These files require no mixing
+   *eam/fs* or *eam/he* styles as described below.  These files require no
-   since they specify alloy interactions explicitly.
+   mixing since they specify alloy interactions explicitly.
 
 .. note::
 
@@ -143,10 +150,6 @@ DYNAMO single-element *funcfl* format.  If the DYNAMO file was created
 by a Fortran program, it cannot have "D" values in it for exponents.
 C only recognizes "e" or "E" for scientific notation.
 
-Note that unlike for other potentials, cutoffs for EAM potentials are
-not set in the pair_style or pair_coeff command; they are specified in
-the EAM potential files themselves.
-
 For style *eam* a potential file must be assigned to each I,I pair of
 atom types by using one or more pair_coeff commands, each with a
 single argument:
@@ -336,8 +339,11 @@ distribution have a ".cdeam" suffix.
 
 Style *eam/fs* computes pairwise interactions for metals and metal
 alloys using a generalized form of EAM potentials due to Finnis and
-Sinclair :ref:`(Finnis) <Finnis1>`.  The total energy Ei of an atom I is
+Sinclair :ref:`(Finnis) <Finnis1>`.  Style *eam/he* is similar to
-given by
+*eam/fs* except that it allows for negative electron density in
+order to capture the behavior of helium in metals :ref:`(Zhou6) <Zhou6>`.
+
+The total energy Ei of an atom I is given by
 
 .. math::
 
@@ -355,36 +361,36 @@ electron density at an atomic site depending on the identity of the
 element at that atomic site.
 
 The associated :doc:`pair_coeff <pair_coeff>` command for style *eam/fs*
-reads a DYNAMO *setfl* file that has been extended to include
+or *eam/he* reads a DYNAMO *setfl* file that has been extended to include
-additional rho_alpha_beta arrays of tabulated values.  A discussion of
+additional :math:`\rho_{\alpha\beta}` arrays of tabulated values.  A
-how FS EAM differs from conventional EAM alloy potentials is given in
+discussion of how FS EAM differs from conventional EAM alloy potentials is
-:ref:`(Ackland1) <Ackland1>`.  An example of such a potential is the same
+given in :ref:`(Ackland1) <Ackland1>`.  An example of such a potential is the
-author's Fe-P FS potential :ref:`(Ackland2) <Ackland2>`.  Note that while FS
+same author's Fe-P FS potential :ref:`(Ackland2) <Ackland2>`.  Note that while
-potentials always specify the embedding energy with a square root
+FS potentials always specify the embedding energy with a square root
 dependence on the total density, the implementation in LAMMPS does not
 require that; the user can tabulate any functional form desired in the
 FS potential files.
 
-For style *eam/fs*\ , the form of the pair_coeff command is exactly the
+For style *eam/fs* and *eam/he* the form of the pair_coeff command is exactly
-same as for style *eam/alloy*\ , e.g.
+the same as for style *eam/alloy*\ , e.g.
 
 .. code-block:: LAMMPS
 
    pair_coeff * * NiAlH_jea.eam.fs Ni Ni Ni Al
 
-where there are N additional arguments after the filename, where N is
+with N additional arguments after the filename, where N is
 the number of LAMMPS atom types.  See the :doc:`pair_coeff <pair_coeff>`
 doc page for alternate ways to specify the path for the potential
 file.  The N values determine the mapping of LAMMPS atom types to EAM
 elements in the file, as described above for style *eam/alloy*\ .  As
 with *eam/alloy*\ , if a mapping value is NULL, the mapping is not
-performed.  This can be used when an *eam/fs* potential is used as
+performed.  This can be used when an *eam/fs* or *eam/he* potential is
-part of the *hybrid* pair style.  The NULL values are used as
+used as part of a *hybrid* pair style.  The NULL values are used as
 placeholders for atom types that will be used with other potentials.
 
-FS EAM files include more information than the DYNAMO *setfl* format
+FS EAM and HE EAM files include more information than the DYNAMO *setfl*
-files read by *eam/alloy*\ , in that i,j density functionals for all
+format files read by *eam/alloy*\ , in that i,j density functionals for
-pairs of elements are included as needed by the Finnis/Sinclair
+all pairs of elements are included as needed by the Finnis/Sinclair
 formulation of the EAM.
 
 FS EAM files in the *potentials* directory of the LAMMPS distribution
@@ -417,6 +423,25 @@ eV-Angstroms) as in EAM *setfl* files.  Note that in Finnis/Sinclair,
 the phi(r) arrays are still symmetric, so only phi arrays for i >= j
 are listed.
 
+HE EAM files in the *potentials* directory of the LAMMPS distribution
+have an ".eam.he" suffix.  They are formatted as follows:
+
+* lines 1,2,3 = comments (ignored)
+* line 4: Nelements Element1 Element2 ... ElementN
+* line 5: Nrho, drho, Nr, dr, cutoff, rhomax
+
+The 5-line header section is identical to an FS EAM file
+except that line 5 lists an additional value, rhomax. Unlike in FS EAM
+files where embedding energies F(rho) are always defined between rho = 0
+and rho = (Nrho -1)drho, F(rho) in HE EAM files are defined between
+rho = rhomin and rho = rhomax.  Since drho = (rhomax - rhomin)/(Nrho - 1),
+rhomin = rhomax - (Nrho - 1)drho.  The embedding energies F(rho) are
+listed for rho = rhomin, rhomin + drho, rhomin + 2drho, ..., rhomax.
+This gives users additional flexibility to define a negative rhomin and
+therefore an embedding energy function that works for both positive and
+negative electron densities.  The format and units of these sections are
+identical to the FS EAM files (see above).
+
 ----------
 
 .. include:: accel_styles.rst
@@ -480,6 +505,10 @@ Daw, Baskes, Phys Rev B, 29, 6443 (1984).
 
 **(Finnis)** Finnis, Sinclair, Philosophical Magazine A, 50, 45 (1984).
 
+.. _Zhou6:
+
+**(Zhou6)** Zhou, Bartelt, Sills, Physical Review B, 103, 014108 (2021).
+
 .. _Stukowski:
 
 **(Stukowski)** Stukowski, Sadigh, Erhart, Caro; Modeling Simulation

doc/src/pair_style.rst

@@ -159,6 +159,7 @@ accelerated styles exist.
 * :doc:`eam/cd <pair_eam>` - concentration-dependent EAM
 * :doc:`eam/cd/old <pair_eam>` - older two-site model for concentration-dependent EAM
 * :doc:`eam/fs <pair_eam>` - Finnis-Sinclair EAM
+* :doc:`eam/he <pair_eam>` - Finnis-Sinclair EAM modified for Helium in metals
 * :doc:`edip <pair_edip>` - three-body EDIP potential
 * :doc:`edip/multi <pair_edip>` - multi-element EDIP potential
 * :doc:`edpd <pair_mesodpd>` - eDPD particle interactions

doc/utils/sphinx-config/false_positives.txt

@@ -200,6 +200,7 @@ barostatted
 barostatting
 Barostatting
 Barrat
+Bartelt
 Bartels
 barycenter
 barye
@@ -2678,6 +2679,8 @@ rhodo
 Rhodo
 rhodopsin
 rhok
+rhomax
+rhomin
 rhorho
 rhosum
 ri

src/KSPACE/fft3d.cpp

@@ -734,7 +734,7 @@ void fft_1d_only(FFT_DATA *data, int nsize, int flag, struct fft_plan_3d *plan)
 
 // fftw3 and Dfti in MKL encode the number of transforms
 // into the plan, so we cannot operate on a smaller data set
-  
+
 #if defined(FFT_MKL) || defined(FFT_FFTW3)
   if ((total1 > nsize) || (total2 > nsize) || (total3 > nsize))
     return;

src/KSPACE/gridcomm.cpp

@@ -525,7 +525,7 @@ void GridComm::setup_regular(int &nbuf1, int &nbuf2)
     ngrid = MAX(ngrid,swap[i].npack);
     ngrid = MAX(ngrid,swap[i].nunpack);
   }
-  
+
   nbuf1 = nbuf2 = ngrid;
 }
 

src/KSPACE/pppm_disp.cpp

@@ -1289,13 +1289,13 @@ void PPPMDisp::init_coeffs()
   int tmp;
   int n = atom->ntypes;
   int converged;
-  
+
   delete [] B;
   B = nullptr;
 
   // no mixing rule or arithmetic
-  
+
-  if (function[3] + function[2]) {                     
+  if (function[3] + function[2]) {
     if (function[2] && me == 0)
       utils::logmesg(lmp,"  Optimizing splitting of Dispersion coefficients\n");
 
@@ -1323,11 +1323,11 @@ void PPPMDisp::init_coeffs()
         error->all(FLERR,
 		   "Matrix factorization to split dispersion coefficients failed");
       }
-      
+
       // determine number of used eigenvalues
       //   based on maximum allowed number or cutoff criterion
       //   sort eigenvalues according to their size with bubble sort
-      
+
       double t;
       for (int i = 0; i < n; i++) {
         for (int j = 0; j < n-1-i; j++) {
@@ -1346,7 +1346,7 @@ void PPPMDisp::init_coeffs()
 
       // check which eigenvalue is the first that is smaller than a specified tolerance
       // check how many are maximum allowed by the user
-      
+
       double amax = fabs(A[0][0]);
       double acrit = amax*splittol;
       double bmax = 0;
@@ -1365,7 +1365,7 @@ void PPPMDisp::init_coeffs()
         error->warning(FLERR,fmt::format("Estimated error in splitting of "
                                          "dispersion coeffs is {}",err));
       // set B
-      
+
       B = new double[nsplit*n+nsplit];
       for (int i = 0; i < nsplit; i++) {
         B[i] = A[i][i];
@@ -1376,11 +1376,11 @@ void PPPMDisp::init_coeffs()
 
       nsplit_alloc = nsplit;
       if (nsplit % 2 == 1) nsplit_alloc = nsplit + 1;
-      
+
     } else nsplit = 1;      // use geometric mixing
 
     // check if the function should preferably be [1] or [2] or [3]
-    
+
     if (nsplit == 1) {
       if (B) delete [] B;
       function[3] = 0;
@@ -1389,21 +1389,21 @@ void PPPMDisp::init_coeffs()
       if (me == 0)
         utils::logmesg(lmp,"  Using geometric mixing for reciprocal space\n");
     }
-    
+
     if (function[2] && nsplit <= 6) {
       if (me == 0)
         utils::logmesg(lmp,fmt::format("  Using {} instead of 7 structure "
                                        "factors\n",nsplit));
       //function[3] = 1;
-      //function[2] = 0; 
+      //function[2] = 0;
       if (B) delete [] B;   // remove this when un-comment previous 2 lines
    }
-    
+
     if (function[2] && (nsplit > 6)) {
       if (me == 0) utils::logmesg(lmp,"  Using 7 structure factors\n");
       if (B) delete [] B;
     }
-    
+
     if (function[3]) {
       if (me == 0)
         utils::logmesg(lmp,fmt::format("  Using {} structure factors\n",
@@ -1416,14 +1416,14 @@ void PPPMDisp::init_coeffs()
     memory->destroy(A);
     memory->destroy(Q);
   }
-  
+
   if (function[1]) {                                    // geometric 1/r^6
     double **b = (double **) force->pair->extract("B",tmp);
     B = new double[n+1];
     B[0] = 0.0;
     for (int i=1; i<=n; ++i) B[i] = sqrt(fabs(b[i][i]));
   }
-  
+
   if (function[2]) {                                    // arithmetic 1/r^6
     double **epsilon = (double **) force->pair->extract("epsilon",tmp);
     double **sigma = (double **) force->pair->extract("sigma",tmp);
@@ -1432,7 +1432,7 @@ void PPPMDisp::init_coeffs()
     double eps_i,sigma_i,sigma_n;
     B = new double[7*n+7];
     double c[7] = {1.0,sqrt(6.0),sqrt(15.0),sqrt(20.0),sqrt(15.0),sqrt(6.0),1.0};
-    
+
     for (int i=1; i<=n; ++i) {
       eps_i = sqrt(epsilon[i][i]);
       sigma_i = sigma[i][i];
@@ -1567,22 +1567,22 @@ void PPPMDisp::qr_tri(double** Qi, double** A, int n)
       Qi[i][j] = 0.0;
   for (int i = 0; i < n; i++)
     Qi[i][i] = 1.0;
-  
+
   // loop over main diagonal and first of diagonal of A
 
   for (int i = 0; i < n-1; i++) {
     j = i+1;
-    
+
     // coefficients of the rotation matrix
-    
+
     a = A[i][i];
     b = A[j][i];
     r = sqrt(a*a + b*b);
     c = a/r;
     s = b/r;
-    
+
     // update the entries of A and Q
-    
+
     k0 = (i-1>0)?i-1:0;   //min(i-1,0);
     kmax = (i+3<n)?i+3:n;  //min(i+3,n);
     for (k = k0; k < kmax; k++) {
@@ -1612,14 +1612,14 @@ void PPPMDisp::mmult(double** A, double** B, double** C, int n)
       C[i][j] = 0.0;
 
   // perform matrix multiplication
-  
+
   for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++)
       for (int k = 0; k < n; k++)
         C[i][j] += A[i][k] * B[k][j];
 
   // copy the result back to matrix A
-  
+
   for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++)
       A[i][j] = C[i][j];
@@ -1638,9 +1638,9 @@ int PPPMDisp::check_convergence(double** A, double** Q, double** A0,
   double epsmax = -1;
   double Bmax = 0.0;
   double diff;
-  
+
   // get the largest eigenvalue of the original matrix
-  
+
   for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++)
       Bmax = (Bmax>A0[i][j])?Bmax:A0[i][j];  //max(Bmax,A0[i][j]);
@@ -1648,36 +1648,36 @@ int PPPMDisp::check_convergence(double** A, double** Q, double** A0,
 
   // reconstruct the original matrix
   // store the diagonal elements in D
-  
+
   for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++)
       D[i][j] = 0.0;
   for (int i = 0; i < n; i++)
     D[i][i] = A[i][i];
-  
+
   // store matrix Q in E
-  
+
   for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++)
       E[i][j] = Q[i][j];
-  
+
   // E = Q*A
-  
+
   mmult(E,D,C,n);
-  
+
   // store transpose of Q in D
-  
+
   for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++)
       D[i][j] = Q[j][i];
 
   // E = Q*A*Q.t
-  
+
   mmult(E,D,C,n);
 
 
   //compare the original matrix and the final matrix
-  
+
   for (int i = 0; i < n; i++) {
     for (int j = 0; j < n; j++) {
       diff = A0[i][j] - E[i][j];
@@ -2690,15 +2690,15 @@ void PPPMDisp::set_grid()
    set the FFT parameters
 ------------------------------------------------------------------------- */
 
-void PPPMDisp::set_fft_parameters(int& nx_p, int& ny_p, int& nz_p, 
+void PPPMDisp::set_fft_parameters(int& nx_p, int& ny_p, int& nz_p,
-                                  int& nxlo_f, int& nylo_f, int& nzlo_f, 
+                                  int& nxlo_f, int& nylo_f, int& nzlo_f,
-                                  int& nxhi_f, int& nyhi_f, int& nzhi_f, 
+                                  int& nxhi_f, int& nyhi_f, int& nzhi_f,
-                                  int& nxlo_i, int& nylo_i, int& nzlo_i, 
+                                  int& nxlo_i, int& nylo_i, int& nzlo_i,
-                                  int& nxhi_i, int& nyhi_i, int& nzhi_i, 
+                                  int& nxhi_i, int& nyhi_i, int& nzhi_i,
-                                  int& nxlo_o, int& nylo_o, int& nzlo_o, 
+                                  int& nxlo_o, int& nylo_o, int& nzlo_o,
-                                  int& nxhi_o, int& nyhi_o, int& nzhi_o, 
+                                  int& nxhi_o, int& nyhi_o, int& nzhi_o,
-                                  int& nlow,  int& nupp, 
+                                  int& nlow,  int& nupp,
-                                  int& ng, int& nf, int& nfb, 
+                                  int& ng, int& nf, int& nfb,
                                   double& sft, double& sftone, int& ord)
 {
   // global indices of PPPM grid range from 0 to N-1
@@ -2924,8 +2924,8 @@ double PPPMDisp::derivf()
   double df,f1,f2,g_ewald_old;
 
   // derivative step-size
-  
+
-  double h = 0.000001;  
+  double h = 0.000001;
 
   f1 = f();
   g_ewald_old = g_ewald;
@@ -3003,7 +3003,7 @@ double PPPMDisp::compute_qopt_6()
   double qopt;
   if (differentiation_flag == 1) qopt = compute_qopt_6_ad();
   else qopt = compute_qopt_6_ik();
-  
+
   double qopt_all;
   MPI_Allreduce(&qopt,&qopt_all,1,MPI_DOUBLE,MPI_SUM,world);
   return qopt_all;
@@ -4414,7 +4414,7 @@ void PPPMDisp::make_rho_g()
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
   // (dx,dy,dz) = distance to "lower left" grid pt
   // (mx,my,mz) = global coords of moving stencil pt
-  
+
   int type;
   double **x = atom->x;
   int nlocal = atom->nlocal;
@@ -4428,7 +4428,7 @@ void PPPMDisp::make_rho_g()
     dz = nz+shiftone_6 - (x[i][2]-boxlo[2])*delzinv_6;
 
     compute_rho1d(dx,dy,dz,order_6,rho_coeff_6,rho1d_6);
-    
+
     type = atom->type[i];
     z0 = delvolinv_6 * B[type];
     for (n = nlower_6; n <= nupper_6; n++) {
@@ -4479,7 +4479,7 @@ void PPPMDisp::make_rho_a()
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
   // (dx,dy,dz) = distance to "lower left" grid pt
   // (mx,my,mz) = global coords of moving stencil pt
-  
+
   int type;
   double **x = atom->x;
   int nlocal = atom->nlocal;
@@ -4491,9 +4491,9 @@ void PPPMDisp::make_rho_a()
     dx = nx+shiftone_6 - (x[i][0]-boxlo[0])*delxinv_6;
     dy = ny+shiftone_6 - (x[i][1]-boxlo[1])*delyinv_6;
     dz = nz+shiftone_6 - (x[i][2]-boxlo[2])*delzinv_6;
-    
+
     compute_rho1d(dx,dy,dz,order_6,rho_coeff_6,rho1d_6);
-    
+
     type = atom->type[i];
     z0 = delvolinv_6;
     for (n = nlower_6; n <= nupper_6; n++) {
@@ -4531,7 +4531,7 @@ void PPPMDisp::make_rho_none()
   FFT_SCALAR dx,dy,dz,x0,y0,z0,w;
 
   // clear 3d density array
-  
+
   for (k = 0; k < nsplit_alloc; k++)
     memset(&(density_brick_none[k][nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
            ngrid_6*sizeof(FFT_SCALAR));
@@ -4540,7 +4540,7 @@ void PPPMDisp::make_rho_none()
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
   // (dx,dy,dz) = distance to "lower left" grid pt
   // (mx,my,mz) = global coords of moving stencil pt
-  
+
   int type;
   double **x = atom->x;
   int nlocal = atom->nlocal;
@@ -4552,11 +4552,11 @@ void PPPMDisp::make_rho_none()
     dx = nx+shiftone_6 - (x[i][0]-boxlo[0])*delxinv_6;
     dy = ny+shiftone_6 - (x[i][1]-boxlo[1])*delyinv_6;
     dz = nz+shiftone_6 - (x[i][2]-boxlo[2])*delzinv_6;
-    
+
     compute_rho1d(dx,dy,dz,order_6,rho_coeff_6,rho1d_6);
-    
+
     type = atom->type[i];
-    z0 = delvolinv_6; 
+    z0 = delvolinv_6;
     for (n = nlower_6; n <= nupper_6; n++) {
       mz = n+nz;
       y0 = z0*rho1d_6[2][n];
@@ -4602,7 +4602,7 @@ void PPPMDisp::poisson_ik(FFT_SCALAR* wk1, FFT_SCALAR* wk2,
   double eng;
 
   // transform charge/dispersion density (r -> k)
-  
+
   n = 0;
   for (i = 0; i < nft; i++) {
     wk1[n++] = dfft[i];
@@ -4682,7 +4682,7 @@ void PPPMDisp::poisson_ik(FFT_SCALAR* wk1, FFT_SCALAR* wk2,
         }
 
     ft2->compute(wk2,wk2,FFT3d::BACKWARD);
-    
+
     n = 0;
     for (k = nzlo_i; k <= nzhi_i; k++)
       for (j = nylo_i; j <= nyhi_i; j++)
@@ -4705,7 +4705,7 @@ void PPPMDisp::poisson_ik(FFT_SCALAR* wk1, FFT_SCALAR* wk2,
         }
 
     ft2->compute(wk2,wk2,FFT3d::BACKWARD);
-    
+
     n = 0;
     for (k = nzlo_i; k <= nzhi_i; k++)
       for (j = nylo_i; j <= nyhi_i; j++)
@@ -4742,7 +4742,7 @@ void PPPMDisp::poisson_ad(FFT_SCALAR* wk1, FFT_SCALAR* wk2,
   double eng;
 
   // transform charge/dispersion density (r -> k)
-  
+
   n = 0;
   for (i = 0; i < nft; i++) {
     wk1[n++] = dfft[i];
@@ -4822,7 +4822,7 @@ void PPPMDisp::poisson_peratom(FFT_SCALAR* wk1, FFT_SCALAR* wk2, LAMMPS_NS::FFT3
 			       FFT_SCALAR*** v5_pa)
 {
  // v0 & v1 term
-  
+
   int n, i, j, k;
   n = 0;
   for (i = 0; i < nft; i++) {
@@ -4905,7 +4905,7 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
 
   // transform charge/dispersion density (r -> k)
   // only one transform when energies and pressures not calculated
-  
+
   if (eflag_global + vflag_global == 0) {
     n = 0;
     for (i = 0; i < nfft_6; i++) {
@@ -4916,7 +4916,7 @@ poisson_2s_ik(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
     fft1_6->compute(work1_6,work1_6,FFT3d::FORWARD);
 
   // two transforms when energies and pressures are calculated
-  
+
   } else {
     n = 0;
     for (i = 0; i < nfft_6; i++) {
@@ -5172,7 +5172,7 @@ poisson_none_ik(int n1, int n2,FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
         work2_6[n+1] = 0.5*(fky_6[j]-fky2_6[j])*work1_6[n];
         n += 2;
       }
- 
+
   fft2_6->compute(work2_6,work2_6,FFT3d::BACKWARD);
 
   n = 0;
@@ -5295,9 +5295,9 @@ poisson_2s_ad(FFT_SCALAR* dfft_1, FFT_SCALAR* dfft_2,
         n += 2;
       }
     }
-    
+
     // unify the two transformed vectors for efficient calculations later
-    
+
     for (i = 0; i < 2*nfft_6; i++)
       work1_6[i] += work2_6[i];
   }
@@ -5578,7 +5578,7 @@ poisson_none_peratom(int n1, int n2,
   int n,i,j,k;
 
   // compute first virial term v0
-  
+
   n = 0;
   for (i = 0; i < nfft_6; i++) {
     work2_6[n] = work1_6[n]*vg_6[i][0];
@@ -5962,7 +5962,7 @@ void PPPMDisp::fieldforce_g_ik()
     }
 
     // convert E-field to force
-    
+
     type = atom->type[i];
     lj = B[type];
     f[i][0] += lj*ekx;
@@ -6041,7 +6041,7 @@ void PPPMDisp::fieldforce_g_ad()
     ekz *= hz_inv;
 
     // convert E-field to force
-    
+
     type = atom->type[i];
     lj = B[type];
 
@@ -6124,7 +6124,7 @@ void PPPMDisp::fieldforce_g_peratom()
     }
 
     // convert E-field to force
-    
+
     type = atom->type[i];
     lj = B[type]*0.5;
 
@@ -6799,13 +6799,13 @@ void PPPMDisp::fieldforce_none_peratom()
         }
       }
     }
-    
+
     // convert D-field to force
- 
+
     type = atom->type[i];
     for (k = 0; k < nsplit; k++) {
       lj = B[nsplit*type + k]*0.5;
- 
+
       if (eflag_atom) {
         eatom[i] += u_pa[k]*lj;
       }
@@ -7878,7 +7878,7 @@ void PPPMDisp::pack_reverse_grid(int flag, void *vbuf, int nlist, int *list)
     FFT_SCALAR *src = &density_brick_g[nzlo_out_6][nylo_out_6][nxlo_out_6];
     for (int i = 0; i < nlist; i++)
       buf[i] = src[list[i]];
-    
+
   // dispersion interactions, arithmetic mixing
 
   } else if (flag == REVERSE_RHO_ARITH) {

src/KSPACE/pppm_disp.h

@@ -96,7 +96,7 @@ class PPPMDisp : public KSpace {
   int ngrid_6,nfft_6,nfft_both_6;
 
   // the following variables are needed for every structure factor
-  
+
   FFT_SCALAR ***density_brick;
   FFT_SCALAR ***vdx_brick,***vdy_brick,***vdz_brick;
   FFT_SCALAR *density_fft;

src/MANYBODY/pair_eam.cpp

@@ -63,6 +63,8 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
   z2r = nullptr;
   scale = nullptr;
 
+  rhomax = rhomin = 0.0;
+
   frho_spline = nullptr;
   rhor_spline = nullptr;
   z2r_spline = nullptr;

src/MANYBODY/pair_eam.h

@@ -42,7 +42,7 @@ class PairEAM : public Pair {
 
   // potentials in spline form used for force computation
 
-  double dr,rdr,drho,rdrho,rhomax;
+  double dr,rdr,drho,rdrho,rhomax,rhomin;
   double ***rhor_spline,***frho_spline,***z2r_spline;
 
   PairEAM(class LAMMPS *);

src/MANYBODY/pair_eam_fs.cpp

@@ -34,6 +34,7 @@ PairEAMFS::PairEAMFS(LAMMPS *lmp) : PairEAM(lmp)
 {
   one_coeff = 1;
   manybody_flag = 1;
+  he_flag = 0;
 }
 
 /* ----------------------------------------------------------------------
@@ -118,7 +119,8 @@ void PairEAMFS::read_file(char *filename)
 
   // read potential file
   if (comm->me == 0) {
-    PotentialFileReader reader(lmp, filename, "eam/fs", unit_convert_flag);
+    PotentialFileReader reader(lmp, filename, he_flag ? "eam/he" : "eam/fs",
+                               unit_convert_flag);
 
     // transparently convert units for supported conversions
 
@@ -149,12 +151,14 @@ void PairEAMFS::read_file(char *filename)
 
       //
 
-      values = reader.next_values(5);
+      if (he_flag) values = reader.next_values(6);
+      else values = reader.next_values(5);
       file->nrho = values.next_int();
       file->drho = values.next_double();
       file->nr   = values.next_int();
       file->dr   = values.next_double();
       file->cut  = values.next_double();
+      if (he_flag) rhomax = values.next_double();
 
       if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
         error->one(FLERR,"Invalid EAM potential file");
@@ -202,6 +206,7 @@ void PairEAMFS::read_file(char *filename)
   MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
   MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
   MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+  MPI_Bcast(&rhomax, 1, MPI_DOUBLE, 0, world);
 
   // allocate memory on other procs
   if (comm->me != 0) {
@@ -255,7 +260,8 @@ void PairEAMFS::file2array()
   nr = fs->nr;
   drho = fs->drho;
   dr = fs->dr;
-  rhomax = (nrho-1) * drho;
+  if (he_flag) rhomin = rhomax - (nrho-1) * drho;
+  else rhomax = (nrho-1) * drho;
 
   // ------------------------------------------------------------------
   // setup frho arrays

src/MANYBODY/pair_eam_fs.h

@@ -35,6 +35,7 @@ class PairEAMFS : virtual public PairEAM {
  protected:
   void read_file(char *);
   void file2array();
+  int he_flag;
 };
 
 }

src/MANYBODY/pair_eam_he.cpp

@@ -0,0 +1,238 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://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: Xiaowng Zhou (Sandia)
+------------------------------------------------------------------------- */
+
+#include "pair_eam_he.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairEAMHE::PairEAMHE(LAMMPS *lmp) : PairEAMFS(lmp), PairEAM(lmp)
+{
+  he_flag = 1;
+}
+
+void PairEAMHE::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,m,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+  double rsq,r,p,rhoip,rhojp,z2,z2p,recip,phip,psip,phi;
+  double *coeff;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = 0.0;
+  ev_init(eflag,vflag);
+
+  // grow energy and fp arrays if necessary
+  // need to be atom->nmax in length
+
+  if (atom->nmax > nmax) {
+    memory->destroy(rho);
+    memory->destroy(fp);
+    memory->destroy(numforce);
+    nmax = atom->nmax;
+    memory->create(rho,nmax,"pair:rho");
+    memory->create(fp,nmax,"pair:fp");
+    memory->create(numforce,nmax,"pair:numforce");
+  }
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // zero out density
+
+  if (newton_pair) {
+    for (i = 0; i < nall; i++) rho[i] = 0.0;
+  } else for (i = 0; i < nlocal; i++) rho[i] = 0.0;
+
+  // rho = density at each atom
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    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 < cutforcesq) {
+        jtype = type[j];
+        p = sqrt(rsq)*rdr + 1.0;
+        m = static_cast<int> (p);
+        m = MIN(m,nr-1);
+        p -= m;
+        p = MIN(p,1.0);
+        coeff = rhor_spline[type2rhor[jtype][itype]][m];
+        rho[i] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+        if (newton_pair || j < nlocal) {
+          coeff = rhor_spline[type2rhor[itype][jtype]][m];
+          rho[j] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+        }
+      }
+    }
+  }
+
+  // communicate and sum densities
+
+  if (newton_pair) comm->reverse_comm_pair(this);
+
+  // fp = derivative of embedding energy at each atom
+  // phi = embedding energy at each atom
+  // if rho > rhomax or rho < rhomin (i.e., table is exceeded),
+  //   add linear term to conserve energy
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    p = (rho[i]-rhomin)*rdrho + 1.0;
+    m = static_cast<int> (p);
+    if (m < 2) {
+      m = 2;
+    } else if (m > nrho-1) {
+      m = nrho-1;
+    }
+    p -= m;
+    if (p < -1.0) {
+      p = -1.0;
+    } else if (p > 1.0) {
+      p = 1.0;
+    }
+    coeff = frho_spline[type2frho[type[i]]][m];
+    fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2];
+    if (eflag) {
+      phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+      if (rho[i] < rhomin) {
+        phi += fp[i] * (rho[i]-rhomin);
+      } else if (rho[i] > rhomax) {
+        phi += fp[i] * (rho[i]-rhomax);
+      }
+      phi *= scale[type[i]][type[i]];
+      if (eflag_global) eng_vdwl += phi;
+      if (eflag_atom) eatom[i] += phi;
+    }
+  }
+
+  // communicate derivative of embedding function
+
+  comm->forward_comm_pair(this);
+  embedstep = update->ntimestep;
+
+  // compute forces on each atom
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    numforce[i] = 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 < cutforcesq) {
+        ++numforce[i];
+        jtype = type[j];
+        r = sqrt(rsq);
+        p = r*rdr + 1.0;
+        m = static_cast<int> (p);
+        m = MIN(m,nr-1);
+        p -= m;
+        p = MIN(p,1.0);
+
+        // rhoip = derivative of (density at atom j due to atom i)
+        // rhojp = derivative of (density at atom i due to atom j)
+        // phi = pair potential energy
+        // phip = phi'
+        // z2 = phi * r
+        // z2p = (phi * r)' = (phi' r) + phi
+        // psip needs both fp[i] and fp[j] terms since r_ij appears in two
+        //   terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
+        //   hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
+        // scale factor can be applied by thermodynamic integration
+
+        coeff = rhor_spline[type2rhor[itype][jtype]][m];
+        rhoip = (coeff[0]*p + coeff[1])*p + coeff[2];
+        coeff = rhor_spline[type2rhor[jtype][itype]][m];
+        rhojp = (coeff[0]*p + coeff[1])*p + coeff[2];
+        coeff = z2r_spline[type2z2r[itype][jtype]][m];
+        z2p = (coeff[0]*p + coeff[1])*p + coeff[2];
+        z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+
+        recip = 1.0/r;
+        phi = z2*recip;
+        phip = z2p*recip - phi*recip;
+        psip = fp[i]*rhojp + fp[j]*rhoip + phip;
+        fpair = -scale[itype][jtype]*psip*recip;
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) evdwl = scale[itype][jtype]*phi;
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,0.0,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+

src/MANYBODY/pair_eam_he.h

@@ -0,0 +1,65 @@
+/* -*- 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 PAIR_CLASS
+
+PairStyle(eam/he,PairEAMHE)
+
+#else
+
+#ifndef LMP_PAIR_EAM_HE_H
+#define LMP_PAIR_EAM_HE_H
+
+#include "pair_eam_fs.h"
+
+namespace LAMMPS_NS {
+
+class PairEAMHE : public PairEAMFS {
+ public:
+  PairEAMHE(class LAMMPS *);
+  virtual ~PairEAMHE() {}
+
+ protected:
+  void compute(int, int);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: No matching element in EAM potential file
+
+The EAM potential file does not contain elements that match the
+requested elements.
+
+E: Cannot open EAM potential file %s
+
+The specified EAM potential file cannot be opened.  Check that the
+path and name are correct.
+
+E: Incorrect element names in EAM potential file
+
+The element names in the EAM file do not match those requested.
+
+E: Invalid EAM potential file
+
+UNDOCUMENTED
+
+*/

unittest/force-styles/tests/atomic-pair-eam_he.yaml

@@ -0,0 +1,89 @@
+---
+lammps_version: 24 Dec 2020
+date_generated: Wed Jan 13 22:16:02 202
+epsilon: 5e-12
+prerequisites: ! |
+  pair eam/he
+pre_commands: ! ""
+post_commands: ! ""
+input_file: in.metal
+pair_style: eam/he
+pair_coeff: ! |
+  * * PdHHe.eam.he Pd He
+extract: ! ""
+natoms: 32
+init_vdwl: -15.3146152609365
+init_coul: 0
+init_stress: ! |2-
+   1.0127679615895376e+02  9.3559970433415444e+01  9.1432412459889193e+01  4.0035576925472673e+00 -9.7686923241135581e-01  2.8980241224200443e+00
+init_forces: ! |2
+    1  1.2005847856522014e+00  2.2040396377071869e+00 -1.0711805842453261e+00
+    2 -1.5935390465573379e-01 -8.3057674140080351e-01  3.5885201320339977e-01
+    3 -8.3922102051696701e-02  4.1053523947712707e+00  3.2379570522514384e-01
+    4  9.6680271004438867e-01 -6.0272944088190550e-01  1.7088406391974090e-01
+    5 -1.5118032750014099e-01  4.8152982355843301e+00 -2.3396251503834162e-01
+    6  6.9769996639420007e-01  2.1986594375389688e+00  3.7782787462812606e-01
+    7 -1.9293747297385899e+00  2.3965302429887498e+00  5.3526894592078611e-01
+    8 -5.7898968354416502e-01 -1.7424269600946102e-01  4.6553500563912831e-01
+    9 -1.9170501143342826e+00 -1.5811474204082574e+00 -1.5153955037081890e+00
+   10 -1.0389874633978984e+00 -1.0136677465869111e+00 -2.8606245342679020e+00
+   11 -6.8937210559650852e-01 -4.9861949626741522e+00  1.7468278589450823e+00
+   12 -1.0247245841335400e+00 -2.9144176183755111e+00  1.0162869145592908e+00
+   13  8.8392380716960872e-01 -2.1032444766660849e-01 -4.5408294970791102e-01
+   14  3.9708226038326155e-01 -7.8680161984030961e-01 -1.9977901194159892e-01
+   15 -7.8687732774064545e-02 -2.1157301146984339e-01  3.4042801915998766e-01
+   16  4.0325457730206029e+00 -2.3548979291247609e+00  9.2949967143894952e-01
+   17  7.4997446543261548e-01  2.0845833390037725e+00  1.7238817466288217e+00
+   18 -2.1412087035667221e-01 -5.6906054172322229e-01 -5.2781467006833294e-01
+   19 -3.0256742141254084e-01  6.0688322127888294e-01 -9.1127162282408275e-02
+   20  2.3174481031809935e-01  3.0892939020181726e-01 -3.7137738763066941e-01
+   21  1.1702211057625094e+00  4.2920154821923315e+00  1.3460541385609648e+00
+   22  3.8027613826247031e-02  4.3239633632972230e-01 -3.8188409283423194e-02
+   23  1.4000432801696054e+00  1.0040601640391840e+00 -2.4122350019076917e+00
+   24 -1.5604155772955447e-01  3.4572668285914510e-01 -2.8556703863036959e-01
+   25 -2.9449464597969849e-01 -3.4630638648128692e-01  1.1805865362173559e-01
+   26 -1.8108866036308173e+00 -1.8950909756776353e+00  3.3196635723271326e+00
+   27 -7.7538420123902196e-01 -1.2937697587184989e+00 -9.6725767253143236e-01
+   28 -3.5707629120214823e-01 -2.8995606245962768e-01 -1.2007211500278167e-01
+   29  3.6987879522593697e-01 -4.9287949481541249e-01  5.4972323630766012e-02
+   30  1.1712105233973889e-01 -6.9110122964840481e-01  9.5437848811806628e-02
+   31  1.9388288555816860e-01 -2.0146460156127194e-01 -2.0863139798712499e-01
+   32 -8.8731897202011578e-01 -3.3482718789714783e+00 -1.5659784019873610e+00
+run_vdwl: -15.3236335310355
+run_coul: 0
+run_stress: ! |2-
+   1.0125543392557182e+02  9.3539230810233988e+01  9.1388997082229878e+01  4.0040941706030253e+00 -9.7826716756924303e-01  2.9018476991088571e+00
+run_forces: ! |2
+    1  1.1949883047712304e+00  2.2060858143622002e+00 -1.0733885545165418e+00
+    2 -1.5963527935293531e-01 -8.3496477000211577e-01  3.6354718487355586e-01
+    3 -8.4458087648033864e-02  4.1335068175946956e+00  3.2562083254074076e-01
+    4  9.7516402928123347e-01 -6.1050193540570252e-01  1.7035859794498676e-01
+    5 -1.5399721579534215e-01  4.8120569088649683e+00 -2.3294021119313149e-01
+    6  6.9813149221746262e-01  2.1977391468237610e+00  3.7752101367086355e-01
+    7 -1.9335938530288574e+00  2.3989898112356141e+00  5.3226592065468015e-01
+    8 -5.8346129298071359e-01 -1.7815682091755050e-01  4.6582930751668622e-01
+    9 -1.9172329351091069e+00 -1.5829627245185132e+00 -1.5163042130781048e+00
+   10 -1.0334176082685609e+00 -1.0148545000176199e+00 -2.8584769425899079e+00
+   11 -6.8756244404382794e-01 -4.9891360565883884e+00  1.7443920243481270e+00
+   12 -1.0236334662680082e+00 -2.9150791430800980e+00  1.0158839648906106e+00
+   13  8.8312645129364897e-01 -2.1071143805836662e-01 -4.5244072317123446e-01
+   14  3.9927871470591697e-01 -7.8979672016550584e-01 -2.0106817979273284e-01
+   15 -7.8923521050882781e-02 -2.1032668862489418e-01  3.3903131601122610e-01
+   16  4.0274555753165444e+00 -2.3555929271854597e+00  9.2684665883544881e-01
+   17  7.5117151275169469e-01  2.0849083749786277e+00  1.7261549167735377e+00
+   18 -2.1494836879728668e-01 -5.7509260272248663e-01 -5.3373034744908598e-01
+   19 -3.0249976756523766e-01  6.0674463278548640e-01 -9.1809428603960672e-02
+   20  2.3156742504597116e-01  3.0887145855490367e-01 -3.7127049070468199e-01
+   21  1.1686499084188677e+00  4.2891627780548545e+00  1.3459230037755703e+00
+   22  3.8179702424449569e-02  4.3243601345801236e-01 -3.8079863481504578e-02
+   23  1.4067770990537660e+00  1.0070141778822215e+00 -2.4069585948142982e+00
+   24 -1.5591693014836788e-01  3.4568121542161234e-01 -2.8546260901671355e-01
+   25 -2.9429104311734711e-01 -3.4618233584978003e-01  1.1804410855998390e-01
+   26 -1.8072996525588660e+00 -1.8933503719653344e+00  3.3208542120004427e+00
+   27 -7.7461160804433704e-01 -1.2927908441630280e+00 -9.6864680654489232e-01
+   28 -3.5780846246953718e-01 -2.8988724628268286e-01 -1.1924240028778682e-01
+   29  3.6970648233559189e-01 -4.9268726775164107e-01  5.5185028378882096e-02
+   30  1.1705695741665428e-01 -6.9122064837750197e-01  9.5592509524696681e-02
+   31  1.9349589373949760e-01 -1.9991899011439837e-01 -2.0661879538790651e-01
+   32 -8.9145801252527535e-01 -3.3499831182258872e+00 -1.5666124396675569e+00
+...

unittest/force-styles/tests/atomic-pair-eam_he_real.yaml

@@ -0,0 +1,90 @@
+---
+lammps_version: 24 Dec 2020
+date_generated: Wed Jan 13 22:16:02 202
+epsilon: 7.5e-12
+prerequisites: ! |
+  pair eam/he
+pre_commands: ! |
+  variable units index real
+post_commands: ! ""
+input_file: in.metal
+pair_style: eam/he
+pair_coeff: ! |
+  * * PdHHe.eam.he Pd He
+extract: ! ""
+natoms: 32
+init_vdwl: -353.163435640974
+init_coul: 0
+init_stress: ! |2-
+   2.3354985203865667e+03  2.1575442826183294e+03  2.1084816277194832e+03  9.2324238343315059e+01 -2.2527140800613445e+01  6.6830027278251166e+01
+init_forces: ! |2
+    1  2.7686144278186653e+01  5.0826364063289034e+01 -2.4702012350837244e+01
+    2 -3.6747885266547895e+00 -1.9153555643332677e+01  8.2753244342250110e+00
+    3 -1.9352897465461538e+00  9.4671680061884146e+01  7.4669067263334208e+00
+    4  2.2295001268309708e+01 -1.3899271805198000e+01  3.9406803293404700e+00
+    5 -3.4863013501526918e+00  1.1104342091130563e+02 -5.3953040422049119e+00
+    6  1.6089344262331657e+01  5.0702293693679884e+01  8.7129182164280419e+00
+    7 -4.4492440494497494e+01  5.5265303098423985e+01  1.2343595755584241e+01
+    8 -1.3351819967863772e+01 -4.0181322292175965e+00  1.0735492808756252e+01
+    9 -4.4208228097061294e+01 -3.6462127564548197e+01 -3.4945852267642280e+01
+   10 -2.3959621310073842e+01 -2.3375734739886113e+01 -6.5967572243087659e+01
+   11 -1.5897299220341502e+01 -1.1498439326030027e+02  4.0282809435768357e+01
+   12 -2.3630711483916183e+01 -6.7208070295011467e+01  2.3436134191253181e+01
+   13  2.0383768267501306e+01 -4.8501972313137021e+00 -1.0471402111803629e+01
+   14  9.1569349225997474e+00 -1.8144077307608963e+01 -4.6070136940518891e+00
+   15 -1.8145823173352931e+00 -4.8789897980772752e+00  7.8504570168111005e+00
+   16  9.2992719393484805e+01 -5.4305239084580101e+01  2.1434772718702309e+01
+   17  1.7294822908857860e+01  4.8071636233680366e+01  3.9753659488339430e+01
+   18 -4.9377448227825411e+00 -1.3122848506373817e+01 -1.2171696062028875e+01
+   19 -6.9773708472871210e+00  1.3995060261579386e+01 -2.1014423910442734e+00
+   20  5.3441625538361990e+00  7.1240813402889511e+00 -8.5641664449490893e+00
+   21  2.6985941150269763e+01  9.8976233335854843e+01  3.1040747418937880e+01
+   22  8.7693765199327156e-01  9.9712969013523196e+00 -8.8064568351232408e-01
+   23  3.2285766664472092e+01  2.3154178611774061e+01 -5.5627463461007146e+01
+   24 -3.5984039880587519e+00  7.9726471106807173e+00 -6.5853326871205784e+00
+   25 -6.7912082138528707e+00 -7.9860153944650030e+00  2.7224973667181924e+00
+   26 -4.1760039256471202e+01 -4.3701838304071202e+01  7.6553264473163750e+01
+   27 -1.7880785366498504e+01 -2.9835040915646299e+01 -2.2305492953037302e+01
+   28 -8.2343753100058397e+00 -6.6865459861976770e+00 -2.7689288915546530e+00
+   29  8.5296080813687887e+00 -1.1366071741285882e+01  1.2676919627309815e+00
+   30  2.7008757664122220e+00 -1.5937173770267322e+01  2.2008491889792485e+00
+   31  4.4710457826753673e+00 -4.6458843160683996e+00 -4.8111545762194012e+00
+   32 -2.0462062632899592e+01 -7.7212987730343386e+01 -3.6112321671970648e+01
+run_vdwl: -353.230598025017
+run_coul: 0
+run_stress: ! |2-
+   2.3352018947327042e+03  2.1575803047701361e+03  2.1079997407908877e+03  9.2286943467903825e+01 -2.2591362673321409e+01  6.6981941401935686e+01
+run_forces: ! |2
+    1  2.7570840184953529e+01  5.0878936563954468e+01 -2.4745193716885819e+01
+    2 -3.6751057733124783e+00 -1.9167341173087838e+01  8.2931871358142466e+00
+    3 -1.9458273195805660e+00  9.4745013001500453e+01  7.4613246126708113e+00
+    4  2.2322999682742939e+01 -1.3910836126089546e+01  3.9444558071725182e+00
+    5 -3.5611755339443869e+00  1.1098583048019351e+02 -5.3589546434886657e+00
+    6  1.6098529359830774e+01  5.0681471667414634e+01  8.7065695410635300e+00
+    7 -4.4532333454110486e+01  5.5329574524982760e+01  1.2275125678150991e+01
+    8 -1.3365798458276837e+01 -4.0311665119688325e+00  1.0738893622536514e+01
+    9 -4.4219909611484340e+01 -3.6506027825529728e+01 -3.4970646977778003e+01
+   10 -2.3838805073020580e+01 -2.3388505371089686e+01 -6.5920271795280030e+01
+   11 -1.5863107119046816e+01 -1.1501474143750076e+02  4.0236861649778703e+01
+   12 -2.3597639915238879e+01 -6.7219256881236532e+01  2.3428611972869568e+01
+   13  2.0372282061050672e+01 -4.8599967847476790e+00 -1.0435173423332810e+01
+   14  9.1868968538149876e+00 -1.8180822851027013e+01 -4.6226879613716969e+00
+   15 -1.8187982708833512e+00 -4.8633248636475823e+00  7.8331412550634454e+00
+   16  9.2883964220226872e+01 -5.4312212635591663e+01  2.1376066838247318e+01
+   17  1.7335571325471630e+01  4.8086629272887443e+01  3.9799611586568851e+01
+   18 -4.9380025631152833e+00 -1.3140883981374724e+01 -1.2192856535881621e+01
+   19 -6.9765790243039625e+00  1.3993117559236460e+01 -2.1082447485751517e+00
+   20  5.3413030155749288e+00  7.1244381338347607e+00 -8.5630059362592803e+00
+   21  2.6939069851048053e+01  9.8929039552947501e+01  3.1025635765431552e+01
+   22  8.7771899094057326e-01  9.9759331856305540e+00 -8.7675403015085085e-01
+   23  3.2451609933852048e+01  2.3224549699639688e+01 -5.5509981817200973e+01
+   24 -3.5964621834607793e+00  7.9732311561843838e+00 -6.5839559665478173e+00
+   25 -6.7879787318782503e+00 -7.9847098239766732e+00  2.7230073490545976e+00
+   26 -4.1688921616984544e+01 -4.3653314583134815e+01  7.6583618010794297e+01
+   27 -1.7863754347437737e+01 -2.9811723761490370e+01 -2.2339924855988965e+01
+   28 -8.2430574229283415e+00 -6.6852368720976481e+00 -2.7583576098041354e+00
+   29  8.5269603298722565e+00 -1.1364005508745327e+01  1.2731276160949609e+00
+   30  2.6994843727982003e+00 -1.5940803816308822e+01  2.2046514286792602e+00
+   31  4.4654048480426765e+00 -4.6271939008948024e+00 -4.7862215560852475e+00
+   32 -2.0559378611212470e+01 -7.7265660088866639e+01 -3.6131658295360147e+01
+...