tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

CTIP pair style and qeq fix implemented

Browse Source
This commit is contained in:
Gabriel Plummer
2024-05-31 13:29:49 -07:00
parent e8ee0d9e71
commit e2e17b1326
7 changed files with 1072 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/QEQ/Install.sh
View File

@ -44,3 +44,5 @@ action fix_qeq_shielded.cpp
action fix_qeq_shielded.h
action fix_qeq_slater.cpp
action fix_qeq_slater.h
action fix_qeq_ctip.cpp
action fix_qeq_ctip.h

43
src/QEQ/fix_qeq.cpp
View File

@ -27,7 +27,9 @@
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "pair.h"
#include "respa.h"
#include "suffix.h"
#include "text_file_reader.h"
#include "update.h"
@ -53,12 +55,12 @@ namespace {
FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), list(nullptr), chi(nullptr), eta(nullptr),
gamma(nullptr), zeta(nullptr), zcore(nullptr), chizj(nullptr), shld(nullptr),
gamma(nullptr), zeta(nullptr), zcore(nullptr), qmin(nullptr), qmax(nullptr), omega(nullptr), chizj(nullptr), shld(nullptr),
s(nullptr), t(nullptr), s_hist(nullptr), t_hist(nullptr), Hdia_inv(nullptr), b_s(nullptr),
b_t(nullptr), p(nullptr), q(nullptr), r(nullptr), d(nullptr),
qf(nullptr), q1(nullptr), q2(nullptr), qv(nullptr)
{
if (narg < 8) utils::missing_cmd_args(FLERR, "fix " + std::string(style), error);
if (narg < 8) error->all(FLERR,"Illegal fix qeq command");
scalar_flag = 1;
extscalar = 0;
@ -74,9 +76,6 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
if ((nevery <= 0) || (cutoff <= 0.0) || (tolerance <= 0.0) || (maxiter <= 0))
error->all(FLERR,"Illegal fix qeq command");
// must have charges
if (!atom->q_flag) error->all(FLERR, "Fix {} requires atom attribute q", style);
alpha = 0.20;
swa = 0.0;
swb = cutoff;
@ -115,6 +114,7 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
q2 = nullptr;
streitz_flag = 0;
reax_flag = 0;
ctip_flag = 0;
qv = nullptr;
comm_forward = comm_reverse = 1;
@ -134,6 +134,8 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
streitz_flag = 1;
} else if (utils::strmatch(arg[7],"^reax..")) {
reax_flag = 1;
} else if (strcmp(arg[7],"coul/ctip") == 0) {
ctip_flag = 1;
} else {
read_file(arg[7]);
}
@ -144,8 +146,7 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
FixQEq::~FixQEq()
{
// unregister callbacks to this fix from Atom class
if (modify->get_fix_by_id(id)) atom->delete_callback(id,Atom::GROW);
atom->delete_callback(id,Atom::GROW);
memory->destroy(s_hist);
memory->destroy(t_hist);
@ -161,6 +162,9 @@ FixQEq::~FixQEq()
memory->destroy(gamma);
memory->destroy(zeta);
memory->destroy(zcore);
memory->destroy(qmin);
memory->destroy(qmax);
memory->destroy(omega);
}
}
@ -339,6 +343,12 @@ void FixQEq::setup_pre_force(int vflag)
if (force->newton_pair == 0)
error->all(FLERR,"QEQ with 'newton pair off' not supported");
if (force->pair) {
if (force->pair->suffix_flag & (Suffix::INTEL|Suffix::GPU))
error->all(FLERR,"QEQ is not compatiple with suffix version "
"of pair style");
}
deallocate_storage();
allocate_storage();
@ -752,6 +762,9 @@ void FixQEq::read_file(char *file)
memory->create(gamma,ntypes+1,"qeq:gamma");
memory->create(zeta,ntypes+1,"qeq:zeta");
memory->create(zcore,ntypes+1,"qeq:zcore");
memory->create(qmin,ntypes+1,"qeq:qmin");
memory->create(qmax,ntypes+1,"qeq:qmax");
memory->create(omega,ntypes+1,"qeq:omega");
// read each line out of file, skipping blank lines or leading '#'
// store line of params if all 3 element tags are in element list
@ -760,7 +773,7 @@ void FixQEq::read_file(char *file)
int *setflag = new int[ntypes+1];
for (int n=0; n <= ntypes; ++n) {
setflag[n] = 0;
chi[n] = eta[n] = gamma[n] = zeta[n] = zcore[n] = 0.0;
chi[n] = eta[n] = gamma[n] = zeta[n] = zcore[n] = qmin[n] = qmax[n] = omega[n] = 0.0;
}
try {
@ -777,7 +790,7 @@ void FixQEq::read_file(char *file)
auto values = reader.next_values(0);
if (values.count() == 0) continue;
if (values.count() < 6)
if (values.count() < 9)
throw qeq_parser_error("Invalid qeq parameter file");
auto word = values.next_string();
@ -795,6 +808,12 @@ void FixQEq::read_file(char *file)
for (int n=nlo; n <= nhi; ++n) zeta[n] = val;
val = values.next_double();
for (int n=nlo; n <= nhi; ++n) zcore[n] = val;
val = values.next_double();
for (int n=nlo; n <= nhi; ++n) qmin[n] = val;
val = values.next_double();
for (int n=nlo; n <= nhi; ++n) qmax[n] = val;
val = values.next_double();
for (int n=nlo; n <= nhi; ++n) omega[n] = val;
for (int n=nlo; n <= nhi; ++n) setflag[n] = 1;
}
} catch (EOFException &) {
@ -805,7 +824,8 @@ void FixQEq::read_file(char *file)
for (int n=1; n <= ntypes; ++n)
if (setflag[n] == 0)
error->one(FLERR,"Parameters for atom type {} missing in qeq parameter file", n);
error->one(FLERR,fmt::format("Parameters for atom type {} missing in "
"qeq parameter file", n));
delete[] setflag;
}
@ -814,4 +834,7 @@ void FixQEq::read_file(char *file)
MPI_Bcast(gamma,ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(zeta,ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(zcore,ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(qmin,ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(qmax,ntypes+1,MPI_DOUBLE,0,world);
MPI_Bcast(omega,ntypes+1,MPI_DOUBLE,0,world);
}

9
src/QEQ/fix_qeq.h
View File

@ -70,10 +70,10 @@ class FixQEq : public Fix {
int maxwarn; // print warning when max iterations was reached
double cutoff, cutoff_sq; // neighbor cutoff
double *chi, *eta, *gamma, *zeta, *zcore; // qeq parameters
double *chi, *eta, *gamma, *zeta, *zcore, *qmin, *qmax, *omega;
double *chizj;
double **shld;
int streitz_flag, reax_flag;
int streitz_flag, reax_flag, ctip_flag;
bigint ngroup;
@ -100,6 +100,11 @@ class FixQEq : public Fix {
double alpha;
// ctip
double cdamp;
int maxrepeat;
// damped dynamics
double *qf, *q1, *q2, qdamp, qstep;

394
src/QEQ/fix_qeq_ctip.cpp Normal file
View File

@ -0,0 +1,394 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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: Gabriel Plummer (NASA)
------------------------------------------------------------------------- */
#include "fix_qeq_ctip.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kspace.h"
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "pair.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace FixConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
FixQEqCTIP::FixQEqCTIP(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg) {
cdamp = 0.30;
maxrepeat = 10;
// optional arg
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg], "cdamp") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/ctip cdamp", error);
cdamp = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "maxrepeat") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/ctip maxrepeat", error);
maxrepeat = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "warn") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/ctip warn", error);
maxwarn = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else
error->all(FLERR, "Unknown fix qeq/ctip keyword: {}", arg[iarg]);
}
if (ctip_flag) extract_ctip();
}
/* ---------------------------------------------------------------------- */
void FixQEqCTIP::init()
{
FixQEq::init();
neighbor->add_request(this, NeighConst::REQ_FULL);
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shielding");
}
/* ---------------------------------------------------------------------- */
void FixQEqCTIP::extract_ctip()
{
Pair *pair = force->pair_match("coul/ctip",1);
if (pair == nullptr) error->all(FLERR,"No pair coul/ctip for fix qeq/ctip");
int tmp;
chi = (double *) pair->extract("chi",tmp);
eta = (double *) pair->extract("eta",tmp);
gamma = (double *) pair->extract("gamma",tmp);
zeta = (double *) pair->extract("zeta",tmp);
zcore = (double *) pair->extract("zcore",tmp);
qmin = (double *) pair->extract("qmin",tmp);
qmax = (double *) pair->extract("qmax",tmp);
omega = (double *) pair->extract("omega",tmp);
if (chi == nullptr || eta == nullptr || gamma == nullptr
|| zeta == nullptr || zcore == nullptr || qmin == nullptr || qmax == nullptr || omega == nullptr)
error->all(FLERR,
"Fix qeq/ctip could not extract params from pair coul/ctip");
}
/* ---------------------------------------------------------------------- */
void FixQEqCTIP::pre_force(int /*vflag*/)
{
int i,n;
if (update->ntimestep % nevery) return;
nlocal = atom->nlocal;
if (atom->nmax > nmax) reallocate_storage();
if (nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
reallocate_matrix();
for (i=1; i <= maxrepeat; i++) {
init_matvec();
matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel
matvecs /= 2;
n=calculate_check_Q();
MPI_Allreduce(&n, &nout, 1, MPI_INTEGER, MPI_SUM, world);
if (nout == 0) break;
}
if (i > maxrepeat && comm->me == 0)
error->all(FLERR,"Fix qeq some charges not bound within the domain");
if (force->kspace) force->kspace->qsum_qsq();
}
/* ---------------------------------------------------------------------- */
void FixQEqCTIP::init_matvec()
{
compute_H();
int inum, ii, i;
int *ilist;
double *q = atom->q, qi;
double r = cutoff;
double rsq = r*r;
double r6 = rsq*rsq*rsq;
double erfcd_cut = exp(-cdamp * cdamp * rsq);
double t_cut = 1.0 / (1.0 + EWALD_P * cdamp * r);
double erfcc_cut = (t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut) / r;
int elt1;
double reff[atom->ntypes], reffsq[atom->ntypes], reff4[atom->ntypes], reff7[atom->ntypes], s2d_self[atom->ntypes];
for (elt1 = 0; elt1 < atom->ntypes; elt1++) {
reff[elt1] = cbrt(rsq*r + 1/(gamma[elt1+1]*gamma[elt1+1]*gamma[elt1+1]));
reffsq[elt1] = reff[elt1]*reff[elt1];
reff4[elt1] = reffsq[elt1]*reffsq[elt1];
reff7[elt1] = reff4[elt1]*reffsq[elt1]*reff[elt1];
s2d_self[elt1] = 2.0*force->qqr2e*(1.5*erfcc_cut + 2.0*cdamp/MY_PIS*erfcd_cut + cdamp*cdamp*cdamp/MY_PIS*rsq*erfcd_cut + 0.5/reff[elt1] - 1.5/r + r6/reff7[elt1] + cdamp/MY_PIS);
}
inum = list->inum;
ilist = list->ilist;
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
qi=q[i];
if (qi < qmin[atom->type[i]]) {
Hdia_inv[i] = 1. / (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1]);
b_s[i] = -((chi[atom->type[i]]-2*qmin[atom->type[i]]*omega[atom->type[i]]) + chizj[i]);
} else if (qi < qmax[atom->type[i]]) {
Hdia_inv[i] = 1. / (eta[atom->type[i]]-s2d_self[atom->type[i]-1]);
b_s[i] = -(chi[atom->type[i]] + chizj[i]);
} else {
Hdia_inv[i] = 1. / (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1]);
b_s[i] = -((chi[atom->type[i]]-2*qmax[atom->type[i]]*omega[atom->type[i]]) + chizj[i]);
}
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]);
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}
pack_flag = 2;
comm->forward_comm(this); //Dist_vector(s);
pack_flag = 3;
comm->forward_comm(this); //Dist_vector(t);
}
/* ---------------------------------------------------------------------- */
void FixQEqCTIP::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int i, j, ii, jj;
double **x;
double dx, dy, dz, r_sqr, r, reff;
double cutoffsq, cutoffcu, cutoff4, cdampcu, erfcc_cut, erfcd_cut, t_cut;
double erfcc, erfcd, t;
int elt1, elt2;
double shield[atom->ntypes][atom->ntypes], shieldcu[atom->ntypes][atom->ntypes], reffc[atom->ntypes][atom->ntypes], reffcsq[atom->ntypes][atom->ntypes], reffc4[atom->ntypes][atom->ntypes], reffc7[atom->ntypes][atom->ntypes];
double s2d_shift[atom->ntypes][atom->ntypes], f_shift[atom->ntypes][atom->ntypes], e_shift[atom->ntypes][atom->ntypes];
x = atom->x;
int *mask = atom->mask;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
cutoffsq = cutoff * cutoff;
cutoffcu = cutoffsq * cutoff;
cutoff4 = cutoffsq * cutoffsq;
cdampcu = cdamp * cdamp * cdamp;
erfcd_cut = exp(-cdamp * cdamp * cutoffsq);
t_cut = 1.0 / (1.0 + EWALD_P * cdamp * cutoff);
erfcc_cut = t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut;
for (elt1 = 0; elt1 < atom->ntypes; elt1++) {
for (elt2 = 0; elt2 < atom->ntypes; elt2++) {
shield[elt1][elt2] = sqrt(gamma[elt1+1] * gamma[elt2+1]);
shieldcu[elt1][elt2] = shield[elt1][elt2] * shield[elt1][elt2] * shield[elt1][elt2];
reffc[elt1][elt2] = cbrt(cutoffcu + 1 / shieldcu[elt1][elt2]);
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cutoffcu + 4.0 * cdamp / MY_PIS * erfcd_cut / cutoffsq + 4.0 * cdampcu / MY_PIS * erfcd_cut - 2 / cutoffcu + 4 * cutoff4 / reffc7[elt1][elt2] - 2 * cutoff / reffc4[elt1][elt2];
f_shift[elt1][elt2] = erfcc_cut / cutoffsq + 2.0 * cdamp / MY_PIS * erfcd_cut / cutoff - 1 / cutoffsq + cutoffsq / reffc4[elt1][elt2];
e_shift[elt1][elt2] = erfcc_cut / cutoff + 1 / reffc[elt1][elt2] - 1 / cutoff;
}
}
// fill in the H matrix
m_fill = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
jnum = numneigh[i];
H.firstnbr[i] = m_fill;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
dz = x[j][2] - x[i][2];
r_sqr = dx*dx + dy*dy + dz*dz;
if (r_sqr <= cutoff_sq) {
H.jlist[m_fill] = j;
r = sqrt(r_sqr);
reff = cbrt(r_sqr * r + 1 / shieldcu[atom->type[i]-1][atom->type[j]-1]);
erfcd = exp(-cdamp * cdamp * r_sqr);
t = 1.0 / (1.0 + EWALD_P * cdamp * r);
erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
H.val[m_fill] = 0.5*force->qqr2e*(erfcc/r+1/reff-1/r-e_shift[atom->type[i]-1][atom->type[j]-1]+f_shift[atom->type[i]-1][atom->type[j]-1]*(r-cutoff)-s2d_shift[atom->type[i]-1][atom->type[j]-1]*0.5*(r-cutoff)*(r-cutoff));
m_fill++;
}
}
H.numnbrs[i] = m_fill - H.firstnbr[i];
}
}
if (m_fill >= H.m)
error->all(FLERR,"Fix qeq/ctip has insufficient H matrix "
"size: m_fill={} H.m={}\n",m_fill, H.m);
}
/* ---------------------------------------------------------------------- */
void FixQEqCTIP::sparse_matvec(sparse_matrix *A, double *x, double *b)
{
int i, j, itr_j;
double *q=atom->q, qi;
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
double r = cutoff;
double rsq = r*r;
double r6 = rsq*rsq*rsq;
double erfcd_cut = exp(-cdamp * cdamp * rsq);
double t_cut = 1.0 / (1.0 + EWALD_P * cdamp * r);
double erfcc_cut = (t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut) / r;
int elt1;
double reff[atom->ntypes], reffsq[atom->ntypes], reff4[atom->ntypes], reff7[atom->ntypes], s2d_self[atom->ntypes];
for (elt1 = 0; elt1 < atom->ntypes; elt1++) {
reff[elt1] = cbrt(rsq*r + 1/(gamma[elt1+1]*gamma[elt1+1]*gamma[elt1+1]));
reffsq[elt1] = reff[elt1]*reff[elt1];
reff4[elt1] = reffsq[elt1]*reffsq[elt1];
reff7[elt1] = reff4[elt1]*reffsq[elt1]*reff[elt1];
s2d_self[elt1] = 2.0*force->qqr2e*(1.5*erfcc_cut + 2.0*cdamp/MY_PIS*erfcd_cut + cdamp*cdamp*cdamp/MY_PIS*rsq*erfcd_cut + 0.5/reff[elt1] - 1.5/r + r6/reff7[elt1] + cdamp/MY_PIS);
}
for (i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit) {
qi=q[i];
if (qi < qmin[atom->type[i]]) {
b[i] = (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1])*x[i];
} else if (qi < qmax[atom->type[i]]) {
b[i] = (eta[atom->type[i]]-s2d_self[atom->type[i]-1]) * x[i];
} else {
b[i] = (eta[atom->type[i]]+2*omega[atom->type[i]]-s2d_self[atom->type[i]-1])*x[i];
}
}
}
for (i = nlocal; i < nall; ++i) {
if (atom->mask[i] & groupbit)
b[i] = 0;
}
for (i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit) {
for (itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
j = A->jlist[itr_j];
b[i] += A->val[itr_j] * x[j];
b[j] += A->val[itr_j] * x[i];
}
}
}
}
/* ---------------------------------------------------------------------- */
int FixQEqCTIP::calculate_check_Q()
{
int i, k, inum, ii;
int *ilist;
double u, s_sum, t_sum;
double *q = atom->q;
double qi_old,qi_new;
double qi_check1,qi_check2;
double qi_check3;
int n;
inum = list->inum;
ilist = list->ilist;
s_sum = parallel_vector_acc( s, inum );
t_sum = parallel_vector_acc( t, inum);
u = s_sum / t_sum;
n = 0;
for( ii = 0; ii < inum; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
qi_old = q[i];
q[i] = s[i] - u * t[i];
for( k = 4; k > 0; --k ) {
s_hist[i][k] = s_hist[i][k-1];
t_hist[i][k] = t_hist[i][k-1];
}
s_hist[i][0] = s[i];
t_hist[i][0] = t[i];
qi_new = q[i];
qi_check1=(qi_new-qmin[atom->type[i]])*(qi_old-qmin[atom->type[i]]);
qi_check2=(qi_new-qmax[atom->type[i]])*(qi_old-qmax[atom->type[i]]);
if ( qi_check1 < 0.0 || qi_check2 < 0.0 ) {
qi_check3=abs(qi_new-qi_old);
if (qi_check3 > tolerance) n++;
}
}
}
pack_flag = 4;
comm->forward_comm( this ); //Dist_vector( atom->q );
return n;
}

44
src/QEQ/fix_qeq_ctip.h Normal file
View File

@ -0,0 +1,44 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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 FIX_CLASS
// clang-format off
FixStyle(qeq/ctip,FixQEqCTIP);
// clang-format on
#else
#ifndef LMP_FIX_QEQ_CTIP_H
#define LMP_FIX_QEQ_CTIP_H
#include "fix_qeq.h"
namespace LAMMPS_NS {
class FixQEqCTIP : public FixQEq {
public:
FixQEqCTIP(class LAMMPS *, int, char **);
void init() override;
void pre_force(int) override;
private:
void init_matvec();
void sparse_matvec(sparse_matrix *, double *, double *) override;
void compute_H();
void extract_ctip();
int calculate_check_Q();
int nout;
};
} // namespace LAMMPS_NS
#endif
#endif

525
src/pair_coul_ctip.cpp Normal file
View File

@ -0,0 +1,525 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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: Gabriel Plummer (NASA)
------------------------------------------------------------------------- */
#include "pair_coul_ctip.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairCoulCTIP::PairCoulCTIP(LAMMPS *lmp) : Pair(lmp)
{
params= nullptr;
}
/* ---------------------------------------------------------------------- */
PairCoulCTIP::~PairCoulCTIP()
{
if (copymode) return;
memory->sfree(params);
memory->destroy(elem1param);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(qeq_x);
memory->destroy(qeq_j);
memory->destroy(qeq_g);
memory->destroy(qeq_z);
memory->destroy(qeq_c);
memory->destroy(qeq_q1);
memory->destroy(qeq_q2);
memory->destroy(qeq_w);
}
}
/* ---------------------------------------------------------------------- */
void PairCoulCTIP::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
int iparam_i, jparam_j;
double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair;
double r, rsq, reff, reffsq, reff4, forcecoul, factor_coul;
double prefactor, erfcc, erfcd, t;
double selfion;
int *ilist, *jlist, *numneigh, **firstneigh;
double cut_coulcu, cut_coul4, alphacu, erfcd_cut, t_cut, erfcc_cut;
int elt1, elt2;
double shield[nelements][nelements], shieldcu[nelements][nelements], reffc[nelements][nelements], reffcsq[nelements][nelements], reffc4[nelements][nelements], reffc7[nelements][nelements];
double s2d_shift[nelements][nelements], f_shift[nelements][nelements], e_shift[nelements][nelements], self_factor[nelements][nelements];
ecoul = 0.0;
selfion = 0.0;
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
cut_coulcu = cut_coulsq * cut_coul;
cut_coul4 = cut_coulsq * cut_coulsq;
alphacu = alpha * alpha * alpha;
erfcd_cut = exp(-alpha * alpha * cut_coulsq);
t_cut = 1.0 / (1.0 + EWALD_P * alpha * cut_coul);
erfcc_cut = t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut;
for (elt1 = 0; elt1 < nelements; elt1++) {
for (elt2 = 0; elt2 < nelements; elt2++) {
shield[elt1][elt2] = sqrt(params[elt1].gamma * params[elt2].gamma);
shieldcu[elt1][elt2] = shield[elt1][elt2] * shield[elt1][elt2] * shield[elt1][elt2];
reffc[elt1][elt2] = cbrt(cut_coulcu + 1 / shieldcu[elt1][elt2]);
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coulcu + 4.0 * alpha / MY_PIS * erfcd_cut / cut_coulsq + 4.0 * alphacu / MY_PIS * erfcd_cut - 2 / cut_coulcu + 4 * cut_coul4 / reffc7[elt1][elt2] - 2 * cut_coul / reffc4[elt1][elt2];
f_shift[elt1][elt2] = erfcc_cut / cut_coulsq + 2.0 * alpha / MY_PIS * erfcd_cut / cut_coul - 1 / cut_coulsq + cut_coulsq / reffc4[elt1][elt2];
e_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coul + 2.0 * alpha / MY_PIS * erfcd_cut - 2 / cut_coul + 1 / reffc[elt1][elt2] + cut_coulcu / reffc4[elt1][elt2] + s2d_shift[elt1][elt2] * cut_coulsq * 0.5;
self_factor[elt1][elt2] = -(e_shift[elt1][elt2] * 0.5 + alpha / MY_PIS) * qqrd2e;
}
}
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = map[type[i]];
jlist = firstneigh[i];
jnum = numneigh[i];
iparam_i = elem1param[itype];
selfion = self(&params[iparam_i],qtmp);
if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
if (eflag) {
double e_self = self_factor[iparam_i][iparam_i] * qtmp * qtmp;
ev_tally(i, i, nlocal, 0, 0.0, e_self, 0.0, 0.0, 0.0, 0.0);
}
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_coulsq) {
jtype = map[type[j]];
jparam_j = elem1param[jtype];
r = sqrt(rsq);
reff = cbrt(rsq * r + 1/shieldcu[iparam_i][jparam_j]);
reffsq = reff * reff;
reff4 = reffsq * reffsq;
prefactor = qqrd2e * qtmp * q[j] / r;
erfcd = exp(-alpha * alpha * rsq);
t = 1.0 / (1.0 + EWALD_P * alpha * r);
erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
forcecoul = prefactor * (erfcc / r + 2.0 * alpha / MY_PIS * erfcd + rsq * r / reff4 - 1 / r - r * f_shift[iparam_i][jparam_j] + r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul)) * r;
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
fpair = forcecoul / rsq;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) {
ecoul = prefactor * (erfcc + r / reff - 1 - r * erfcc_cut / cut_coul - r / reffc[iparam_i][jparam_j] + r / cut_coul + r * f_shift[iparam_i][jparam_j] * (r - cut_coul) - r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul) * (r - cut_coul)* 0.5);
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
} else
ecoul = 0.0;
if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
double PairCoulCTIP::self(Param *param, double qi)
{
double s1=param->chi, s2=param->eta, s3=param->qmin, s4=param->qmax, s5=param->omega;
if ( qi < s3 ) {
return qi*((s1-2*s3*s5)+qi*(0.50*s2+s5))+s3*s3*s5;
} else if ( qi < s4 ) {
return qi*(s1+qi*(0.50*s2));
} else {
return qi*((s1-2*s4*s5)+qi*(0.50*s2+s5))+s4*s4*s5;
}
return 0.0;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairCoulCTIP::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag, n + 1, n + 1, "pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++) setflag[i][j] = 0;
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(qeq_x,n+1,"pair:qeq_x");
memory->create(qeq_j,n+1,"pair:qeq_j");
memory->create(qeq_g,n+1,"pair:qeq_g");
memory->create(qeq_z,n+1,"pair:qeq_z");
memory->create(qeq_c,n+1,"pair:qeq_c");
memory->create(qeq_q1,n+1,"pair:qeq_q1");
memory->create(qeq_q2,n+1,"pair:qeq_q2");
memory->create(qeq_w,n+1,"pair:qeq_w");
map = new int[n+1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairCoulCTIP::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR, "Illegal pair_style command");
alpha = utils::numeric(FLERR, arg[0], false, lmp);
cut_coul = utils::numeric(FLERR, arg[1], false, lmp);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairCoulCTIP::coeff(int narg, char **arg)
{
if (!allocated) allocate();
map_element2type(narg-3, arg+3);
// read potential file and initialize potential parameters
read_file(arg[2]);
setup_params();
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairCoulCTIP::init_style()
{
if (!atom->q_flag) error->all(FLERR, "Pair style coul/ctip requires atom attribute q");
neighbor->add_request(this);
cut_coulsq = cut_coul * cut_coul;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCoulCTIP::init_one(int /*i*/, int /*j*/)
{
return cut_coul;
}
/* ---------------------------------------------------------------------- */
void PairCoulCTIP::read_file(char *file)
{
memory->sfree(params);
params = nullptr;
nparams = 0;
maxparam = 0;
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "coul/ctip");
char * line;
while ((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
std::string iname = values.next_string();
// ielement = 1st args
int ielement;
for (ielement = 0; ielement < nelements; ielement++)
if (iname == elements[ielement]) break;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
// make certain all addional allocated storage is initialized
// to avoid false positives when checking with valgrind
memset(params + nparams, 0, DELTA*sizeof(Param));
}
params[nparams].ielement = ielement;
params[nparams].chi = values.next_double();
params[nparams].eta = values.next_double();
params[nparams].gamma = values.next_double();
params[nparams].zeta = values.next_double();
params[nparams].zcore = values.next_double();
params[nparams].qmin = values.next_double();
params[nparams].qmax = values.next_double();
params[nparams].omega = values.next_double();
} catch (TokenizerException &e) {
error->one(FLERR, e.what());
}
// parameter sanity check
if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 ||
params[nparams].zcore < 0.0 || params[nparams].gamma < 0.0 || params[nparams].qmin > params[nparams].qmax || params[nparams].omega < 0.0 )
error->one(FLERR,"Illegal coul/ctip parameter");
nparams++;
}
}
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if (comm->me != 0) {
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
}
/* ---------------------------------------------------------------------- */
void PairCoulCTIP::setup_params()
{
int i,m,n;
// set elem1param
memory->destroy(elem1param);
memory->create(elem1param,nelements,"pair:elem1param");
for (i = 0; i < nelements; i++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
elem1param[i] = n;
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulCTIP::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j], sizeof(int), 1, fp); }
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulCTIP::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulCTIP::write_restart_settings(FILE *fp)
{
fwrite(&alpha, sizeof(double), 1, fp);
fwrite(&cut_coul, sizeof(double), 1, fp);
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulCTIP::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR, &alpha, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
}
MPI_Bcast(&alpha, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */
void *PairCoulCTIP::extract(const char *str, int &dim)
{
if (strcmp(str, "cut_coul") == 0) {
dim = 0;
return (void *) &cut_coul;
}
if (strcmp(str,"chi") == 0 && qeq_x) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_x[i] = params[map[i]].chi;
else qeq_x[i] = 0.0;
return (void *) qeq_x;
}
if (strcmp(str,"eta") == 0 && qeq_j) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_j[i] = params[map[i]].eta;
else qeq_j[i] = 0.0;
return (void *) qeq_j;
}
if (strcmp(str,"gamma") == 0 && qeq_g) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_g[i] = params[map[i]].gamma;
else qeq_g[i] = 0.0;
return (void *) qeq_g;
}
if (strcmp(str,"zeta") == 0 && qeq_z) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_z[i] = params[map[i]].zeta;
else qeq_z[i] = 0.0;
return (void *) qeq_z;
}
if (strcmp(str,"zcore") == 0 && qeq_c) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_c[i] = params[map[i]].zcore;
else qeq_c[i] = 0.0;
return (void *) qeq_c;
}
if (strcmp(str,"qmin") == 0 && qeq_q1) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_q1[i] = params[map[i]].qmin;
else qeq_q1[i] = 0.0;
return (void *) qeq_q1;
}
if (strcmp(str,"qmax") == 0 && qeq_q2) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_q2[i] = params[map[i]].qmax;
else qeq_q2[i] = 0.0;
return (void *) qeq_q2;
}
if (strcmp(str,"omega") == 0 && qeq_w) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_w[i] = params[map[i]].omega;
else qeq_w[i] = 0.0;
return (void *) qeq_w;
}
return nullptr;
}

67
src/pair_coul_ctip.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
LAMMPS development team: developers@lammps.org
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(coul/ctip,PairCoulCTIP);
// clang-format on
#else
#ifndef LMP_PAIR_COUL_CTIP_H
#define LMP_PAIR_COUL_CTIP_H
#include "pair.h"
namespace LAMMPS_NS {
class PairCoulCTIP : public Pair {
public:
PairCoulCTIP(class LAMMPS *);
~PairCoulCTIP() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void *extract(const char *, int &) override;
static constexpr int NPARAMS_PER_LINE = 9;
protected:
double cut_coul, cut_coulsq;
double alpha;
struct Param {
double chi, eta, gamma, zeta, zcore, qmin, qmax, omega;
int ielement;
};
Param *params;
double *qeq_x, *qeq_j, *qeq_g, *qeq_z, *qeq_c, *qeq_q1, *qeq_q2, *qeq_w;
void allocate();
void read_file(char *);
void setup_params();
double self(Param *, double);
};
} // namespace LAMMPS_NS
#endif
#endif