1 /* ----------------------------------------------------------------------
2    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3    http://lammps.sandia.gov, Sandia National Laboratories
4    Steve Plimpton, sjplimp@sandia.gov
5 
6    Copyright (2003) Sandia Corporation.  Under the terms of Contract
7    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8    certain rights in this software.  This software is distributed under
9    the GNU General Public License.
10 
11    See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13 
14 /* ----------------------------------------------------------------------
15    Modified for use with KIM by Tobias Brink (2020).
16 ------------------------------------------------------------------------- */
17 
18 
19 #ifndef LMP_PAIR_TERSOFF_ZBL_H
20 #define LMP_PAIR_TERSOFF_ZBL_H
21 
22 #include <string>
23 #include <map>
24 #include <istream>
25 
26 #include "KIM_ModelDriverHeaders.hpp"
27 
28 #include "pair_tersoff.hpp"
29 #include "ndarray.hpp"
30 
31 namespace model_driver_Tersoff {
32 
33 
34 class PairTersoffZBL : public PairTersoff {
35  public:
36   // Parameters are stored together in a struct because I found that
37   // using the layout that conforms to KIM (one array per parameter)
38   // can lead to a slowdown of up to 5% because the data is not in
39   // adjacent memory.  The parameter data is not so much, so we will
40   // copy between the KIM-published parameters and this internal data.
41   struct ParamsZBL2 {
42     // Two-body parameters.
43     double ZBLcut;
44     double ZBLexpscale;
45     // Pre-computed.
46     double a;       // 0.8854 * a0 / (Z_i^0.23 + Z_j^0.23)
47     double premult; // Z_i * Z_j * e^2 / (4 * pi * epsilon0)
48     // Not used in computation.
49     //double Z_i, Z_j;
50   };
51   struct KIMParamsZBL {
KIMParamsZBLmodel_driver_Tersoff::PairTersoffZBL::KIMParamsZBL52     explicit KIMParamsZBL(int N) // Number of particle types
53       : Z_i(N,N), Z_j(N,N),
54         ZBLcut(N,N), ZBLexpscale(N,N) {}
55     // Copy data from a Params array.
from_paramsmodel_driver_Tersoff::PairTersoffZBL::KIMParamsZBL56     void from_params(const Array2D<ParamsZBL2>& p2) {
57       for (int i = 0; i < Z_i.extent(0); ++i)
58         for (int j = 0; j < Z_i.extent(1); ++j) {
59           //Z_i(i,j) = p2(i,j).Z_i; // Those are not kept there,
60           //Z_j(i,j) = p2(i,j).Z_j; // but only derived in "a" and "premult"
61           ZBLcut(i,j) = p2(i,j).ZBLcut;
62           ZBLexpscale(i,j) = p2(i,j).ZBLexpscale;
63         }
64     };
65     // Copy data to a Params array.
to_paramsmodel_driver_Tersoff::PairTersoffZBL::KIMParamsZBL66     void to_params(Array2D<ParamsZBL2>& p2) const {
67       for (int i = 0; i < Z_i.extent(0); ++i)
68         for (int j = 0; j < Z_i.extent(1); ++j) {
69           //p2(i,j).Z_i = Z_i(i,j); // Those are not kept there,
70           //p2(i,j).Z_i = Z_i(i,j); // but only derived in "a" and "premult"
71           p2(i,j).ZBLcut = ZBLcut(i,j);
72           p2(i,j).ZBLexpscale = ZBLexpscale(i,j);
73         }
74     }
75     Array2D<double> Z_i, Z_j;
76     Array2D<double> ZBLcut, ZBLexpscale;
77     // The size of all parameter arrays. Needed for calls to
78     // KIM::ModelDriverCreate.SetParameterPointer().
79     //const int size2; -> can take this from PairTersoff::KIMParams!
80   };
81 
82   PairTersoffZBL(const std::string& parameter_file,
83                  int n_spec,
84                  std::map<std::string,int> type_map,
85                  // Conversion factors.
86                  double energy_conv,
87                  double inv_energy_conv,
88                  double length_conv,
89                  double inv_length_conv,
90                  double charge_conv);
91   virtual ~PairTersoffZBL();
92 
93   virtual void update_params(); // Copy from KIM-published parameters to internal.
94 
95   KIMParamsZBL kim_params_zbl; // ZBL parameters published to KIM, see
96                                // above why we keep two copies.
97 
98   void write_params(std::ofstream&); // Write parameters in the correct format.
99 
100  protected:
101   Array2D<ParamsZBL2> params_zbl_2; // n_spec*n_spec array of ZBL parameters
102 
103   void read_params(std::istream&, std::map<std::string,int>,
104                    double, double, double);
105   void prepare_params();
106 
107   virtual double repulsive(double, double, double, int, int,
108                            bool, double&) const;
109   virtual double ters_fa(double, double, int, int) const;
110   virtual double ters_fa_d(double, double, double, int, int) const;
111 
112  private:
113   const double global_a_0;       // Bohr radius for Coulomb repulsion
114   const double global_epsilon_0; // permittivity of vacuum for Coulomb repulsion
115   const double global_e;         // proton charge (negative of electron charge)
116 
117   // Derived, cached values.
118   const double global_e_sq;
119 
120   // Fermi-like smoothing from ZBL to Tersoff.
121   double F_fermi(double, double, double) const;
122   double F_fermi_d(double, double, double) const;
123 };
124 
125 }
126 
127 #endif
128