From 81a497adcd6d1fa52ac3bf0931bec76a7c8caaf4 Mon Sep 17 00:00:00 2001 From: oywg11 Date: Thu, 18 May 2023 21:51:53 +0800 Subject: [PATCH] add standard version of ilp_water_2dm --- src/INTERLAYER/pair_ilp_graphene_hbn.cpp | 25 +-- src/INTERLAYER/pair_ilp_graphene_hbn.h | 38 +++- src/INTERLAYER/pair_ilp_tmd.cpp | 213 +++++++++++++++++++---- src/INTERLAYER/pair_ilp_tmd.h | 2 +- src/INTERLAYER/pair_ilp_water_2dm.cpp | 80 +++++++++ src/INTERLAYER/pair_ilp_water_2dm.h | 67 +++++++ 6 files changed, 373 insertions(+), 52 deletions(-) create mode 100644 src/INTERLAYER/pair_ilp_water_2dm.cpp create mode 100644 src/INTERLAYER/pair_ilp_water_2dm.h diff --git a/src/INTERLAYER/pair_ilp_graphene_hbn.cpp b/src/INTERLAYER/pair_ilp_graphene_hbn.cpp index 0ff2811339..956d7ba2be 100644 --- a/src/INTERLAYER/pair_ilp_graphene_hbn.cpp +++ b/src/INTERLAYER/pair_ilp_graphene_hbn.cpp @@ -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 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; diff --git a/src/INTERLAYER/pair_ilp_graphene_hbn.h b/src/INTERLAYER/pair_ilp_graphene_hbn.h index 9987830b1d..6381fb8f35 100644 --- a/src/INTERLAYER/pair_ilp_graphene_hbn.h +++ b/src/INTERLAYER/pair_ilp_graphene_hbn.h @@ -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. + +*/ diff --git a/src/INTERLAYER/pair_ilp_tmd.cpp b/src/INTERLAYER/pair_ilp_tmd.cpp index c024e23079..d855b9423d 100644 --- a/src/INTERLAYER/pair_ilp_tmd.cpp +++ b/src/INTERLAYER/pair_ilp_tmd.cpp @@ -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 #include @@ -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. diff --git a/src/INTERLAYER/pair_ilp_tmd.h b/src/INTERLAYER/pair_ilp_tmd.h index 8381c2e830..8d5446d49c 100644 --- a/src/INTERLAYER/pair_ilp_tmd.h +++ b/src/INTERLAYER/pair_ilp_tmd.h @@ -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 diff --git a/src/INTERLAYER/pair_ilp_water_2dm.cpp b/src/INTERLAYER/pair_ilp_water_2dm.cpp new file mode 100644 index 0000000000..58e56894f5 --- /dev/null +++ b/src/INTERLAYER/pair_ilp_water_2dm.cpp @@ -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 +#include + +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); +} diff --git a/src/INTERLAYER/pair_ilp_water_2dm.h b/src/INTERLAYER/pair_ilp_water_2dm.h new file mode 100644 index 0000000000..c35d757ea0 --- /dev/null +++ b/src/INTERLAYER/pair_ilp_water_2dm.h @@ -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. + +*/