Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65195225
pair_multi_lucy.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
Sat, Jun 1, 18:02
Size
24 KB
Mime Type
text/x-c
Expires
Mon, Jun 3, 18:02 (2 d)
Engine
blob
Format
Raw Data
Handle
17995342
Attached To
rLAMMPS lammps
pair_multi_lucy.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 authors:
James Larentzos and Joshua Moore (U.S. Army Research Laboratory)
Please cite the related publications:
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include "math_const.h"
#include <stdlib.h>
#include <string.h>
#include "pair_multi_lucy.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
using namespace LAMMPS_NS;
enum{NONE,RLINEAR,RSQ};
#define MAXLINE 1024
static const char cite_pair_multi_lucy[] =
"pair_style multi/lucy command:\n\n"
"@Article{Moore16,\n"
" author = {J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor and J. K. Brennan},\n"
" title = {A coarse-grain force field for RDX: Density dependent and energy conserving},\n"
" journal = {J. Chem. Phys.},\n"
" year = 2016,\n"
" volume = 144\n"
" pages = {104501}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairMultiLucy::PairMultiLucy(LAMMPS *lmp) : Pair(lmp),
ntables(0), tables(NULL), tabindex(NULL)
{
if (lmp->citeme) lmp->citeme->add(cite_pair_multi_lucy);
if (atom->rho_flag != 1) error->all(FLERR,"Pair multi/lucy command requires atom_style with density (e.g. dpd, meso)");
ntables = 0;
tables = NULL;
comm_forward = 1;
comm_reverse = 1;
}
/* ---------------------------------------------------------------------- */
PairMultiLucy::~PairMultiLucy()
{
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
memory->sfree(tables);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(tabindex);
}
}
/* ---------------------------------------------------------------------- */
void PairMultiLucy::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype,itable;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
Table *tb;
int tlm1 = tablength - 1;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double pi = MathConst::MY_PI;
double A_i;
double A_j;
double fraction_i,fraction_j;
int jtable;
double *rho = atom->rho;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
computeLocalDensity();
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
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 &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
tb = &tables[tabindex[itype][jtype]];
if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq)
error->one(FLERR,"Density < table inner cutoff");
if (tabstyle == LOOKUP) {
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1)
error->one(FLERR,"Density > table outer cutoff");
A_i = tb->f[itable];
A_j = tb->f[jtable];
fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair/sqrt(rsq);
} else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1)
error->one(FLERR,"Density > table outer cutoff");
fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
fraction_j = (((rho[j]*rho[j]) - tb->rsq[jtable]) * tb->invdelta);
A_i = tb->f[itable] + fraction_i*tb->df[itable];
A_j = tb->f[jtable] + fraction_j*tb->df[jtable];
fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair / sqrt(rsq);
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy");
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 (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,0.0,fpair,delx,dely,delz);
}
}
tb = &tables[tabindex[itype][itype]];
if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq)
error->one(FLERR,"Density < table inner cutoff");
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
if (tabstyle == LOOKUP) evdwl = tb->e[itable];
else if (tabstyle == LINEAR){
if (itable >= tlm1) error->one(FLERR,"Density > table outer cutoff");
if(itable==0) fraction_i=0.0;
else fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
evdwl,0.0,0.0,0.0,0.0,0.0);
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairMultiLucy::allocate()
{
allocated = 1;
const int nt = atom->ntypes + 1;
memory->create(setflag,nt,nt,"pair:setflag");
memory->create(cutsq,nt,nt,"pair:cutsq");
memory->create(tabindex,nt,nt,"pair:tabindex");
memset(&setflag[0][0],0,nt*nt*sizeof(int));
memset(&cutsq[0][0],0,nt*nt*sizeof(double));
memset(&tabindex[0][0],0,nt*nt*sizeof(int));
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairMultiLucy::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal pair_style command");
// new settings
if (strcmp(arg[0],"lookup") == 0) tabstyle = LOOKUP;
else if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
else error->all(FLERR,"Unknown table style in pair_style command");
tablength = force->inumeric(FLERR,arg[1]);
if (tablength < 2) error->all(FLERR,"Illegal number of pair table entries");
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
memory->sfree(tables);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(tabindex);
}
allocated = 0;
ntables = 0;
tables = NULL;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairMultiLucy::coeff(int narg, char **arg)
{
if (narg != 4 && narg != 5) error->all(FLERR,"Illegal pair_coeff command");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
int me;
MPI_Comm_rank(world,&me);
tables = (Table *)
memory->srealloc(tables,(ntables+1)*sizeof(Table),"pair:tables");
Table *tb = &tables[ntables];
null_table(tb);
if (me == 0) read_table(tb,arg[2],arg[3]);
bcast_table(tb);
// set table cutoff
if (narg == 5) tb->cut = force->numeric(FLERR,arg[4]);
else if (tb->rflag) tb->cut = tb->rhi;
else tb->cut = tb->rfile[tb->ninput-1];
// error check on table parameters
// insure cutoff is within table
if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length");
double rlo;
if (tb->rflag == 0) {
rlo = tb->rfile[0];
} else {
rlo = tb->rlo;
}
rho_0 = rlo;
tb->match = 0;
if (tabstyle == LINEAR && tb->ninput == tablength &&
tb->rflag == RSQ) tb->match = 1;
// spline read-in values and compute r,e,f vectors within table
if (tb->match == 0) spline_table(tb);
compute_table(tb);
// store ptr to table in tabindex
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
tabindex[i][j] = ntables;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Illegal pair_coeff command");
ntables++;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairMultiLucy::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
tabindex[j][i] = tabindex[i][j];
return tables[tabindex[i][j]].cut;
}
/* ----------------------------------------------------------------------
read a table section from a tabulated potential file
only called by proc 0
this function sets these values in Table:
ninput,rfile,efile,ffile,rflag,rlo,rhi,fpflag,fplo,fphi
------------------------------------------------------------------------- */
void PairMultiLucy::read_table(Table *tb, char *file, char *keyword)
{
char line[MAXLINE];
// open file
FILE *fp = fopen(file,"r");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
// loop until section found with matching keyword
while (1) {
if (fgets(line,MAXLINE,fp) == NULL)
error->one(FLERR,"Did not find keyword in table file");
if (strspn(line," \t\n\r") == strlen(line)) continue; // blank line
if (line[0] == '#') continue; // comment
char *word = strtok(line," \t\n\r");
if (strcmp(word,keyword) == 0) break; // matching keyword
fgets(line,MAXLINE,fp); // no match, skip section
param_extract(tb,line);
fgets(line,MAXLINE,fp);
for (int i = 0; i < tb->ninput; i++) fgets(line,MAXLINE,fp);
}
// read args on 2nd line of section
// allocate table arrays for file values
fgets(line,MAXLINE,fp);
param_extract(tb,line);
memory->create(tb->rfile,tb->ninput,"pair:rfile");
memory->create(tb->efile,tb->ninput,"pair:efile");
memory->create(tb->ffile,tb->ninput,"pair:ffile");
// read r,e,f table values from file
// if rflag set, compute r
// if rflag not set, use r from file
int itmp;
double rtmp;
fgets(line,MAXLINE,fp);
for (int i = 0; i < tb->ninput; i++) {
fgets(line,MAXLINE,fp);
sscanf(line,"%d %lg %lg %lg",&itmp,&rtmp,&tb->efile[i],&tb->ffile[i]);
if (tb->rflag == RLINEAR)
rtmp = tb->rlo + (tb->rhi - tb->rlo)*i/(tb->ninput-1);
else if (tb->rflag == RSQ) {
rtmp = tb->rlo*tb->rlo +
(tb->rhi*tb->rhi - tb->rlo*tb->rlo)*i/(tb->ninput-1);
rtmp = sqrt(rtmp);
}
tb->rfile[i] = rtmp;
}
// close file
fclose(fp);
}
/* ----------------------------------------------------------------------
broadcast read-in table info from proc 0 to other procs
this function communicates these values in Table:
ninput,rfile,efile,ffile,rflag,rlo,rhi,fpflag,fplo,fphi
------------------------------------------------------------------------- */
void PairMultiLucy::bcast_table(Table *tb)
{
MPI_Bcast(&tb->ninput,1,MPI_INT,0,world);
int me;
MPI_Comm_rank(world,&me);
if (me > 0) {
memory->create(tb->rfile,tb->ninput,"pair:rfile");
memory->create(tb->efile,tb->ninput,"pair:efile");
memory->create(tb->ffile,tb->ninput,"pair:ffile");
}
MPI_Bcast(tb->rfile,tb->ninput,MPI_DOUBLE,0,world);
MPI_Bcast(tb->efile,tb->ninput,MPI_DOUBLE,0,world);
MPI_Bcast(tb->ffile,tb->ninput,MPI_DOUBLE,0,world);
MPI_Bcast(&tb->rflag,1,MPI_INT,0,world);
if (tb->rflag) {
MPI_Bcast(&tb->rlo,1,MPI_DOUBLE,0,world);
MPI_Bcast(&tb->rhi,1,MPI_DOUBLE,0,world);
}
MPI_Bcast(&tb->fpflag,1,MPI_INT,0,world);
if (tb->fpflag) {
MPI_Bcast(&tb->fplo,1,MPI_DOUBLE,0,world);
MPI_Bcast(&tb->fphi,1,MPI_DOUBLE,0,world);
}
}
/* ----------------------------------------------------------------------
build spline representation of e,f over entire range of read-in table
this function sets these values in Table: e2file,f2file
------------------------------------------------------------------------- */
void PairMultiLucy::spline_table(Table *tb)
{
memory->create(tb->e2file,tb->ninput,"pair:e2file");
memory->create(tb->f2file,tb->ninput,"pair:f2file");
double ep0 = - tb->ffile[0];
double epn = - tb->ffile[tb->ninput-1];
spline(tb->rfile,tb->efile,tb->ninput,ep0,epn,tb->e2file);
if (tb->fpflag == 0) {
tb->fplo = (tb->ffile[1] - tb->ffile[0]) / (tb->rfile[1] - tb->rfile[0]);
tb->fphi = (tb->ffile[tb->ninput-1] - tb->ffile[tb->ninput-2]) /
(tb->rfile[tb->ninput-1] - tb->rfile[tb->ninput-2]);
}
double fp0 = tb->fplo;
double fpn = tb->fphi;
spline(tb->rfile,tb->ffile,tb->ninput,fp0,fpn,tb->f2file);
}
/* ----------------------------------------------------------------------
extract attributes from parameter line in table section
format of line: N value R/RSQ lo hi FP fplo fphi
N is required, other params are optional
------------------------------------------------------------------------- */
void PairMultiLucy::param_extract(Table *tb, char *line)
{
tb->ninput = 0;
tb->rflag = NONE;
tb->fpflag = 0;
char *word = strtok(line," \t\n\r\f");
while (word) {
if (strcmp(word,"N") == 0) {
word = strtok(NULL," \t\n\r\f");
tb->ninput = atoi(word);
} else if (strcmp(word,"R") == 0 || strcmp(word,"RSQ") == 0) {
if (strcmp(word,"R") == 0) tb->rflag = RLINEAR;
else if (strcmp(word,"RSQ") == 0) tb->rflag = RSQ;
word = strtok(NULL," \t\n\r\f");
tb->rlo = atof(word);
word = strtok(NULL," \t\n\r\f");
tb->rhi = atof(word);
} else if (strcmp(word,"FP") == 0) {
tb->fpflag = 1;
word = strtok(NULL," \t\n\r\f");
tb->fplo = atof(word);
word = strtok(NULL," \t\n\r\f");
tb->fphi = atof(word);
} else {
printf("WORD: %s\n",word);
error->one(FLERR,"Invalid keyword in pair table parameters");
}
word = strtok(NULL," \t\n\r\f");
}
if (tb->ninput == 0) error->one(FLERR,"Pair table parameters did not set N");
}
/* ----------------------------------------------------------------------
compute r,e,f vectors from splined values
------------------------------------------------------------------------- */
void PairMultiLucy::compute_table(Table *tb)
{
int tlm1 = tablength-1;
// inner = inner table bound
// cut = outer table bound
// delta = table spacing in rsq for N-1 bins
double inner;
if (tb->rflag) inner = tb->rlo;
else inner = tb->rfile[0];
tb->innersq = inner*inner;
tb->delta = (tb->rhi*tb->rhi - tb->innersq) / tlm1;
tb->invdelta = 1.0/tb->delta;
// direct lookup tables
// N-1 evenly spaced bins in rsq from inner to cut
// e,f = value at midpt of bin
// e,f are N-1 in length since store 1 value at bin midpt
// f is converted to f/r when stored in f[i]
// e,f are never a match to read-in values, always computed via spline interp
if (tabstyle == LOOKUP) {
memory->create(tb->e,tlm1,"pair:e");
memory->create(tb->f,tlm1,"pair:f");
double r,rsq;
for (int i = 0; i < tlm1; i++) {
rsq = tb->innersq + (i+0.5)*tb->delta;
r = sqrt(rsq);
tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r);
tb->f[i] = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,r);
}
}
// linear tables
// N-1 evenly spaced bins in rsq from inner to cut
// rsq,e,f = value at lower edge of bin
// de,df values = delta from lower edge to upper edge of bin
// rsq,e,f are N in length so de,df arrays can compute difference
// f is converted to f/r when stored in f[i]
// e,f can match read-in values, else compute via spline interp
if (tabstyle == LINEAR) {
memory->create(tb->rsq,tablength,"pair:rsq");
memory->create(tb->e,tablength,"pair:e");
memory->create(tb->f,tablength,"pair:f");
memory->create(tb->de,tlm1,"pair:de");
memory->create(tb->df,tlm1,"pair:df");
double r,rsq;
for (int i = 0; i < tablength; i++) {
rsq = tb->innersq + i*tb->delta;
r = sqrt(rsq);
tb->rsq[i] = rsq;
if (tb->match) {
tb->e[i] = tb->efile[i];
tb->f[i] = tb->ffile[i];
} else {
tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r);
tb->f[i] = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,r);
}
}
for (int i = 0; i < tlm1; i++) {
tb->de[i] = tb->e[i+1] - tb->e[i];
tb->df[i] = tb->f[i+1] - tb->f[i];
}
}
}
/* ----------------------------------------------------------------------
set all ptrs in a table to NULL, so can be freed safely
------------------------------------------------------------------------- */
void PairMultiLucy::null_table(Table *tb)
{
tb->rfile = tb->efile = tb->ffile = NULL;
tb->e2file = tb->f2file = NULL;
tb->rsq = tb->drsq = tb->e = tb->de = NULL;
tb->f = tb->df = tb->e2 = tb->f2 = NULL;
}
/* ----------------------------------------------------------------------
free all arrays in a table
------------------------------------------------------------------------- */
void PairMultiLucy::free_table(Table *tb)
{
memory->destroy(tb->rfile);
memory->destroy(tb->efile);
memory->destroy(tb->ffile);
memory->destroy(tb->e2file);
memory->destroy(tb->f2file);
memory->destroy(tb->rsq);
memory->destroy(tb->drsq);
memory->destroy(tb->e);
memory->destroy(tb->de);
memory->destroy(tb->f);
memory->destroy(tb->df);
memory->destroy(tb->e2);
memory->destroy(tb->f2);
}
/* ----------------------------------------------------------------------
spline and splint routines modified from Numerical Recipes
------------------------------------------------------------------------- */
void PairMultiLucy::spline(double *x, double *y, int n,
double yp1, double ypn, double *y2)
{
int i,k;
double p,qn,sig,un;
double *u = new double[n];
if (yp1 > 0.99e30) y2[0] = u[0] = 0.0;
else {
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
}
for (i = 1; i < n-1; i++) {
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p = sig*y2[i-1] + 2.0;
y2[i] = (sig-1.0) / p;
u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
}
if (ypn > 0.99e30) qn = un = 0.0;
else {
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2] + 1.0);
for (k = n-2; k >= 0; k--) y2[k] = y2[k]*y2[k+1] + u[k];
delete [] u;
}
/* ---------------------------------------------------------------------- */
double PairMultiLucy::splint(double *xa, double *ya, double *y2a, int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while (khi-klo > 1) {
k = (khi+klo) >> 1;
if (xa[k] > x) khi = k;
else klo = k;
}
h = xa[khi]-xa[klo];
a = (xa[khi]-x) / h;
b = (x-xa[klo]) / h;
y = a*ya[klo] + b*ya[khi] +
((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return y;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairMultiLucy::write_restart(FILE *fp)
{
write_restart_settings(fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairMultiLucy::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairMultiLucy::write_restart_settings(FILE *fp)
{
fwrite(&tabstyle,sizeof(int),1,fp);
fwrite(&tablength,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairMultiLucy::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&tabstyle,sizeof(int),1,fp);
fread(&tablength,sizeof(int),1,fp);
}
MPI_Bcast(&tabstyle,1,MPI_INT,0,world);
MPI_Bcast(&tablength,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
void PairMultiLucy::computeLocalDensity()
{
int i,j,m,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
double **x = atom->x;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double factor;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double pi = MathConst::MY_PI;
double *rho = atom->rho;
// zero out density
if (newton_pair) {
m = nlocal + atom->nghost;
for (i = 0; i < m; i++) rho[i] = 0.0;
} else for (i = 0; i < nlocal; i++) rho[i] = 0.0;
// rho = density at each atom
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
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 &= NEIGHMASK;
jtype = 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 < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;
}
}
}
}
if (newton_pair) comm->reverse_comm_pair(this);
comm->forward_comm_pair(this);
}
/* ---------------------------------------------------------------------- */
int PairMultiLucy::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
double *rho = atom->rho;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = rho[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairMultiLucy::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
double *rho = atom->rho;
m = 0;
last = first + n;
for (i = first; i < last; i++) rho[i] = buf[m++];
}
/* ---------------------------------------------------------------------- */
int PairMultiLucy::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
double *rho = atom->rho;
m = 0;
last = first + n;
for (i = first; i < last; i++) buf[m++] = rho[i];
return m;
}
/* ---------------------------------------------------------------------- */
void PairMultiLucy::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
double *rho = atom->rho;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
}
}
Event Timeline
Log In to Comment