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

changed some comments

Browse Source
This commit is contained in:
Daniele
2021-11-11 16:12:26 +01:00
parent 4f4b18ab7e
commit 1c0e6578b6
4 changed files with 578 additions and 614 deletions
Download Patch File Download Diff File Expand all files Collapse all files

537
src/pair_smatb.cpp
View File

@ -1,58 +1,58 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This pair style is written by Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#include "pair_smatb.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_smatb.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
//#define DR_DEBUG
#ifdef DR_DEBUG
#include <iostream>
#include <iomanip>
#endif //DR_DEBUG
using namespace LAMMPS_NS;
#define MAXLINE 1024
PairSMATB::PairSMATB(LAMMPS *lmp)
: Pair(lmp),
nmax(0),
on_eb(nullptr),
r0(nullptr),
p(nullptr),
A(nullptr),
q(nullptr),
QSI(nullptr),
cutOffStart(nullptr),
cutOffEnd(nullptr),
cutOffEnd2(nullptr),
a3(nullptr), a4(nullptr), a5(nullptr),
x3(nullptr), x4(nullptr), x5(nullptr) {
single_enable = 0; // 1 if single() routine exists
restartinfo = 1; // 1 if pair style writes restart info
respa_enable = 0; // 1 if inner/middle/outer rRESPA routines
one_coeff = 0; // 1 if allows only one coeff * * call
manybody_flag = 1; // 1 if a manybody potential
no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
writedata = 1; // 1 if writes coeffs to data file
ghostneigh = 0; // 1 if pair style needs neighbors of ghosts
PairSMATB::PairSMATB(LAMMPS *lmp) :
Pair(lmp), nmax(0), on_eb(nullptr), r0(nullptr), p(nullptr), A(nullptr), q(nullptr),
QSI(nullptr), cutOffStart(nullptr), cutOffEnd(nullptr), cutOffEnd2(nullptr), a3(nullptr),
a4(nullptr), a5(nullptr), x3(nullptr), x4(nullptr), x5(nullptr)
{
single_enable = 0; // 1 if single() routine exists
restartinfo = 1; // 1 if pair style writes restart info
respa_enable = 0; // 1 if inner/middle/outer rRESPA routines
one_coeff = 0; // 1 if allows only one coeff * * call
manybody_flag = 1; // 1 if a manybody potential
no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
writedata = 1; // 1 if writes coeffs to data file
ghostneigh = 0; // 1 if pair style needs neighbors of ghosts
// set comm size needed by this Pair
comm_forward = 1;
comm_reverse = 1;
}
PairSMATB::~PairSMATB() {
if (copymode) {
return;
}
PairSMATB::~PairSMATB()
{
if (copymode) { return; }
memory->destroy(on_eb);
if (allocated) {
memory->destroy(setflag);
@ -75,17 +75,18 @@ PairSMATB::~PairSMATB() {
}
}
void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes pairwise interactions
void PairSMATB::compute(int eflag, int vflag)
{ //workhorse routine that computes pairwise interactions
//eflag means compute energy
//vflag means compute virial
int i,j,ii,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,del[3],fpair;
double dijsq,dij;
int i, j, ii, jj, jnum, itype, jtype;
double xtmp, ytmp, ztmp, del[3], fpair;
double dijsq, dij;
double espo, aexpp, qsiexpq, eb_i, Fb, Fr;
double polyval, polyval2, polyval3, polyval4, polyval5;
//sets up the flags for energy caclulations
if (eflag || vflag) {
ev_setup(eflag,vflag);
ev_setup(eflag, vflag);
eng_vdwl = 0;
} else {
evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
@ -96,14 +97,14 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
if (atom->nmax > nmax) {
nmax = atom->nmax;
memory->grow(on_eb,nmax,"pair_smatb:on_eb");
memory->grow(on_eb, nmax, "pair_smatb:on_eb");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair;
@ -122,39 +123,41 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
//FIRST LOOP: CALCULATES the squared bounding energy and accumulate it in on_eb for each atom
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; ++jj) {
j = jlist[jj];
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
del[0] = xtmp - x[j][0];
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if ( dijsq < cutOffEnd2[itype][jtype] ) {
dij = sqrt (dijsq);
if ( dij < cutOffStart[itype][jtype] ) {
qsiexpq = (QSI[itype][jtype]*QSI[itype][jtype]) * exp(2.0*q[itype][jtype]*(1.0 - dij/r0[itype][jtype]));
} else {
polyval = dij-cutOffEnd[itype][jtype];
polyval3 = polyval*polyval*polyval;
polyval4 = polyval3*polyval;
polyval5 = polyval4*polyval;
qsiexpq = x5[itype][jtype]*polyval5+x4[itype][jtype]*polyval4+x3[itype][jtype]*polyval3;
qsiexpq = qsiexpq* qsiexpq;
}
on_eb[i]+=qsiexpq;
//if (newton_pair || j < nlocal) {
on_eb[j]+=qsiexpq;
//}
if (dijsq < cutOffEnd2[itype][jtype]) {
dij = sqrt(dijsq);
if (dij < cutOffStart[itype][jtype]) {
qsiexpq = (QSI[itype][jtype] * QSI[itype][jtype]) *
exp(2.0 * q[itype][jtype] * (1.0 - dij / r0[itype][jtype]));
} else {
polyval = dij - cutOffEnd[itype][jtype];
polyval3 = polyval * polyval * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
qsiexpq = x5[itype][jtype] * polyval5 + x4[itype][jtype] * polyval4 +
x3[itype][jtype] * polyval3;
qsiexpq = qsiexpq * qsiexpq;
}
on_eb[i] += qsiexpq;
//if (newton_pair || j < nlocal) {
on_eb[j] += qsiexpq;
//}
}
}
}
@ -167,22 +170,18 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
// the store the reciprocal in on_eb in order to not do it in the SECOND LOOP
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (i < nlocal){
eb_i=sqrt(on_eb[i]);
if (eb_i!=0.0) {
on_eb[i]=1.0/eb_i;
i = ilist[ii];
if (i < nlocal) {
eb_i = sqrt(on_eb[i]);
if (eb_i != 0.0) {
on_eb[i] = 1.0 / eb_i;
} else {
on_eb[i] = 0.0;
}
//if needed the bounding energy is accumulated:
if (eflag_either) {
if (eflag_atom) {
eatom[i] -= eb_i;
}
if (eflag_global) {
eng_vdwl -= eb_i;
}
if (eflag_atom) { eatom[i] -= eb_i; }
if (eflag_global) { eng_vdwl -= eb_i; }
}
}
}
@ -207,60 +206,63 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
if ( dijsq < cutOffEnd2[itype][jtype] ) {
dij = sqrt (dijsq);
if ( dij < cutOffStart[itype][jtype] ) {
espo = 1.0 - dij/r0[itype][jtype];
aexpp = exp(p[itype][jtype]*espo)*A[itype][jtype];
Fr = (2.0*aexpp)*(p[itype][jtype]/r0[itype][jtype]);
qsiexpq = (QSI[itype][jtype]*QSI[itype][jtype]) * exp(2.0*q[itype][jtype]*espo);
Fb = -qsiexpq * q[itype][jtype]/r0[itype][jtype];
} else {
polyval = dij-cutOffEnd[itype][jtype];
polyval2 = polyval*polyval;
polyval3 = polyval2*polyval;
polyval4 = polyval3*polyval;
polyval5 = polyval4*polyval;
aexpp = a5[itype][jtype]*polyval5+a4[itype][jtype]*polyval4+a3[itype][jtype]*polyval3;
Fr = -2.0*(5.0*a5[itype][jtype]*polyval4+4.0*a4[itype][jtype]*polyval3+3.0*a3[itype][jtype]*polyval2);
qsiexpq = x5[itype][jtype]*polyval5+x4[itype][jtype]*polyval4+x3[itype][jtype]*polyval3;
Fb = ((5.0*x5[itype][jtype]*polyval4+4.0*x4[itype][jtype]*polyval3+3.0*x3[itype][jtype]*polyval2))*qsiexpq;
}
//if needed the repulsive energy is accumulated:
if (eflag_either) {
if (eflag_atom) {
eatom[i] += aexpp;
if (newton_pair || j < nlocal) {
eatom[j] += aexpp;
}
}
if (eflag_global) {
if (newton_pair || j < nlocal) {
eng_vdwl += 2.0*(aexpp);
} else {
eng_vdwl += aexpp;
}
}
}
//calculates the module of the pair energy between i and j
fpair = (Fb*(on_eb[i]+on_eb[j]) + Fr)/dij;
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if (dijsq < cutOffEnd2[itype][jtype]) {
dij = sqrt(dijsq);
if (dij < cutOffStart[itype][jtype]) {
espo = 1.0 - dij / r0[itype][jtype];
aexpp = exp(p[itype][jtype] * espo) * A[itype][jtype];
Fr = (2.0 * aexpp) * (p[itype][jtype] / r0[itype][jtype]);
qsiexpq = (QSI[itype][jtype] * QSI[itype][jtype]) * exp(2.0 * q[itype][jtype] * espo);
Fb = -qsiexpq * q[itype][jtype] / r0[itype][jtype];
} else {
polyval = dij - cutOffEnd[itype][jtype];
polyval2 = polyval * polyval;
polyval3 = polyval2 * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
aexpp = a5[itype][jtype] * polyval5 + a4[itype][jtype] * polyval4 +
a3[itype][jtype] * polyval3;
Fr = -2.0 *
(5.0 * a5[itype][jtype] * polyval4 + 4.0 * a4[itype][jtype] * polyval3 +
3.0 * a3[itype][jtype] * polyval2);
qsiexpq = x5[itype][jtype] * polyval5 + x4[itype][jtype] * polyval4 +
x3[itype][jtype] * polyval3;
Fb = ((5.0 * x5[itype][jtype] * polyval4 + 4.0 * x4[itype][jtype] * polyval3 +
3.0 * x3[itype][jtype] * polyval2)) *
qsiexpq;
}
//if needed the repulsive energy is accumulated:
if (eflag_either) {
if (eflag_atom) {
eatom[i] += aexpp;
if (newton_pair || j < nlocal) { eatom[j] += aexpp; }
}
if (eflag_global) {
if (newton_pair || j < nlocal) {
eng_vdwl += 2.0 * (aexpp);
} else {
eng_vdwl += aexpp;
}
}
}
//calculates the module of the pair energy between i and j
fpair = (Fb * (on_eb[i] + on_eb[j]) + Fr) / dij;
f[i][0] += del[0]*fpair;
f[i][1] += del[1]*fpair;
f[i][2] += del[2]*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= del[0]*fpair;
f[j][1] -= del[1]*fpair;
f[j][2] -= del[2]*fpair;
}
if (vflag_atom) {
ev_tally(i, j, nlocal, newton_pair,
0.0, 0.0,//Energy is tally'd in the other parts of the potential
fpair, del[0], del[1], del[2]);
}
f[i][0] += del[0] * fpair;
f[i][1] += del[1] * fpair;
f[i][2] += del[2] * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= del[0] * fpair;
f[j][1] -= del[1] * fpair;
f[j][2] -= del[2] * fpair;
}
if (vflag_atom) {
ev_tally(i, j, nlocal, newton_pair, 0.0,
0.0, //Energy is tally'd in the other parts of the potential
fpair, del[0], del[1], del[2]);
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
@ -270,43 +272,42 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
global settings
------------------------------------------------------------------------- */
void PairSMATB::settings(int narg, char **) {//reads the input script line with arguments you define
if (narg > 0) error->all(FLERR,"Illegal pair_style command: smatb accepts no options");
void PairSMATB::settings(int narg, char **)
{ //reads the input script line with arguments you define
if (narg > 0) error->all(FLERR, "Illegal pair_style command: smatb accepts no options");
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSMATB::allocate() {
void PairSMATB::allocate()
{
int n = atom->ntypes;
int natoms=atom->natoms;
int natoms = atom->natoms;
memory->create(setflag, n+1, n+1, "pair_smatb:setflag");
memory->create(setflag, n + 1, n + 1, "pair_smatb:setflag");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
}
for (int j = i; j <= n; j++) { setflag[i][j] = 0; }
}
memory->create(cutsq, n+1, n+1, "pair_smatb:cutsq");
memory->create(cutsq, n + 1, n + 1, "pair_smatb:cutsq");
//memory->create is needed to make a false nxn array on a n^2x1 line of data
memory->create(r0, n+1, n+1, "pair_smatb:r0");
memory->create(p, n+1, n+1, "pair_smatb:p");
memory->create(A, n+1, n+1, "pair_smatb:A");
memory->create(q, n+1, n+1, "pair_smatb:q");
memory->create(QSI, n+1, n+1, "pair_smatb:QSI");
memory->create(cutOffStart, n+1, n+1, "pair_smatb:cutOffStart");
memory->create(cutOffEnd, n+1, n+1, "pair_smatb:cutOffEnd");
memory->create(cutOffEnd2, n+1, n+1, "pair_smatb:cutOffEnd2");
memory->create(a3, n+1, n+1, "pair_smatb:a1");
memory->create(a4, n+1, n+1, "pair_smatb:a2");
memory->create(a5, n+1, n+1, "pair_smatb:a5");
memory->create(x3, n+1, n+1, "pair_smatb:x1");
memory->create(x4, n+1, n+1, "pair_smatb:x2");
memory->create(x5, n+1, n+1, "pair_smatb:x3");
memory->create(r0, n + 1, n + 1, "pair_smatb:r0");
memory->create(p, n + 1, n + 1, "pair_smatb:p");
memory->create(A, n + 1, n + 1, "pair_smatb:A");
memory->create(q, n + 1, n + 1, "pair_smatb:q");
memory->create(QSI, n + 1, n + 1, "pair_smatb:QSI");
memory->create(cutOffStart, n + 1, n + 1, "pair_smatb:cutOffStart");
memory->create(cutOffEnd, n + 1, n + 1, "pair_smatb:cutOffEnd");
memory->create(cutOffEnd2, n + 1, n + 1, "pair_smatb:cutOffEnd2");
memory->create(a3, n + 1, n + 1, "pair_smatb:a1");
memory->create(a4, n + 1, n + 1, "pair_smatb:a2");
memory->create(a5, n + 1, n + 1, "pair_smatb:a5");
memory->create(x3, n + 1, n + 1, "pair_smatb:x1");
memory->create(x4, n + 1, n + 1, "pair_smatb:x2");
memory->create(x5, n + 1, n + 1, "pair_smatb:x3");
allocated = 1;
}
@ -315,28 +316,29 @@ void PairSMATB::allocate() {
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSMATB::coeff(int narg, char **arg) {//set coefficients for one i,j type pair
if (!allocated) {
allocate();
}
void PairSMATB::coeff(int narg, char **arg)
{ //set coefficients for one i,j type pair
if (!allocated) { allocate(); }
if (narg != 9) {
error->all(FLERR,"Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
error->all(
FLERR,
"Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
}
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
//reading parameters from input
double myr0 = utils::numeric(FLERR,arg[2],false,lmp),
myp = utils::numeric(FLERR,arg[3],false,lmp),
myq = utils::numeric(FLERR,arg[4],false,lmp),
myA = utils::numeric(FLERR,arg[5],false,lmp),
myQSI = utils::numeric(FLERR,arg[6],false,lmp),
mycutOffStart = utils::numeric(FLERR,arg[7],false,lmp),
mycutOffEnd = utils::numeric(FLERR,arg[8],false,lmp);
double myr0 = utils::numeric(FLERR, arg[2], false, lmp),
myp = utils::numeric(FLERR, arg[3], false, lmp),
myq = utils::numeric(FLERR, arg[4], false, lmp),
myA = utils::numeric(FLERR, arg[5], false, lmp),
myQSI = utils::numeric(FLERR, arg[6], false, lmp),
mycutOffStart = utils::numeric(FLERR, arg[7], false, lmp),
mycutOffEnd = utils::numeric(FLERR, arg[8], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
r0[i][j] = myr0;
p[i][j] = myp;
A[i][j] = myA;
@ -351,8 +353,7 @@ void PairSMATB::coeff(int narg, char **arg) {//set coefficients for one i,j type
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -368,41 +369,42 @@ void PairSMATB::coeff(int narg, char **arg) {//set coefficients for one i,j type
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSMATB::init_one(int i, int j) {//perform initialization for one i,j type pair
double PairSMATB::init_one(int i, int j)
{ //perform initialization for one i,j type pair
if (setflag[i][j] == 0) {
///@todo implement smatb mixing rules
cutOffStart[i][j] = MIN(cutOffStart[i][i],cutOffStart[j][j]);
cutOffEnd[i][j] = MAX(cutOffEnd[i][i],cutOffEnd[j][j]);
cutOffStart[i][j] = MIN(cutOffStart[i][i], cutOffStart[j][j]);
cutOffEnd[i][j] = MAX(cutOffEnd[i][i], cutOffEnd[j][j]);
error->all(FLERR,"All pair coeffs are not set");
error->all(FLERR, "All pair coeffs are not set");
}
//calculating the polynomial linking to zero
double es = cutOffEnd[i][j] - cutOffStart[i][j];
double es2 = es*es;
double es3 = es2*es;
double es2 = es * es;
double es3 = es2 * es;
//variables for poly for p and A
double expp = A[i][j] * exp(p[i][j]*(1.-cutOffStart[i][j]/r0[i][j]));
double ap = -1./es3;
double bp = p[i][j]/(r0[i][j]*es2);
double cp = -(p[i][j]*p[i][j])/(es*r0[i][j]*r0[i][j]);
double expp = A[i][j] * exp(p[i][j] * (1. - cutOffStart[i][j] / r0[i][j]));
double ap = -1. / es3;
double bp = p[i][j] / (r0[i][j] * es2);
double cp = -(p[i][j] * p[i][j]) / (es * r0[i][j] * r0[i][j]);
a5[i][j]= expp * (12.*ap + 6.*bp + cp)/(2.*es2);
a4[i][j]= expp * (15.*ap + 7.*bp + cp)/es;
a3[i][j]= expp * (20.*ap + 8.*bp + cp)/2.;
a5[i][j] = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
a4[i][j] = expp * (15. * ap + 7. * bp + cp) / es;
a3[i][j] = expp * (20. * ap + 8. * bp + cp) / 2.;
//variables for poly for q and qsi
double expq = QSI[i][j]*exp(q[i][j]*(1.-cutOffStart[i][j]/r0[i][j]));
double aq = -1/es3;
double bq = q[i][j]/(es2*r0[i][j]);
double cq = -(q[i][j]*q[i][j])/(es*r0[i][j]*r0[i][j]);
double expq = QSI[i][j] * exp(q[i][j] * (1. - cutOffStart[i][j] / r0[i][j]));
double aq = -1 / es3;
double bq = q[i][j] / (es2 * r0[i][j]);
double cq = -(q[i][j] * q[i][j]) / (es * r0[i][j] * r0[i][j]);
x5[i][j] = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
x4[i][j] = expq * (15. * aq + 7. * bq + cq) / es;
x3[i][j] = expq * (20. * aq + 8. * bq + cq) / 2.;
x5[i][j] = expq * (12.*aq + 6.*bq + cq)/(2.*es2);
x4[i][j] = expq * (15.*aq + 7.*bq + cq)/es;
x3[i][j] = expq * (20.*aq + 8.*bq + cq)/2.;
cutOffEnd2[i][j] = cutOffEnd[i][j] * cutOffEnd[i][j];
if ( i!=j ) {
if (i != j) {
setflag[j][i] = 1;
cutOffEnd2[j][i] = cutOffEnd2[i][j];
@ -421,26 +423,15 @@ double PairSMATB::init_one(int i, int j) {//perform initialization for one i,j t
x4[j][i] = x4[i][j];
x5[j][i] = x5[i][j];
}
#ifdef DR_DEBUG
std::cout << i << " " << j <<std::endl;
std::cout<< "r0\tp\tq\tA\tXi\tCS\tCE\tCE^2\n";
std::cout <<r0[j][i]<<"\t"<< p[j][i]<<"\t"<< q[j][i]
<<"\t"<< A[j][i]<<"\t"<< QSI[j][i]<<"\t"
<< cutOffStart[j][i]<<"\t"<< cutOffEnd[j][i]<<"\t"<< cutOffEnd2[j][i]<<"\n"
<< "Polynomial link\n"
<< "a: " << a5[i][j] <<"\t"<< a4[i][j]<<"\t" << a3[i][j]<<"\n"
<< "q: " << x5[i][j] <<"\t"<< x4[i][j]<<"\t" << x3[i][j]<<"\n";
#endif //DR_DEBUG
//returns cutOffEnd to calculate cutforce and cutsq
return cutOffEnd[i][j];
}
/* ---------------------------------------------------------------------- */
int PairSMATB::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int PairSMATB::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
int i, j, m;
m = 0;
for (i = 0; i < n; ++i) {
@ -452,33 +443,32 @@ int PairSMATB::pack_forward_comm(int n, int *list, double *buf,
/* ---------------------------------------------------------------------- */
void PairSMATB::unpack_forward_comm(int n, int first, double *buf) {
int i,m,last;
void PairSMATB::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) {
on_eb[i] = buf[m++];
}
for (i = first; i < last; ++i) { on_eb[i] = buf[m++]; }
}
/* ---------------------------------------------------------------------- */
int PairSMATB::pack_reverse_comm(int n, int first, double *buf) {
int i,m,last;
int PairSMATB::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) {
buf[m++] = on_eb[i];
}
for (i = first; i < last; ++i) { buf[m++] = on_eb[i]; }
return m;
}
/* ---------------------------------------------------------------------- */
void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf) {
int i,j,m;
void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
@ -487,93 +477,96 @@ void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf) {
}
//write binary data of this simulation:
void PairSMATB::write_restart_settings(FILE *fp) {
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag, sizeof(int),1,fp);
fwrite(&tail_flag, sizeof(int),1,fp);
void PairSMATB::write_restart_settings(FILE *fp)
{
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag, sizeof(int), 1, fp);
}
void PairSMATB::read_restart_settings(FILE *fp) {
void PairSMATB::read_restart_settings(FILE *fp)
{
int me = comm->me;
size_t result;
if (me == 0) {
result = fread(&offset_flag,sizeof(int),1,fp);
result = fread(&mix_flag, sizeof(int),1,fp);
result = fread(&tail_flag, sizeof(int),1,fp);
result = fread(&offset_flag, sizeof(int), 1, fp);
result = fread(&mix_flag, sizeof(int), 1, fp);
result = fread(&tail_flag, sizeof(int), 1, fp);
}
MPI_Bcast(&offset_flag, 1,MPI_INT,0,world);
MPI_Bcast(&mix_flag, 1,MPI_INT,0,world);
MPI_Bcast(&tail_flag, 1,MPI_INT,0,world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
}
void PairSMATB::write_restart(FILE *fp) {
void PairSMATB::write_restart(FILE *fp)
{
write_restart_settings(fp);
//"I J r0 p q A QSI CO_start CO_end"
int i,j;
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);
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&r0[i][j], sizeof(double),1,fp);
fwrite(&p[i][j], sizeof(double),1,fp);
fwrite(&q[i][j], sizeof(double),1,fp);
fwrite(&A[i][j], sizeof(double),1,fp);
fwrite(&QSI[i][j], sizeof(double),1,fp);
fwrite(&cutOffStart[i][j],sizeof(double),1,fp);
fwrite(&cutOffEnd[i][j], sizeof(double),1,fp);
fwrite(&r0[i][j], sizeof(double), 1, fp);
fwrite(&p[i][j], sizeof(double), 1, fp);
fwrite(&q[i][j], sizeof(double), 1, fp);
fwrite(&A[i][j], sizeof(double), 1, fp);
fwrite(&QSI[i][j], sizeof(double), 1, fp);
fwrite(&cutOffStart[i][j], sizeof(double), 1, fp);
fwrite(&cutOffEnd[i][j], sizeof(double), 1, fp);
}
}
//maybe we need to save also the values of the various polynomials
}
void PairSMATB::read_restart(FILE *fp) {
void PairSMATB::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
size_t result;
int i,j;
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) {
result = fread(&setflag[i][j],sizeof(int),1,fp);
}
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (me == 0) { result = fread(&setflag[i][j], sizeof(int), 1, fp); }
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
result = fread(&r0[i][j], sizeof(double),1,fp);
result = fread(&p[i][j], sizeof(double),1,fp);
result = fread(&q[i][j], sizeof(double),1,fp);
result = fread(&A[i][j], sizeof(double),1,fp);
result = fread(&QSI[i][j], sizeof(double),1,fp);
result = fread(&cutOffStart[i][j],sizeof(double),1,fp);
result = fread(&cutOffEnd[i][j], sizeof(double),1,fp);
result = fread(&r0[i][j], sizeof(double), 1, fp);
result = fread(&p[i][j], sizeof(double), 1, fp);
result = fread(&q[i][j], sizeof(double), 1, fp);
result = fread(&A[i][j], sizeof(double), 1, fp);
result = fread(&QSI[i][j], sizeof(double), 1, fp);
result = fread(&cutOffStart[i][j], sizeof(double), 1, fp);
result = fread(&cutOffEnd[i][j], sizeof(double), 1, fp);
}
MPI_Bcast(&r0[i][j], 1,MPI_DOUBLE,0,world);
MPI_Bcast(&p[i][j], 1,MPI_DOUBLE,0,world);
MPI_Bcast(&q[i][j], 1,MPI_DOUBLE,0,world);
MPI_Bcast(&A[i][j], 1,MPI_DOUBLE,0,world);
MPI_Bcast(&QSI[i][j], 1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutOffStart[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutOffEnd[i][j], 1,MPI_DOUBLE,0,world);
MPI_Bcast(&r0[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&p[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&q[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&A[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&QSI[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffStart[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffEnd[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
void PairSMATB::write_data(FILE *fp) {
void PairSMATB::write_data(FILE *fp)
{
//smatb needs I J r0 p q A QSI CO_start CO_end
for (int i = 1; i <= atom->ntypes; i++) {
fprintf(fp,"%d %g %g %g %g %g %g %g\n",
i, r0[i][i], p[i][i], q[i][i], A[i][i], QSI[i][i], cutOffStart[i][i], cutOffEnd[i][i]);
fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, r0[i][i], p[i][i], q[i][i], A[i][i], QSI[i][i],
cutOffStart[i][i], cutOffEnd[i][i]);
}
}
void PairSMATB::write_data_all(FILE *fp) {
void PairSMATB::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
fprintf(fp,"%d %d %g %g %g %g %g %g %g\n",
i, j, r0[i][j], p[i][j], q[i][j], A[i][j], QSI[i][j], cutOffStart[i][j], cutOffEnd[i][j]);
fprintf(fp, "%d %d %g %g %g %g %g %g %g\n", i, j, r0[i][j], p[i][j], q[i][j], A[i][j],
QSI[i][j], cutOffStart[i][j], cutOffEnd[i][j]);
}
}
}

80
src/pair_smatb.h
View File

@ -26,44 +26,45 @@ PairStyle(smatb,PairSMATB)
namespace LAMMPS_NS {
class PairSMATB : public Pair {
public:
PairSMATB(class LAMMPS *);
virtual ~PairSMATB();
void compute(int, int);//workhorse routine that computes pairwise interactions
/*
class PairSMATB : public Pair {
public:
PairSMATB(class LAMMPS *);
virtual ~PairSMATB();
void compute(int, int); //workhorse routine that computes pairwise interactions
/*
void compute_inner();
void compute_middle();
void compute_outer(int, int);*/
void settings(int, char **);//reads the input script line with arguments you define
void coeff(int, char **);//set coefficients for one i,j type pair
//void init_style();//initialization specific to this pair style
double init_one(int, int);//perform initialization for one i,j type pair
//double single(int, int, int, int, double, double, double, double &);//force and energy of a single pairwise interaction between 2 atoms
void compute_outer(int, int);
*/
void settings(int, char **); //reads the input script line with arguments you define
void coeff(int, char **); //set coefficients for one i,j type pair
//void init_style();//initialization specific to this pair style
double init_one(int, int); //perform initialization for one i,j type pair
//double single(int, int, int, int, double, double, double, double &);//force and energy of a single pairwise interaction between 2 atoms
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
virtual int pack_forward_comm(int, int *, double *, int, int *);
virtual void unpack_forward_comm(int, int, double *);
virtual int pack_reverse_comm(int, int, double *);
virtual void unpack_reverse_comm(int, int *, double *);
protected:
virtual void allocate();
int nmax; // allocated size of per-atom arrays
double *on_eb; //allocated to store up caclulation values
double **r0; // interaction radius, user-given
double **p, **A, **q, **QSI; // parameters user-given
double **cutOffStart, **cutOffEnd;//cut offs, user given
double **cutOffEnd2; //squared cut off end, calculated
double **a3, **a4, **a5; //polynomial for cutoff linking to zero: Ae^p substitution
double **x3, **x4, **x5; //polynomial for cutoff linking to zero: QSIe^q substitution
/* latex form of the potential (R_c is cutOffEnd, \Xi is QSI):
virtual int pack_forward_comm(int, int *, double *, int, int *);
virtual void unpack_forward_comm(int, int, double *);
virtual int pack_reverse_comm(int, int, double *);
virtual void unpack_reverse_comm(int, int *, double *);
protected:
virtual void allocate();
int nmax; // allocated size of per-atom arrays
double *on_eb; //allocated to store up caclulation values
double **r0; // interaction radius, user-given
double **p, **A, **q, **QSI; // parameters user-given
double **cutOffStart, **cutOffEnd; //cut offs, user given
double **cutOffEnd2; //squared cut off end, calculated
double **a3, **a4, **a5; //polynomial for cutoff linking to zero: Ae^p substitution
double **x3, **x4, **x5; //polynomial for cutoff linking to zero: QSIe^q substitution
/* latex form of the potential (R_c is cutOffEnd, \Xi is QSI):
E_i =
\sum_{j,R_{ij}\leq R_c} A e^{-p \lrt{\frac{R_{ij}}{R_{0}}-1}}
@ -71,14 +72,14 @@ namespace LAMMPS_NS {
NB::this form does not have the polynomial link to 0 for the cut off
*/
};
};
}
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
/* ERROR/WARNING messages:
E: Illegal ... command
@ -90,9 +91,4 @@ namespace LAMMPS_NS {
Self-explanatory. Check the input script or data file.
E: Cannot open EAM potential file %s
The specified EAM potential file cannot be opened. Check that the
path and name are correct.
*/

493
src/pair_smatb_single.cpp
View File

@ -1,58 +1,57 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This pair style is written by Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#include "pair_smatb_single.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_smatb_single.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
//#define DR_DEBUG
#ifdef DR_DEBUG
#include <iostream>
#include <iomanip>
#endif //DR_DEBUG
using namespace LAMMPS_NS;
#define MAXLINE 1024
PairSMATBSingle::PairSMATBSingle(LAMMPS *lmp)
: Pair(lmp),
nmax(0),
on_eb(nullptr),
r0(0),
p(0),
A(0),
q(0),
QSI(0),
cutOffStart(0),
cutOffEnd(0),
cutOffEnd2(0),
a3(0), a4(0), a5(0),
x3(0), x4(0), x5(0) {
single_enable = 0; // 1 if single() routine exists
restartinfo = 1; // 1 if pair style writes restart info
respa_enable = 0; // 1 if inner/middle/outer rRESPA routines
one_coeff = 0; // 1 if allows only one coeff * * call
manybody_flag = 1; // 1 if a manybody potential
no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
writedata = 1; // 1 if writes coeffs to data file
ghostneigh = 0; // 1 if pair style needs neighbors of ghosts
PairSMATBSingle::PairSMATBSingle(LAMMPS *lmp) :
Pair(lmp), nmax(0), on_eb(nullptr), r0(0), p(0), A(0), q(0), QSI(0), cutOffStart(0),
cutOffEnd(0), cutOffEnd2(0), a3(0), a4(0), a5(0), x3(0), x4(0), x5(0)
{
single_enable = 0; // 1 if single() routine exists
restartinfo = 1; // 1 if pair style writes restart info
respa_enable = 0; // 1 if inner/middle/outer rRESPA routines
one_coeff = 0; // 1 if allows only one coeff * * call
manybody_flag = 1; // 1 if a manybody potential
no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
writedata = 1; // 1 if writes coeffs to data file
ghostneigh = 0; // 1 if pair style needs neighbors of ghosts
// set comm size needed by this Pair
comm_forward = 1;
comm_reverse = 1;
}
PairSMATBSingle::~PairSMATBSingle() {
if (copymode) {
return;
}
PairSMATBSingle::~PairSMATBSingle()
{
if (copymode) { return; }
memory->destroy(on_eb);
if (allocated) {
memory->destroy(setflag);
@ -60,17 +59,18 @@ PairSMATBSingle::~PairSMATBSingle() {
}
}
void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that computes pairwise interactions
void PairSMATBSingle::compute(int eflag, int vflag)
{ //workhorse routine that computes pairwise interactions
//eflag means compute energy
//vflag means compute virial
int i,j,ii,jj,jnum;//,itype,jtype;
double xtmp,ytmp,ztmp,del[3],fpair;
double dijsq,dij;
int i, j, ii, jj, jnum; //,itype,jtype;
double xtmp, ytmp, ztmp, del[3], fpair;
double dijsq, dij;
double espo, aexpp, qsiexpq, eb_i, Fb, Fr;
double polyval, polyval2, polyval3, polyval4, polyval5;
//sets up the flags for energy caclulations
if (eflag || vflag) {
ev_setup(eflag,vflag);
ev_setup(eflag, vflag);
eng_vdwl = 0;
} else {
evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
@ -81,14 +81,14 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
if (atom->nmax > nmax) {
nmax = atom->nmax;
memory->grow(on_eb,nmax,"pair_smatb:on_eb");
memory->grow(on_eb, nmax, "pair_smatb:on_eb");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair;
@ -107,39 +107,39 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
//FIRST LOOP: CALCULATES the squared bounding energy and accumulate it in on_eb for each atom
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
//itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; ++jj) {
j = jlist[jj];
j = jlist[jj];
j &= NEIGHMASK;
//jtype = type[j];
del[0] = xtmp - x[j][0];
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if ( dijsq < cutOffEnd2 ) {
dij = sqrt (dijsq);
if ( dij < cutOffStart ) {
qsiexpq = (QSI*QSI) * exp(2.0*q*(1.0 - dij/r0));
} else {
polyval = dij-cutOffEnd;
polyval3 = polyval*polyval*polyval;
polyval4 = polyval3*polyval;
polyval5 = polyval4*polyval;
qsiexpq = x5*polyval5+x4*polyval4+x3*polyval3;
qsiexpq = qsiexpq* qsiexpq;
}
on_eb[i]+=qsiexpq;
//if (newton_pair || j < nlocal) {
on_eb[j]+=qsiexpq;
//}
if (dijsq < cutOffEnd2) {
dij = sqrt(dijsq);
if (dij < cutOffStart) {
qsiexpq = (QSI * QSI) * exp(2.0 * q * (1.0 - dij / r0));
} else {
polyval = dij - cutOffEnd;
polyval3 = polyval * polyval * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
qsiexpq = x5 * polyval5 + x4 * polyval4 + x3 * polyval3;
qsiexpq = qsiexpq * qsiexpq;
}
on_eb[i] += qsiexpq;
//if (newton_pair || j < nlocal) {
on_eb[j] += qsiexpq;
//}
}
}
}
@ -152,22 +152,18 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
// the store the reciprocal in on_eb in order to not do it in the SECOND LOOP
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (i < nlocal){
eb_i=sqrt(on_eb[i]);
if (eb_i!=0.0) {
on_eb[i]=1.0/eb_i;
i = ilist[ii];
if (i < nlocal) {
eb_i = sqrt(on_eb[i]);
if (eb_i != 0.0) {
on_eb[i] = 1.0 / eb_i;
} else {
on_eb[i] = 0.0;
}
//if needed the bounding energy is accumulated:
if (eflag_either) {
if (eflag_atom) {
eatom[i] -= eb_i;
}
if (eflag_global) {
eng_vdwl -= eb_i;
}
if (eflag_atom) { eatom[i] -= eb_i; }
if (eflag_global) { eng_vdwl -= eb_i; }
}
}
}
@ -192,60 +188,57 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
if ( dijsq < cutOffEnd2 ) {
dij = sqrt (dijsq);
if ( dij < cutOffStart ) {
espo = 1.0 - dij/r0;
aexpp = exp(p*espo)*A;
Fr = (2.0*aexpp)*(p/r0);
qsiexpq = (QSI*QSI) * exp(2.0*q*espo);
Fb = -qsiexpq * q/r0;
} else {
polyval = dij-cutOffEnd;
polyval2 = polyval*polyval;
polyval3 = polyval2*polyval;
polyval4 = polyval3*polyval;
polyval5 = polyval4*polyval;
aexpp = a5*polyval5+a4*polyval4+a3*polyval3;
Fr = -2.0*(5.0*a5*polyval4+4.0*a4*polyval3+3.0*a3*polyval2);
qsiexpq = x5*polyval5+x4*polyval4+x3*polyval3;
Fb = ((5.0*x5*polyval4+4.0*x4*polyval3+3.0*x3*polyval2))*qsiexpq;
}
//if needed the repulsive energy is accumulated:
if (eflag_either) {
if (eflag_atom) {
eatom[i] += aexpp;
if (newton_pair || j < nlocal) {
eatom[j] += aexpp;
}
}
if (eflag_global) {
if (newton_pair || j < nlocal) {
eng_vdwl += 2.0*(aexpp);
} else {
eng_vdwl += aexpp;
}
}
}
//calculates the module of the pair energy between i and j
fpair = (Fb*(on_eb[i]+on_eb[j]) + Fr)/dij;
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if (dijsq < cutOffEnd2) {
dij = sqrt(dijsq);
if (dij < cutOffStart) {
espo = 1.0 - dij / r0;
aexpp = exp(p * espo) * A;
Fr = (2.0 * aexpp) * (p / r0);
qsiexpq = (QSI * QSI) * exp(2.0 * q * espo);
Fb = -qsiexpq * q / r0;
} else {
polyval = dij - cutOffEnd;
polyval2 = polyval * polyval;
polyval3 = polyval2 * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
aexpp = a5 * polyval5 + a4 * polyval4 + a3 * polyval3;
Fr = -2.0 * (5.0 * a5 * polyval4 + 4.0 * a4 * polyval3 + 3.0 * a3 * polyval2);
qsiexpq = x5 * polyval5 + x4 * polyval4 + x3 * polyval3;
Fb = ((5.0 * x5 * polyval4 + 4.0 * x4 * polyval3 + 3.0 * x3 * polyval2)) * qsiexpq;
}
//if needed the repulsive energy is accumulated:
if (eflag_either) {
if (eflag_atom) {
eatom[i] += aexpp;
if (newton_pair || j < nlocal) { eatom[j] += aexpp; }
}
if (eflag_global) {
if (newton_pair || j < nlocal) {
eng_vdwl += 2.0 * (aexpp);
} else {
eng_vdwl += aexpp;
}
}
}
//calculates the module of the pair energy between i and j
fpair = (Fb * (on_eb[i] + on_eb[j]) + Fr) / dij;
f[i][0] += del[0]*fpair;
f[i][1] += del[1]*fpair;
f[i][2] += del[2]*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= del[0]*fpair;
f[j][1] -= del[1]*fpair;
f[j][2] -= del[2]*fpair;
}
if (vflag_atom) {
ev_tally(i, j, nlocal, newton_pair,
0.0, 0.0,//Energy is tally'd in the other parts of the potential
fpair, del[0], del[1], del[2]);
}
f[i][0] += del[0] * fpair;
f[i][1] += del[1] * fpair;
f[i][2] += del[2] * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= del[0] * fpair;
f[j][1] -= del[1] * fpair;
f[j][2] -= del[2] * fpair;
}
if (vflag_atom) {
ev_tally(i, j, nlocal, newton_pair, 0.0,
0.0, //Energy is tally'd in the other parts of the potential
fpair, del[0], del[1], del[2]);
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
@ -255,26 +248,26 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
global settings
------------------------------------------------------------------------- */
void PairSMATBSingle::settings(int narg, char **) {//reads the input script line with arguments you define
if (narg > 0) error->all(FLERR,"Illegal pair_style command: smatb accepts no options");
void PairSMATBSingle::settings(int narg, char **)
{ //reads the input script line with arguments you define
if (narg > 0) error->all(FLERR, "Illegal pair_style command: smatb accepts no options");
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSMATBSingle::allocate() {
void PairSMATBSingle::allocate()
{
int n = atom->ntypes;
int natoms=atom->natoms;
int natoms = atom->natoms;
memory->create(setflag, n+1, n+1, "pair_smatb:setflag");
memory->create(setflag, n + 1, n + 1, "pair_smatb:setflag");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
}
for (int j = i; j <= n; j++) { setflag[i][j] = 0; }
}
memory->create(cutsq, n+1, n+1, "pair_smatb:cutsq");
memory->create(cutsq, n + 1, n + 1, "pair_smatb:cutsq");
allocated = 1;
}
@ -283,28 +276,29 @@ void PairSMATBSingle::allocate() {
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSMATBSingle::coeff(int narg, char **arg) {//set coefficients for one i,j type pair
if (!allocated) {
allocate();
}
void PairSMATBSingle::coeff(int narg, char **arg)
{ //set coefficients for one i,j type pair
if (!allocated) { allocate(); }
if (narg != 9) {
error->all(FLERR,"Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
error->all(
FLERR,
"Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
}
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
//reading parameters from input
double myr0 = utils::numeric(FLERR,arg[2],false,lmp),
myp = utils::numeric(FLERR,arg[3],false,lmp),
myq = utils::numeric(FLERR,arg[4],false,lmp),
myA = utils::numeric(FLERR,arg[5],false,lmp),
myQSI = utils::numeric(FLERR,arg[6],false,lmp),
mycutOffStart = utils::numeric(FLERR,arg[7],false,lmp),
mycutOffEnd = utils::numeric(FLERR,arg[8],false,lmp);
double myr0 = utils::numeric(FLERR, arg[2], false, lmp),
myp = utils::numeric(FLERR, arg[3], false, lmp),
myq = utils::numeric(FLERR, arg[4], false, lmp),
myA = utils::numeric(FLERR, arg[5], false, lmp),
myQSI = utils::numeric(FLERR, arg[6], false, lmp),
mycutOffStart = utils::numeric(FLERR, arg[7], false, lmp),
mycutOffEnd = utils::numeric(FLERR, arg[8], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
r0 = myr0;
p = myp;
A = myA;
@ -318,8 +312,7 @@ void PairSMATBSingle::coeff(int narg, char **arg) {//set coefficients for one i,
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -335,35 +328,36 @@ void PairSMATBSingle::coeff(int narg, char **arg) {//set coefficients for one i,
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSMATBSingle::init_one(int i, int j) {//perform initialization for one i,j type pair
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
double PairSMATBSingle::init_one(int i, int j)
{ //perform initialization for one i,j type pair
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
//calculating the polynomial linking to zero
double es = cutOffEnd - cutOffStart;
double es2 = es*es;
double es3 = es2*es;
double es2 = es * es;
double es3 = es2 * es;
//variables for poly for p and A
double expp = A * exp(p*(1.-cutOffStart/r0));
double ap = -1./es3;
double bp = p/(r0*es2);
double cp = -(p*p)/(es*r0*r0);
double expp = A * exp(p * (1. - cutOffStart / r0));
double ap = -1. / es3;
double bp = p / (r0 * es2);
double cp = -(p * p) / (es * r0 * r0);
a5= expp * (12.*ap + 6.*bp + cp)/(2.*es2);
a4= expp * (15.*ap + 7.*bp + cp)/es;
a3= expp * (20.*ap + 8.*bp + cp)/2.;
a5 = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
a4 = expp * (15. * ap + 7. * bp + cp) / es;
a3 = expp * (20. * ap + 8. * bp + cp) / 2.;
//variables for poly for q and qsi
double expq = QSI*exp(q*(1.-cutOffStart/r0));
double aq = -1/es3;
double bq = q/(es2*r0);
double cq = -(q*q)/(es*r0*r0);
double expq = QSI * exp(q * (1. - cutOffStart / r0));
double aq = -1 / es3;
double bq = q / (es2 * r0);
double cq = -(q * q) / (es * r0 * r0);
x5 = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
x4 = expq * (15. * aq + 7. * bq + cq) / es;
x3 = expq * (20. * aq + 8. * bq + cq) / 2.;
x5 = expq * (12.*aq + 6.*bq + cq)/(2.*es2);
x4 = expq * (15.*aq + 7.*bq + cq)/es;
x3 = expq * (20.*aq + 8.*bq + cq)/2.;
cutOffEnd2 = cutOffEnd * cutOffEnd;
if ( i!=j ) {
if (i != j) {
setflag[j][i] = 1;
cutOffEnd2 = cutOffEnd2;
@ -382,26 +376,15 @@ double PairSMATBSingle::init_one(int i, int j) {//perform initialization for one
x4 = x4;
x5 = x5;
}
#ifdef DR_DEBUG
std::cout << i << " " << j <<std::endl;
std::cout<< "r0\tp\tq\tA\tXi\tCS\tCE\n";
std::cout <<r0<<"\t"<< p<<"\t"<< q
<<"\t"<< A<<"\t"<< QSI<<"\t"
<< cutOffStart<<"\t"<< cutOffEnd<<"\n"
<< "Polynomial link\n"
<< "a: " << a5 <<"\t"<< a4<<"\t" << a3<<"\n"
<< "q: " << x5 <<"\t"<< x4<<"\t" << x3<<"\n";
#endif //DR_DEBUG
//returns cutOffEnd to calculate cutforce and cutsq
return cutOffEnd;
}
/* ---------------------------------------------------------------------- */
int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
int i, j, m;
m = 0;
for (i = 0; i < n; ++i) {
@ -413,33 +396,32 @@ int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf,
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::unpack_forward_comm(int n, int first, double *buf) {
int i,m,last;
void PairSMATBSingle::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) {
on_eb[i] = buf[m++];
}
for (i = first; i < last; ++i) { on_eb[i] = buf[m++]; }
}
/* ---------------------------------------------------------------------- */
int PairSMATBSingle::pack_reverse_comm(int n, int first, double *buf) {
int i,m,last;
int PairSMATBSingle::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) {
buf[m++] = on_eb[i];
}
for (i = first; i < last; ++i) { buf[m++] = on_eb[i]; }
return m;
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf) {
int i,j,m;
void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
@ -448,93 +430,94 @@ void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf) {
}
//write binary data of this simulation:
void PairSMATBSingle::write_restart_settings(FILE *fp) {
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag, sizeof(int),1,fp);
fwrite(&tail_flag, sizeof(int),1,fp);
void PairSMATBSingle::write_restart_settings(FILE *fp)
{
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag, sizeof(int), 1, fp);
}
void PairSMATBSingle::read_restart_settings(FILE *fp) {
void PairSMATBSingle::read_restart_settings(FILE *fp)
{
int me = comm->me;
size_t result;
if (me == 0) {
result = fread(&offset_flag,sizeof(int),1,fp);
result = fread(&mix_flag, sizeof(int),1,fp);
result = fread(&tail_flag, sizeof(int),1,fp);
result = fread(&offset_flag, sizeof(int), 1, fp);
result = fread(&mix_flag, sizeof(int), 1, fp);
result = fread(&tail_flag, sizeof(int), 1, fp);
}
MPI_Bcast(&offset_flag, 1,MPI_INT,0,world);
MPI_Bcast(&mix_flag, 1,MPI_INT,0,world);
MPI_Bcast(&tail_flag, 1,MPI_INT,0,world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
}
void PairSMATBSingle::write_restart(FILE *fp) {
void PairSMATBSingle::write_restart(FILE *fp)
{
write_restart_settings(fp);
//"I J r0 p q A QSI CO_start CO_end"
int i,j;
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);
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&r0, sizeof(double),1,fp);
fwrite(&p, sizeof(double),1,fp);
fwrite(&q, sizeof(double),1,fp);
fwrite(&A, sizeof(double),1,fp);
fwrite(&QSI, sizeof(double),1,fp);
fwrite(&cutOffStart,sizeof(double),1,fp);
fwrite(&cutOffEnd, sizeof(double),1,fp);
fwrite(&r0, sizeof(double), 1, fp);
fwrite(&p, sizeof(double), 1, fp);
fwrite(&q, sizeof(double), 1, fp);
fwrite(&A, sizeof(double), 1, fp);
fwrite(&QSI, sizeof(double), 1, fp);
fwrite(&cutOffStart, sizeof(double), 1, fp);
fwrite(&cutOffEnd, sizeof(double), 1, fp);
}
}
//maybe we need to save also the values of the various polynomials
}
void PairSMATBSingle::read_restart(FILE *fp) {
void PairSMATBSingle::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
size_t result;
int i,j;
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) {
result = fread(&setflag[i][j],sizeof(int),1,fp);
}
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (me == 0) { result = fread(&setflag[i][j], sizeof(int), 1, fp); }
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
result = fread(&r0, sizeof(double),1,fp);
result = fread(&p, sizeof(double),1,fp);
result = fread(&q, sizeof(double),1,fp);
result = fread(&A, sizeof(double),1,fp);
result = fread(&QSI, sizeof(double),1,fp);
result = fread(&cutOffStart,sizeof(double),1,fp);
result = fread(&cutOffEnd, sizeof(double),1,fp);
result = fread(&r0, sizeof(double), 1, fp);
result = fread(&p, sizeof(double), 1, fp);
result = fread(&q, sizeof(double), 1, fp);
result = fread(&A, sizeof(double), 1, fp);
result = fread(&QSI, sizeof(double), 1, fp);
result = fread(&cutOffStart, sizeof(double), 1, fp);
result = fread(&cutOffEnd, sizeof(double), 1, fp);
}
MPI_Bcast(&r0, 1,MPI_DOUBLE,0,world);
MPI_Bcast(&p, 1,MPI_DOUBLE,0,world);
MPI_Bcast(&q, 1,MPI_DOUBLE,0,world);
MPI_Bcast(&A, 1,MPI_DOUBLE,0,world);
MPI_Bcast(&QSI, 1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutOffStart,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cutOffEnd, 1,MPI_DOUBLE,0,world);
MPI_Bcast(&r0, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&p, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&q, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&A, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&QSI, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffStart, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffEnd, 1, MPI_DOUBLE, 0, world);
}
}
}
void PairSMATBSingle::write_data(FILE *fp) {
void PairSMATBSingle::write_data(FILE *fp)
{
//smatb needs I J r0 p q A QSI CO_start CO_end
for (int i = 1; i <= atom->ntypes; i++) {
fprintf(fp,"%d %g %g %g %g %g %g %g\n",
i, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
}
}
void PairSMATBSingle::write_data_all(FILE *fp) {
void PairSMATBSingle::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
fprintf(fp,"%d %d %g %g %g %g %g %g %g\n",
i, j, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
fprintf(fp, "%d %d %g %g %g %g %g %g %g\n", i, j, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
}
}
}

82
src/pair_smatb_single.h
View File

@ -10,7 +10,7 @@
See the README file in the top-level LAMMPS directory.
This style is written by Daniele Rapetti (iximiel@gmail.com)
This pair style is written by Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
@ -26,53 +26,51 @@ PairStyle(smatb/single,PairSMATBSingle)
namespace LAMMPS_NS {
class PairSMATBSingle : public Pair {
public:
// public variables so USER-ATC package can access them
PairSMATBSingle(class LAMMPS *);
virtual ~PairSMATBSingle();
void compute(int, int);//workhorse routine that computes pairwise interactions
/*
class PairSMATBSingle : public Pair {
public:
PairSMATBSingle(class LAMMPS *);
virtual ~PairSMATBSingle();
void compute(int, int); //workhorse routine that computes pairwise interactions
/*
void compute_inner();
void compute_middle();
void compute_outer(int, int);*/
void settings(int, char **);//reads the input script line with arguments you define
void coeff(int, char **);//set coefficients for one i,j type pair
//void init_style();//initialization specific to this pair style
double init_one(int, int);//perform initialization for one i,j type pair
//double single(int, int, int, int, double, double, double, double &);//force and energy of a single pairwise interaction between 2 atoms
void compute_outer(int, int);
*/
void settings(int, char **); //reads the input script line with arguments you define
void coeff(int, char **); //set coefficients for one i,j type pair
//void init_style();//initialization specific to this pair style
double init_one(int, int); //perform initialization for one i,j type pair
virtual void write_restart(FILE *);
virtual void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual void write_data(FILE *);
virtual void write_data_all(FILE *);
virtual void write_restart(FILE *);
virtual void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual void write_data(FILE *);
virtual void write_data_all(FILE *);
virtual int pack_forward_comm(int, int *, double *, int, int *);
virtual void unpack_forward_comm(int, int, double *);
virtual int pack_reverse_comm(int , int , double *);
virtual void unpack_reverse_comm(int , int *, double *);
protected:
virtual void allocate();
int nmax; // allocated size of per-atom arrays
double *on_eb; //allocated to store up caclulation values
double r0; // interaction radius, user-given
double p, A, q, QSI; // parameters user-given
double cutOffStart, cutOffEnd;//cut offs, user given
double cutOffEnd2; //squared cut off end, calculated
double a3, a4, a5; //polynomial for cutoff linking to zero: Ae^p substitution
double x3, x4, x5; //polynomial for cutoff linking to zero: QSIe^q substitution
};
virtual int pack_forward_comm(int, int *, double *, int, int *);
virtual void unpack_forward_comm(int, int, double *);
virtual int pack_reverse_comm(int, int, double *);
virtual void unpack_reverse_comm(int, int *, double *);
}
protected:
virtual void allocate();
int nmax; // allocated size of per-atom arrays
double *on_eb; //allocated to store up caclulation values
double r0; // interaction radius, user-given
double p, A, q, QSI; // parameters user-given
double cutOffStart, cutOffEnd; //cut offs, user given
double cutOffEnd2; //squared cut off end, calculated
double a3, a4, a5; //polynomial for cutoff linking to zero: Ae^p substitution
double x3, x4, x5; //polynomial for cutoff linking to zero: QSIe^q substitution
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
/* ERROR/WARNING messages:
E: Illegal ... command
@ -83,10 +81,4 @@ namespace LAMMPS_NS {
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Cannot open EAM potential file %s
The specified EAM potential file cannot be opened. Check that the
path and name are correct.
*/