ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add standard version of ilp_water_2dm

Browse Source
This commit is contained in:
oywg11
2023-05-18 21:51:53 +08:00
parent 10880332ed
commit 81a497adcd
6 changed files with 373 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
src/INTERLAYER/pair_ilp_graphene_hbn.cpp
View File

@ -1,7 +1,7 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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
@ -47,11 +47,11 @@ using namespace InterLayer;
static const char cite_ilp[] =
"ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848\n"
"@Article{Ouyang2018\n"
" author = {W. Ouyang and D. Mandelli and M. Urbakh and O. Hod},\n"
" author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},\n"
" title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},\n"
" journal = {Nano Letters},\n"
" volume = 18,\n"
" pages = 6009,\n"
" pages = {6009}\n"
" year = 2018,\n"
"}\n\n";
@ -59,7 +59,10 @@ static const char cite_ilp[] =
static std::map<int, std::string> variant_map = {
{PairILPGrapheneHBN::ILP_GrhBN, "ilp/graphene/hbn"},
{PairILPGrapheneHBN::ILP_TMD, "ilp/tmd"},
{PairILPGrapheneHBN::SAIP_METAL, "saip/metal"}};
{PairILPGrapheneHBN::ILP_PHOSPHORUS, "ilp/phosphorus"},
{PairILPGrapheneHBN::ILP_WATER_2DM, "ilp/water/2dm"},
{PairILPGrapheneHBN::SAIP_METAL, "saip/metal"},
{PairILPGrapheneHBN::SAIP_METAL_TMD, "saip/metal/tmd"}};
/* ---------------------------------------------------------------------- */
@ -313,14 +316,14 @@ void PairILPGrapheneHBN::read_file(char *filename)
for (int m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement) {
if (n >= 0)
error->all(FLERR, "{} potential file {} has a duplicate entry for: {} {}",
variant_map[variant], filename, elements[i], elements[j]);
error->all(FLERR, "{} potential file {} has a duplicate entry", variant_map[variant],
filename);
n = m;
}
}
if (n < 0)
error->all(FLERR, "{} potential file {} is missing an entry for: {} {}",
variant_map[variant], filename, elements[i], elements[j]);
error->all(FLERR, "{} potential file {} is missing an entry", variant_map[variant],
filename);
elem2param[i][j] = n;
cutILPsq[i][j] = params[n].rcut * params[n].rcut;
}
@ -632,7 +635,7 @@ void PairILPGrapheneHBN::calc_FRep(int eflag, int /* vflag */)
void PairILPGrapheneHBN::ILP_neigh()
{
int i, j, ii, jj, n, allnum, jnum, itype, jtype;
int i, j, ii, jj, n, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
int *ilist, *jlist, *numneigh, **firstneigh;
int *neighptr;
@ -649,7 +652,7 @@ void PairILPGrapheneHBN::ILP_neigh()
(int **) memory->smalloc(maxlocal * sizeof(int *), "ILPGrapheneHBN:firstneigh");
}
allnum = list->inum + list->gnum;
inum = list->inum; // + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -659,7 +662,7 @@ void PairILPGrapheneHBN::ILP_neigh()
ipage->reset();
for (ii = 0; ii < allnum; ii++) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
n = 0;

38
src/INTERLAYER/pair_ilp_graphene_hbn.h
View File

@ -1,7 +1,7 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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
@ -35,11 +35,19 @@ class PairILPGrapheneHBN : public Pair {
double init_one(int, int) override;
void init_style() override;
void calc_FvdW(int, int);
virtual void calc_FRep(int, int);
virtual void ILP_neigh();
virtual void calc_normal();
void read_file(char *);
void allocate();
double single(int, int, int, int, double, double, double, double &) override;
static constexpr int NPARAMS_PER_LINE = 13;
enum { ILP_GrhBN, ILP_TMD, SAIP_METAL }; // for telling class variants apart in shared code
// for telling class variants apart in shared code
enum { ILP_GrhBN, ILP_TMD, SAIP_METAL,
ILP_WATER_2DM, SAIP_METAL_TMD, ILP_PHOSPHORUS };
protected:
int me;
@ -51,6 +59,7 @@ class PairILPGrapheneHBN : public Pair {
int *ILP_numneigh; // # of pair neighbors for each atom
int **ILP_firstneigh; // ptr to 1st neighbor of each atom
int tap_flag; // flag to turn on/off taper function
double ncf; // for ilp/phosphorus, coefficients for calcualting the normals
struct Param {
double z0, alpha, epsilon, C, delta, d, sR, reff, C6, S;
@ -76,15 +85,28 @@ class PairILPGrapheneHBN : public Pair {
double ***dpvet1;
double ***dpvet2;
double ***dNave;
virtual void ILP_neigh();
virtual void calc_normal();
virtual void calc_FRep(int, int);
void read_file(char *);
void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/

213
src/INTERLAYER/pair_ilp_tmd.cpp
View File

@ -1,7 +1,7 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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
@ -22,12 +22,14 @@
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "interlayer_taper.h"
#include "memory.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
@ -35,22 +37,20 @@
using namespace LAMMPS_NS;
using namespace InterLayer;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
static const char cite_ilp_tmd[] =
"ilp/tmd potential doi:10.1021/acs.jctc.1c00782\n"
"@Article{Ouyang2021\n"
" author = {W. Ouyang and R. Sofer and X. Gao and J. Hermann and\n"
" A. Tkatchenko and L. Kronik and M. Urbakh and O. Hod},\n"
" title = {Anisotropic Interlayer Force Field for Transition\n"
" Metal Dichalcogenides: The Case of Molybdenum Disulfide},\n"
" journal = {J.~Chem.\\ Theory Comput.},\n"
" volume = 17,\n"
" pages = {7237--7245}\n"
" year = 2021,\n"
"}\n\n";
static const char cite_ilp_tmd[] = "ilp/tmd potential doi/10.1021/acs.jctc.1c00782\n"
"@Article{Ouyang2021\n"
" author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. "
"Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},\n"
" title = {Anisotropic Interlayer Force Field for Transition "
"Metal Dichalcogenides: The Case of Molybdenum Disulfide},\n"
" journal = {J. Chem. Theory Comput.},\n"
" volume = 17,\n"
" pages = {7237–7245}\n"
" year = 2021,\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
@ -232,7 +232,7 @@ void PairILPTMD::calc_FRep(int eflag, int /* vflag */)
void PairILPTMD::ILP_neigh()
{
int i, j, l, ii, jj, ll, n, allnum, jnum, itype, jtype, ltype, imol, jmol, count;
int i, j, l, ii, jj, ll, n, inum, jnum, itype, jtype, ltype, imol, jmol, count;
double xtmp, ytmp, ztmp, delx, dely, delz, deljx, deljy, deljz, rsq, rsqlj;
int *ilist, *jlist, *numneigh, **firstneigh;
int *neighsort;
@ -249,7 +249,7 @@ void PairILPTMD::ILP_neigh()
ILP_firstneigh = (int **) memory->smalloc(maxlocal * sizeof(int *), "ILPTMD:firstneigh");
}
allnum = list->inum + list->gnum;
inum = list->inum; //+ list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -259,7 +259,7 @@ void PairILPTMD::ILP_neigh()
ipage->reset();
for (ii = 0; ii < allnum; ii++) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
//initialize varibles
@ -288,21 +288,21 @@ void PairILPTMD::ILP_neigh()
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// check if the atom i is TMD, i.e., Mo/S/W/Se
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
// check if the atom i is a TMD atom, i.e., Mo/S/W/Se
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
strcmp(elements[itype], "Te") == 0) {
if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol && type[i] == type[j]) {
neighptr[n++] = j;
}
} else { // atom i is C, B, N or H.
} else { // atom i can be P, C, B, N or H.
if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol) { neighptr[n++] = j; }
}
} // loop over jj
// if atom i is Mo/W/S/Se/Te, then sorting the orders of neighbors
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
strcmp(elements[itype], "Te") == 0) {
// initialize neighsort
for (ll = 0; ll < n; ll++) {
@ -367,11 +367,11 @@ void PairILPTMD::ILP_neigh()
ll++;
}
} // end of sorting the order of neighbors
} else { // for B/N/C/H atoms
} else { // for P/B/N/C/H atoms
if (n > 3)
error->one(
FLERR,
"There are too many neighbors for B/N/C/H atoms, please check your configuration");
"There are too many neighbors for P/B/N/C/H atoms, please check your configuration");
for (ll = 0; ll < n; ll++) { neighsort[ll] = neighptr[ll]; }
}
@ -390,12 +390,14 @@ void PairILPTMD::calc_normal()
{
int i, j, ii, jj, inum, jnum;
int cont, id, ip, m, k, itype;
double nn, xtp, ytp, ztp, delx, dely, delz, nn2;
int *ilist, *jlist;
int iH1,iH2,jH1,jH2;
double Nave[3], dni[3], dpvdri[3][3];
double nn, xtp, ytp, ztp, delx, dely, delz, nn2;
double **x = atom->x;
int *type = atom->type;
tagint *tag = atom->tag;
memory->destroy(dnn);
memory->destroy(vect);
@ -479,9 +481,63 @@ void PairILPTMD::calc_normal()
for (m = 0; m < Nnei; m++) { dnormal[i][id][m][ip] = 0.0; }
}
}
// for hydrogen in water molecule
if (strcmp(elements[itype], "Hw") == 0) {
if(cont == 0) {
jH1 = atom->map(tag[i] - 1);
jH2 = atom->map(tag[i] - 2);
iH1 = map[type[jH1]];
iH2 = map[type[jH2]];
if (strcmp(elements[iH1], "Ow") == 0 ) {
vect[0][0] = x[jH1][0] - xtp;
vect[0][1] = x[jH1][1] - ytp;
vect[0][2] = x[jH1][2] - ztp;
} else if (strcmp(elements[iH2], "Ow") == 0 ) {
vect[0][0] = x[jH2][0] - xtp;
vect[0][1] = x[jH2][1] - ytp;
vect[0][2] = x[jH2][2] - ztp;
} else {
fprintf(screen, "jH1 jH2 = %d %d\n", jH1,jH2);
fprintf(screen, "Atom Type = %d %d\n", type[jH1],type[jH2]);
fprintf(screen, "For atom i = %d %d\n", tag[i],type[i]);
error->one(FLERR, "The order of atoms in water molecule should be O H H !");
}
}
Nave[0] = vect[0][0];
Nave[1] = vect[0][1];
Nave[2] = vect[0][2];
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// Calculte dNave/dri, defined as dpvdri
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
if (ip == id) { dpvdri[id][ip] = -1.0;}
else {dpvdri[id][ip] = 0.0;}
}
}
// derivatives of nn, dnn:3x1 vector
dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn;
dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn;
dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn;
// derivatives of unit vector ni respect to ri, the result is 3x3 matrix
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2;
dnormal[i][id][0][ip] = -dnormdri[i][id][ip];
}
}
}
}
//############################ For the edge atoms of TMD ################################
else if (cont < Nnei) {
else if (cont > 1 && cont < Nnei) {
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0) {
// derivatives of Ni[l] respect to the cont neighbors
@ -557,8 +613,7 @@ void PairILPTMD::calc_normal()
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
Nave[2] * dNave[2][m][id]) / nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
@ -589,12 +644,107 @@ void PairILPTMD::calc_normal()
}
}
} // for TMD
//############################ For Oxygen in the water molecule #######################
else if (strcmp(elements[itype], "Ow") == 0) {
if(cont == 0) {
jH1 = atom->map(tag[i] + 1);
jH2 = atom->map(tag[i] + 2);
iH1 = map[type[jH1]];
iH2 = map[type[jH2]];
if (strcmp(elements[iH1], "Hw") == 0 && strcmp(elements[iH2], "Hw") == 0) {
vect[0][0] = x[jH1][0] - xtp;
vect[0][1] = x[jH1][1] - ytp;
vect[0][2] = x[jH1][2] - ztp;
vect[1][0] = x[jH2][0] - xtp;
vect[1][1] = x[jH2][1] - ytp;
vect[1][2] = x[jH2][2] - ztp;
cont = 2;
} else {
fprintf(screen, "jH1 jH2 = %d %d\n", jH1,jH2);
fprintf(screen, "Atom Type = %d %d\n", type[jH1],type[jH2]);
fprintf(screen, "ID and type of atom i = %d %d\n", tag[i],type[i]);
error->one(FLERR, "The order of atoms in water molecule should be O H H !");
}
}
if (cont == 2) {
Nave[0] = (vect[0][0] + vect[1][0])/cont;
Nave[1] = (vect[0][1] + vect[1][1])/cont;
Nave[2] = (vect[0][2] + vect[1][2])/cont;
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// derivatives of non-normalized normal vector, dNave:3xcontx3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave
// respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (m = 0; m < cont; m++) {
if (ip == id) { dNave[id][m][ip] = 0.5;}
else {dNave[id][m][ip] = 0.0;}
}
}
}
// derivatives of nn, dnn:contx3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) / nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
// Calculte dNave/dri, defined as dpvdri
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dpvdri[id][ip] = 0.0;
for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; }
}
}
// derivatives of nn, dnn:3x1 vector
dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn;
dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn;
dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn;
// derivatives of unit vector ni respect to ri, the result is 3x3 matrix
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2;
}
}
// printf("%4d:\t%e %e %e\n\t%e %e %e\n\t%e %e %e\n",
// i, dnormdri[i][0][0], dnormdri[i][0][1], dnormdri[i][0][2],
// dnormdri[i][1][0],dnormdri[i][1][1],dnormdri[i][1][2],
// dnormdri[i][2][0],dnormdri[i][2][1],dnormdri[i][2][2]);
// exit(0);
}
else if (cont >= 3) {
error->one(FLERR,
"There are too many neighbors for calculating normals of water molecules");
}
}
//############################ For the edge & bulk atoms of GrhBN ################################
else {
if (cont == 2) {
for (ip = 0; ip < 3; ip++) {
pvet[0][ip] = vect[0][modulo(ip + 1, 3)] * vect[1][modulo(ip + 2, 3)] -
vect[0][modulo(ip + 2, 3)] * vect[1][modulo(ip + 1, 3)];
vect[0][modulo(ip + 2, 3)] * vect[1][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
@ -657,8 +807,7 @@ void PairILPTMD::calc_normal()
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
Nave[2] * dNave[2][m][id]) / nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.

2
src/INTERLAYER/pair_ilp_tmd.h
View File

@ -1,7 +1,7 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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

80
src/INTERLAYER/pair_ilp_water_2dm.cpp Normal file
View File

@ -0,0 +1,80 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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 author: Wengen Ouyang (Wuhan University)
e-mail: w.g.ouyang at gmail dot com
This is a full version of the potential described in
[Feng and Ouyang et al, J. Phys. Chem. C 127, 8704-8713 (2023).]
------------------------------------------------------------------------- */
#include "pair_ilp_water_2dm.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "interlayer_taper.h"
#include "memory.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace InterLayer;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
static const char cite_ilp_water[] = "ilp/water/2dm potential doi/10.1021/acs.jpcc.2c08464\n"
"@Article{Feng2023\n"
" author = {Z. Feng, Y. Yao, J. Liu, B. Wu, Z. Liu, and W. Ouyang},\n"
" title = {Registry-Dependent Potential for Interfaces of Water with Graphene},\n"
" journal = {J. Phys. Chem. C},\n"
" volume = 127,\n"
" pages = {8704-8713}\n"
" year = 2023,\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairILPWATER2DM::PairILPWATER2DM(LAMMPS *lmp) : PairILPGrapheneHBN(lmp), PairILPTMD(lmp)
{
variant = ILP_WATER_2DM;
single_enable = 0;
// for TMD, each atom have six neighbors
Nnei = 6;
if (lmp->citeme) lmp->citeme->add(cite_ilp_water);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairILPWATER2DM::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command");
if (!utils::strmatch(force->pair_style, "^hybrid/overlay"))
error->all(FLERR, "Pair style ilp/water/2dm must be used as sub-style with hybrid/overlay");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp);
}

67
src/INTERLAYER/pair_ilp_water_2dm.h Normal file
View File

@ -0,0 +1,67 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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
// clang-format off
PairStyle(ilp/water/2dm,PairILPWATER2DM);
// clang-format on
#else
#ifndef LMP_PAIR_ILP_WATER_2DM_H
#define LMP_PAIR_ILP_WATER_2DM_H
#include "pair_ilp_tmd.h"
namespace LAMMPS_NS {
class PairILPWATER2DM : virtual public PairILPTMD {
public:
PairILPWATER2DM(class LAMMPS *);
protected:
void settings(int, char **) override;
/**************************************************************/
/* modulo operation with cycling around range */
inline int modulo(int k, int range)
{
if (k < 0) k += range;
return k % range;
}
/**************************************************************/
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/