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
/* -*- 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(edip/multi,PairEDIPMulti)
#else
#ifndef LMP_PAIR_EDIP_MULTI_H
#define LMP_PAIR_EDIP_MULTI_H
#include "pair.h"
namespace LAMMPS_NS {
class PairEDIPMulti : public Pair {
public:
PairEDIPMulti(class LAMMPS *);
virtual ~PairEDIPMulti();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
protected:
struct Param {
double A, B;//coefficients for pair interaction I-J
double cutoffA;//cut-off distance for pair interaction I-J
double cutoffC;//lower cut-off distance for calculating Z_I
double alpha;//coefficient for calculating Z_I
double beta;//attractive term for pair I-J
double sigma;//cut-off coefficient for pair I-J
double rho;//pair I-J
double gamma;//coefficient for three-body interaction I-J-K
double eta, lambda;//coefficients for function h(l,Z)
double mu, Q0;//coefficients for function Q(Z)
double u1, u2, u3, u4;//coefficients for function tau(Z)
double cutsq;
int ielement,jelement,kelement;
};
double *preForceCoord;
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements
int ***elem2param; // mapping from element triplets to parameters
int *map; // mapping from atom types to elements
int nparams; // # of stored parameter sets
int maxparam; // max # of parameter sets
Param *params; // parameter set for an I-J-K interaction
// max number of interaction per atom for f(Z) environment potential
static const int leadDimInteractionList = 64;
void allocate();
void allocatePreLoops(void);
void deallocatePreLoops(void);
void read_file(char *);
void setup();
void edip_pair(double, double, Param *, double &, double &, double &);
void edip_fc(double, Param *, double &, double &);
void edip_fcut2(double, Param *, double &, double &);
void edip_tau(double, Param *, double &, double &);
void edip_h(double, double, Param *, double &, double &, double &);
void edip_fcut3(double, Param *, double &, double &);
double vec3_dot(double x[3], double y[3])
{
return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
}
void vec3_add(double k1, double x[3], double k2, double y[3], double *z)
{
z[0] = k1 * x[0] + k2 * y[0];
z[1] = k1 * x[1] + k2 * y[1];
z[2] = k1 * x[2] + k2 * y[2];
}
//dr_ij=r_j - r_i
//dr_ik=r_k - r_i
void costheta_d(double *dr_ij, double r_ij, double *dr_ik, double r_ik,
double *dri, double *drj, double *drk)
{
double costheta;
costheta = vec3_dot(dr_ij, dr_ik) / r_ij / r_ik;
vec3_add(1 / r_ij / r_ik, dr_ik, -costheta / r_ij / r_ij, dr_ij, drj);
vec3_add(1 / r_ij / r_ik, dr_ij, -costheta / r_ik / r_ik, dr_ik, drk);
vec3_add(-1, drj, -1, drk, dri);
}
};
}
#endif
#endif