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