Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80791791
pair_polymorphic.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 2, 15:01
Size
28 KB
Mime Type
text/x-c
Expires
Wed, Sep 4, 15:01 (2 d)
Engine
blob
Format
Raw Data
Handle
20286258
Attached To
rLAMMPS lammps
pair_polymorphic.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Reese Jones, Xiaowang Zhou (SNL)
This modifies from pair_tersoff.cpp by Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_polymorphic.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
one_coeff = 1;
nelements = 0;
elements = NULL;
pairParameters = NULL;
tripletParameters = NULL;
elem2param = NULL;
elem3param = NULL;
map = NULL;
epsilon = 0.0;
neighsize = 0;
firstneighV = NULL;
firstneighW = NULL;
firstneighW1 = NULL;
delxV = NULL;
delyV = NULL;
delzV = NULL;
drV = NULL;
delxW = NULL;
delyW = NULL;
delzW = NULL;
drW = NULL;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairPolymorphic::~PairPolymorphic()
{
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] match;
memory->destroy(pairParameters);
memory->destroy(tripletParameters);
memory->destroy(elem2param);
memory->destroy(elem3param);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
delete [] firstneighV;
delete [] firstneighW;
delete [] firstneighW1;
delete [] delxV;
delete [] delyV;
delete [] delzV;
delete [] drV;
delete [] delxW;
delete [] delyW;
delete [] delzW;
delete [] drW;
}
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::compute(int eflag, int vflag)
{
tagint itag,jtag;
int i,j,k,ii,jj,kk,kk1,inum,jnum;
int itype,jtype,ktype;
int iparam_ii,iparam_jj,iparam_kk,iparam_ij,iparam_ik,iparam_ijk;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r0,r1,r2;
double delr1[3],delr2[3],fi[3],fj[3],fk[3];
double zeta_ij,prefactor,wfac,pfac,gfac,fa,fa_d,bij,bij_d;
double costheta;
int *ilist,*jlist,*numneigh,**firstneigh;
double emb;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = vflag_atom = 0;
double **x = atom->x;
double **f = atom->f;
tagint *tag = atom->tag;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over full neighbor list of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
if (neighsize < jnum) {
delete [] firstneighV;
delete [] delxV;
delete [] delyV;
delete [] delzV;
delete [] drV;
delete [] firstneighW;
delete [] delxW;
delete [] delyW;
delete [] delzW;
delete [] drW;
delete [] firstneighW1;
neighsize = jnum + 20;
firstneighV = new int[neighsize];
delxV = new double[neighsize];
delyV = new double[neighsize];
delzV = new double[neighsize];
drV = new double[neighsize];
firstneighW = new int[neighsize];
delxW = new double[neighsize];
delyW = new double[neighsize];
delzW = new double[neighsize];
drW = new double[neighsize];
firstneighW1 = new int[neighsize];
}
if (eta) {
iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii];
emb = (p.F)->get_vmax();
}
numneighV = -1;
numneighW = -1;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq >= cutmaxsq) continue;
r0 = sqrt(rsq);
iparam_ij = elem2param[itype][jtype];
PairParameters & p = pairParameters[iparam_ij];
// do not include the neighbor if get_vmax() <= epsilon because the function is near zero
if (eta) {
if (emb > epsilon) {
iparam_jj = elem2param[jtype][jtype];
PairParameters & q = pairParameters[iparam_jj];
if (rsq < (q.W)->get_xmaxsq() && (q.W)->get_vmax() > epsilon) {
numneighW = numneighW + 1;
firstneighW[numneighW] = j;
delxW[numneighW] = delx;
delyW[numneighW] = dely;
delzW[numneighW] = delz;
drW[numneighW] = r0;
}
}
} else {
if ((p.F)->get_vmax() > epsilon) {
if (rsq < (p.V)->get_xmaxsq() && (p.V)->get_vmax() > epsilon) {
numneighV = numneighV + 1;
firstneighV[numneighV] = j;
delxV[numneighV] = delx;
delyV[numneighV] = dely;
delzV[numneighV] = delz;
drV[numneighV] = r0;
}
if (rsq < (p.W)->get_xmaxsq() && (p.W)->get_vmax() > epsilon) {
numneighW = numneighW + 1;
firstneighW[numneighW] = j;
delxW[numneighW] = delx;
delyW[numneighW] = dely;
delzW[numneighW] = delz;
drW[numneighW] = r0;
}
}
}
// two-body interactions, skip half of them
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < x[i][2]) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
if (rsq >= (p.U)->get_xmaxsq() || (p.U)->get_vmax() <= epsilon) continue;
(p.U)->value(r0,evdwl,eflag,fpair,1);
fpair = -fpair/r0;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
if (eta) {
if (emb > epsilon) {
iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii];
// accumulate bondorder zeta for each i-j interaction via loop over k
zeta_ij = 0.0;
for (kk = 0; kk <= numneighW; kk++) {
k = firstneighW[kk];
ktype = map[type[k]];
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,1,fpair,0);
zeta_ij += wfac;
}
// pairwise force due to zeta
(p.F)->value(zeta_ij,bij,1,bij_d,1);
prefactor = 0.5* bij_d;
if (eflag) evdwl = -0.5*bij;
if (evflag) ev_tally(i,i,nlocal,newton_pair,evdwl,0.0,0.0,delx,dely,delz);
// attractive term via loop over k
for (kk = 0; kk <= numneighW; kk++) {
k = firstneighW[kk];
ktype = map[type[k]];
delr2[0] = -delxW[kk];
delr2[1] = -delyW[kk];
delr2[2] = -delzW[kk];
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,0,fpair,1);
fpair = -prefactor*fpair/drW[kk];
f[i][0] += delr2[0]*fpair;
f[i][1] += delr2[1]*fpair;
f[i][2] += delr2[2]*fpair;
f[k][0] -= delr2[0]*fpair;
f[k][1] -= delr2[1]*fpair;
f[k][2] -= delr2[2]*fpair;
if (vflag_atom) v_tally2(i, k, -fpair, delr2);
}
}
} else {
for (jj = 0; jj <= numneighV; jj++) {
j = firstneighV[jj];
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype];
PairParameters & p = pairParameters[iparam_ij];
delr1[0] = -delxV[jj];
delr1[1] = -delyV[jj];
delr1[2] = -delzV[jj];
r1 = drV[jj];
// accumulate bondorder zeta for each i-j interaction via loop over k
zeta_ij = 0.0;
numneighW1 = -1;
for (kk = 0; kk <= numneighW; kk++) {
k = firstneighW[kk];
if (j == k) continue;
ktype = map[type[k]];
iparam_ijk = elem3param[jtype][itype][ktype];
TripletParameters & trip = tripletParameters[iparam_ijk];
if ((trip.G)->get_vmax() <= epsilon) continue;
numneighW1 = numneighW1 + 1;
firstneighW1[numneighW1] = kk;
delr2[0] = -delxW[kk];
delr2[1] = -delyW[kk];
delr2[2] = -delzW[kk];
r2 = drW[kk];
costheta = (delr1[0]*delr2[0] + delr1[1]*delr2[1] +
delr1[2]*delr2[2]) / (r1*r2);
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
(q.W)->value(r2,wfac,1,fpair,0);
(q.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0);
(trip.G)->value(costheta,gfac,1,fpair,0);
zeta_ij += wfac*pfac*gfac;
}
// pairwise force due to zeta
(p.V)->value(r1,fa,1,fa_d,1);
(p.F)->value(zeta_ij,bij,1,bij_d,1);
fpair = -0.5*bij*fa_d / r1;
prefactor = 0.5* fa * bij_d;
if (eflag) evdwl = -0.5*bij*fa;
f[i][0] += delr1[0]*fpair;
f[i][1] += delr1[1]*fpair;
f[i][2] += delr1[2]*fpair;
f[j][0] -= delr1[0]*fpair;
f[j][1] -= delr1[1]*fpair;
f[j][2] -= delr1[2]*fpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,-fpair,-delr1[0],-delr1[1],-delr1[2]);
// attractive term via loop over k
for (kk1 = 0; kk1 <= numneighW1; kk1++) {
kk = firstneighW1[kk1];
k = firstneighW[kk];
ktype = map[type[k]];
iparam_ijk = elem3param[jtype][itype][ktype];
TripletParameters & trip = tripletParameters[iparam_ijk];
delr2[0] = -delxW[kk];
delr2[1] = -delyW[kk];
delr2[2] = -delzW[kk];
r2 = drW[kk];
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
attractive(&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk);
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][2] += fj[2];
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
map = new int[n+1];
neighsize = 40;
firstneighV = new int[neighsize];
delxV = new double[neighsize];
delyV = new double[neighsize];
delzV = new double[neighsize];
drV = new double[neighsize];
firstneighW = new int[neighsize];
delxW = new double[neighsize];
delyW = new double[neighsize];
delzW = new double[neighsize];
drW = new double[neighsize];
firstneighW1 = new int[neighsize];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairPolymorphic::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairPolymorphic::coeff(int narg, char **arg)
{
int i,j,n;
if (!allocated) allocate();
if (narg == 4 + atom->ntypes) {
narg--;
epsilon = atof(arg[narg]);
} else if (narg != 3 + atom->ntypes) {
error->all(FLERR,"Incorrect args for pair coefficients");
}
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
// nelements = # of unique elements
// elements = list of element names
if (elements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
}
elements = new char*[atom->ntypes];
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
nelements = 0;
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
map[i-2] = j;
if (j == nelements) {
n = strlen(arg[i]) + 1;
elements[j] = new char[n];
strcpy(elements[j],arg[i]);
nelements++;
}
}
// read potential file and initialize potential parameters
read_file(arg[2]);
setup_params();
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairPolymorphic::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style polymorphic requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style polymorphic requires newton pair on");
// need a full neighbor list
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairPolymorphic::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::read_file(char *file)
{
char line[MAXLINE],*ptr;
int n;
// open file on proc 0
FILE *fp=NULL;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open polymorphic potential file %s",file);
error->one(FLERR,str);
}
// move past comments to first data line
fgets(line,MAXLINE,fp);
while (line == strchr(line,'#')) fgets(line,MAXLINE,fp);
n = strlen(line) + 1;
}
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
ptr = strtok(line," \t\n\r\f"); // 1st line, 1st token
int ntypes = atoi(ptr);
if (ntypes != nelements)
error->all(FLERR,"Incorrect number of elements in potential file");
match = new int[nelements];
ptr = strtok(NULL," \t\n\r\f"); // 1st line, 2nd token
eta = (atoi(ptr)>0) ? true:false;
// map the elements in the potential file to LAMMPS atom types
for (int i = 0; i < nelements; i++) {
if (comm->me == 0) {
fgets(line,MAXLINE,fp);
n = strlen(line) + 1;
}
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
ptr = strtok(line," \t\n\r\f"); // 1st token
ptr = strtok(NULL," \t\n\r\f"); // 2st token
ptr = strtok(NULL," \t\n\r\f"); // 3st token
int j;
for (j = 0; j < nelements; j++) {
if (strcmp(ptr,elements[j]) == 0) break;
}
if (j == nelements)
error->all(FLERR,"Element not defined in potential file");
match[i] = j;
}
// sizes
if (comm->me == 0) {
fgets(line,MAXLINE,fp);
n = strlen(line) + 1;
}
// Note: the format of this line has changed between the
// 2015-06-06 and 2015-12-09 versions of the pair style.
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
nr = ng = nx = 0;
ptr = strtok(line," \t\n\r\f"); // 1st token
if (ptr) nr = atoi(ptr);
ptr = strtok(NULL," \t\n\r\f"); // 2nd token
if (ptr) ng = atoi(ptr);
ptr = strtok(NULL," \t\n\r\f"); // 3rd token
if (ptr) nx = atoi(ptr);
ptr = strtok(NULL," \t\n\r\f"); // 4th token
if (ptr) maxX = atof(ptr);
if (ptr == NULL)
error->all(FLERR,"Potential file incompatible with this pair style version");
if ((ng == 0) || (nr == 0) || (nx == 0))
error->all(FLERR,"Error reading potential file header");
npair = nelements*(nelements+1)/2;
ntriple = nelements*nelements*nelements;
pairParameters = (PairParameters*)
memory->srealloc(pairParameters,npair*sizeof(PairParameters),
"pair:pairParameters");
tripletParameters = (TripletParameters*)
memory->srealloc(tripletParameters,ntriple*sizeof(TripletParameters),
"pair:tripletParameters");
// cutoffs
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
fgets(line,MAXLINE,fp);
n = strlen(line) + 1;
}
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
ptr = strtok(line," \t\n\r\f"); // 1st token
p.cut = atof(ptr);
p.cutsq = p.cut*p.cut;
ptr = strtok(NULL," \t\n\r\f"); // 2nd token
p.xi = (atoi(ptr)>0) ? true:false;
}
// set cutmax to max of all params
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
if (p.cut > cutmax) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
// start reading tabular functions
double * singletable = new double[nr];
for (int i = 0; i < npair; i++) { // U
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.U = new tabularFunction(nr,0.0,p.cut);
(p.U)->set_values(nr,0.0,p.cut,singletable,epsilon);
}
for (int i = 0; i < npair; i++) { // V
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.V = new tabularFunction(nr,0.0,p.cut);
(p.V)->set_values(nr,0.0,p.cut,singletable,epsilon);
}
for (int i = 0; i < npair; i++) { // W
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.W = new tabularFunction(nr,0.0,p.cut);
(p.W)->set_values(nr,0.0,p.cut,singletable,epsilon);
}
for (int i = 0; i < npair; i++) { // P
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
}
delete[] singletable;
singletable = new double[ng];
for (int i = 0; i < ntriple; i++) { // G
TripletParameters & p = tripletParameters[i];
if (comm->me == 0) {
grab(fp,ng,singletable);
}
MPI_Bcast(singletable,ng,MPI_DOUBLE,0,world);
p.G = new tabularFunction(ng,-1.0,1.0);
(p.G)->set_values(ng,-1.0,1.0,singletable,epsilon);
}
delete[] singletable;
singletable = new double[nx];
for (int i = 0; i < npair; i++) { // F
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
grab(fp,nx,singletable);
}
MPI_Bcast(singletable,nx,MPI_DOUBLE,0,world);
p.F = new tabularFunction(nx,0.0,maxX);
(p.F)->set_values(nx,0.0,maxX,singletable,epsilon);
}
delete[] singletable;
if (comm->me == 0) {
fclose(fp);
}
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::setup_params()
{
int i,j,k,n;
memory->destroy(elem2param);
memory->create(elem2param,nelements,nelements,"pair:elem2param");
memory->destroy(elem3param);
memory->create(elem3param,nelements,nelements,nelements,"pair:elem3param");
// map atom pair to parameter index
n = 0;
for (i = 0; i < nelements; i++) {
elem2param[match[i]][match[i]] = n;
n++;
}
for (i = 0; i < nelements-1; i++) {
for (j = i+1; j < nelements; j++) {
elem2param[match[i]][match[j]] = n;
elem2param[match[j]][match[i]] = n;
n++;
}
}
// map atom triplet to parameter index
n = 0;
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
elem3param[match[i]][match[j]][match[k]] = n;
n++;
}
// for debugging, call write_tables() to check the tabular functions
// if (comm->me == 0) {
// write_tables(51);
// errorX->all(FLERR,"Test potential tables");
// }
}
/* ----------------------------------------------------------------------
attractive term
------------------------------------------------------------------------- */
void PairPolymorphic::attractive(PairParameters *p, TripletParameters *trip,
double prefactor, double rij, double rik,
double *delrij, double *delrik,
double *fi, double *fj, double *fk)
{
double rij_hat[3],rik_hat[3];
double rijinv,rikinv;
rijinv = 1.0/rij;
vec3_scale(rijinv,delrij,rij_hat);
rikinv = 1.0/rik;
vec3_scale(rikinv,delrik,rik_hat);
ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,p,trip);
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk,
PairParameters *p, TripletParameters *trip)
{
double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta;
double dcosdri[3],dcosdrj[3],dcosdrk[3];
cos_theta = vec3_dot(rij_hat,rik_hat);
(p->W)->value(rik,fc,1,dfc,1);
(p->P)->value(rij-(p->xi)*rik,ex_delr,1,ex_delr_d,1);
(trip->G)->value(cos_theta,gijk,1,gijk_d,1);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Ri
// dri = -dfc*gijk*ex_delr*rik_hat;
// dri += fc*gijk_d*ex_delr*dcosdri;
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
vec3_scale(prefactor,dri,dri);
// compute the derivative wrt Rj
// drj = fc*gijk_d*ex_delr*dcosdrj;
// drj += fc*gijk*ex_delr_d*rij_hat;
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
vec3_scale(prefactor,drj,drj);
// compute the derivative wrt Rk
// drk = dfc*gijk*ex_delr*rik_hat;
// drk += fc*gijk_d*ex_delr*dcosdrk;
// drk += -fc*gijk*ex_delr_d*rik_hat;
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
vec3_scale(prefactor,drk,drk);
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::costheta_d(double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk)
{
// first element is devative wrt Ri, second wrt Rj, third wrt Rk
double cos_theta = vec3_dot(rij_hat,rik_hat);
vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
vec3_scale(1.0/rij,drj,drj);
vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
vec3_scale(1.0/rik,drk,drk);
vec3_add(drj,drk,dri);
vec3_scale(-1.0,dri,dri);
}
/* ----------------------------------------------------------------------
* grab n values from file fp and put them in list
* values can be several to a line
* only called by proc 0
* ------------------------------------------------------------------------- */
void PairPolymorphic::grab(FILE *fp, int n, double *list)
{
char *ptr;
char line[MAXLINE];
int i = 0;
while (i < n) {
fgets(line,MAXLINE,fp);
ptr = strtok(line," \t\n\r\f");
list[i++] = atof(ptr);
while ((ptr = strtok(NULL," \t\n\r\f")))
list[i++] = atof(ptr);
}
}
/* ---------------------------------------------------------------------- */
void PairPolymorphic::write_tables(int npts)
{
char tag[6] = "";
if (comm->me != 0) sprintf(tag,"%d",comm->me);
FILE* fp = NULL;
double xmin,xmax,x,uf,vf,wf,pf,gf,ff,ufp,vfp,wfp,pfp,gfp,ffp;
char line[MAXLINE];
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_UVW");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.U)->get_xmin();
xmax = (pair.U)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.U)->value(x, uf, 1, ufp, 1);
(pair.V)->value(x, vf, 1, vfp, 1);
(pair.W)->value(x, wf, 1, wfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f \n",x,uf,vf,wf,ufp,vfp,wfp);
}
fclose(fp);
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_P");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.P)->get_xmin();
xmax = (pair.P)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.P)->value(x, pf, 1, pfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp);
}
fclose(fp);
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,elements[k]);
strcat(line,"_G");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.G)->get_xmin();
xmax = (pair.G)->get_xmax();
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.G)->value(x, gf, 1, gfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,gf,gfp);
}
fclose(fp);
}
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_F");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.F)->get_xmin();
xmax = (pair.F)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.F)->value(x, ff, 1, ffp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp);
}
fclose(fp);
}
}
}
Event Timeline
Log In to Comment