Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84067109
msm.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
Fri, Sep 20, 13:20
Size
75 KB
Mime Type
text/x-c
Expires
Sun, Sep 22, 13:20 (2 d)
Engine
blob
Format
Raw Data
Handle
20935391
Attached To
rLAMMPS lammps
msm.cpp
View Options
This document is not UTF8. It was detected as ISO-8859-1 (Latin 1) and converted to UTF8 for display.
/* ----------------------------------------------------------------------
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: Paul Crozier, Stan Moore, Stephen Bond, (all SNL)
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "string.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "msm.h"
#include "atom.h"
#include "comm.h"
#include "commgrid.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAX_LEVELS 10
#define OFFSET 16384
#define SMALL 0.00001
#define LARGE 10000.0
enum{REVERSE_RHO,REVERSE_AD,REVERSE_AD_PERATOM};
enum{FORWARD_RHO,FORWARD_AD,FORWARD_AD_PERATOM};
/* ---------------------------------------------------------------------- */
MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
{
if (narg < 1) error->all(FLERR,"Illegal kspace_style msm command");
msmflag = 1;
accuracy_relative = atof(arg[0]);
nfactors = 1;
factors = new int[nfactors];
factors[0] = 2;
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
phi1d = dphi1d = NULL;
nmax = 0;
part2grid = NULL;
g_direct = NULL;
g_direct_top = NULL;
v0_direct = v1_direct = v2_direct = NULL;
v3_direct = v4_direct = v5_direct = NULL;
v0_direct_top = v1_direct_top = v2_direct_top = NULL;
v3_direct_top = v4_direct_top = v5_direct_top = NULL;
cg = cg_peratom = NULL;
levels = 0;
peratom_allocate_flag = 0;
order = 8;
}
/* ----------------------------------------------------------------------
free all memory
------------------------------------------------------------------------- */
MSM::~MSM()
{
delete [] factors;
deallocate();
deallocate_peratom();
memory->destroy(part2grid);
memory->destroy(g_direct);
memory->destroy(g_direct_top);
memory->destroy(v0_direct);
memory->destroy(v1_direct);
memory->destroy(v2_direct);
memory->destroy(v3_direct);
memory->destroy(v4_direct);
memory->destroy(v5_direct);
memory->destroy(v0_direct_top);
memory->destroy(v1_direct_top);
memory->destroy(v2_direct_top);
memory->destroy(v3_direct_top);
memory->destroy(v4_direct_top);
memory->destroy(v5_direct_top);
deallocate_levels();
}
/* ----------------------------------------------------------------------
called once before run
------------------------------------------------------------------------- */
void MSM::init()
{
if (me == 0) {
if (screen) fprintf(screen,"MSM initialization ...\n");
if (logfile) fprintf(logfile,"MSM initialization ...\n");
}
// error check
if (domain->triclinic)
error->all(FLERR,"Cannot (yet) use MSM with triclinic box");
if (domain->dimension == 2)
error->all(FLERR,"Cannot (yet) use MSM with 2d simulation");
if (!atom->q_flag) error->all(FLERR,"Kspace style requires atom attribute q");
if (slabflag == 1)
error->all(FLERR,"Cannot use slab correction with MSM");
if (order < 4 || order > 10) {
char str[128];
sprintf(str,"MSM order must be 4, 6, 8, or 10");
error->all(FLERR,str);
}
if (order%2 != 0) error->all(FLERR,"MSM order must be 4, 6, 8, or 10");
if (sizeof(FFT_SCALAR) != 8) error->all(FLERR,"Cannot (yet) use single precision with MSM (remove -DFFT_SINGLE from Makefile and recompile)");
// extract short-range Coulombic cutoff from pair style
qqrd2e = force->qqrd2e;
scale = 1.0;
pair_check();
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_msm",itmp);
if (p_cutoff == NULL)
error->all(FLERR,"KSpace style is incompatible with Pair style");
cutoff = *p_cutoff;
// compute qsum & qsqsum and give error if not charge-neutral
qsum = qsqsum = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
qsum += atom->q[i];
qsqsum += atom->q[i]*atom->q[i];
}
double tmp;
MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum = tmp;
MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsqsum = tmp;
q2 = qsqsum * force->qqrd2e / force->dielectric;
if (qsqsum == 0.0)
error->all(FLERR,"Cannot use kspace solver on system with no charge");
// not yet sure of the correction needed for non-neutral systems
if (fabs(qsum) > SMALL) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->all(FLERR,str);
}
// set accuracy (force units) from accuracy_relative or accuracy_absolute
if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute;
else accuracy = accuracy_relative * two_charge_force;
// setup MSM grid resolution
set_grid_global();
setup();
double estimated_error = estimate_total_error();
// output grid stats
int ngrid_max;
MPI_Allreduce(&ngrid[0],&ngrid_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
if (screen) {
fprintf(screen," 3d grid size/proc = %d\n",
ngrid_max);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_error);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_error/two_charge_force);
}
if (logfile) {
fprintf(logfile," 3d grid size/proc = %d\n",
ngrid_max);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_error);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_error/two_charge_force);
}
}
if (me == 0) {
if (screen) {
fprintf(screen," grid = %d %d %d\n",nx_msm[0],ny_msm[0],nz_msm[0]);
fprintf(screen," order = %d\n",order);
}
if (logfile) {
fprintf(logfile," grid = %d %d %d\n",nx_msm[0],ny_msm[0],nz_msm[0]);
fprintf(logfile," order = %d\n",order);
}
}
}
/* ----------------------------------------------------------------------
estimate cutoff for a given grid spacing and error
------------------------------------------------------------------------- */
double MSM::estimate_cutoff(double h, double prd)
{
double a;
int p = order - 1;
double Mp,cprime,error_scaling;
Mp = cprime = error_scaling = 1;
// Mp values from Table 5.1 of Hardy's thesis
// cprime values from equation 4.17 of Hardy's thesis
// error scaling from empirical fitting to convert to rms force errors
if (p == 3) {
Mp = 9;
cprime = 1.0/6.0;
error_scaling = 0.39189561;
} else if (p == 5) {
Mp = 825;
cprime = 1.0/30.0;
error_scaling = 0.150829428;
} else if (p == 7) {
Mp = 130095;
cprime = 1.0/140.0;
error_scaling = 0.049632967;
} else if (p == 9) {
Mp = 34096545;
cprime = 1.0/630.0;
error_scaling = 0.013520855;
} else {
error->all(FLERR,"MSM order must be 4, 6, 8, or 10");
}
// equation 4.1 from Hardy's thesis
double C_p = 4.0*cprime*Mp/3.0;
// use empirical parameters to convert to rms force errors
C_p *= error_scaling;
// equation 3.200 from Hardy's thesis
a = C_p*pow(h,(p-1))/accuracy;
// include dependency of error on other terms
a *= q2/(prd*sqrt(atom->natoms));
a = pow(a,1.0/double(p));
return a;
}
/* ----------------------------------------------------------------------
estimate 1d grid RMS force error for MSM
------------------------------------------------------------------------- */
double MSM::estimate_1d_error(double h, double prd)
{
double a = cutoff;
int p = order - 1;
double Mp,cprime,error_scaling;
Mp = cprime = error_scaling = 1;
// Mp values from Table 5.1 of Hardy's thesis
// cprime values from equation 4.17 of Hardy's thesis
// error scaling from empirical fitting to convert to rms force errors
if (p == 3) {
Mp = 9;
cprime = 1.0/6.0;
error_scaling = 0.39189561;
} else if (p == 5) {
Mp = 825;
cprime = 1.0/30.0;
error_scaling = 0.150829428;
} else if (p == 7) {
Mp = 130095;
cprime = 1.0/140.0;
error_scaling = 0.049632967;
} else if (p == 9) {
Mp = 34096545;
cprime = 1.0/630.0;
error_scaling = 0.013520855;
} else {
error->all(FLERR,"MSM order must be 4, 6, 8, or 10");
}
// equation 4.1 from Hardy's thesis
double C_p = 4.0*cprime*Mp/3.0;
// use empirical parameters to convert to rms force errors
C_p *= error_scaling;
// equation 3.197 from Hardy's thesis
double error_1d = C_p*pow(h,(p-1))/pow(a,(p+1));
// include dependency of error on other terms
error_1d *= q2*a/(prd*sqrt(atom->natoms));
return error_1d;
}
/* ----------------------------------------------------------------------
estimate 3d grid RMS force error
------------------------------------------------------------------------- */
double MSM::estimate_3d_error()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double error_x = estimate_1d_error(xprd/nx_msm[0],xprd);
double error_y = estimate_1d_error(yprd/ny_msm[0],yprd);
double error_z = estimate_1d_error(zprd/nz_msm[0],zprd);
double error_3d =
sqrt(error_x*error_x + error_y*error_y + error_z*error_z) / sqrt(3.0);
return error_3d;
}
/* ----------------------------------------------------------------------
estimate total RMS force error
------------------------------------------------------------------------- */
double MSM::estimate_total_error()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
bigint natoms = atom->natoms;
double grid_error = estimate_3d_error();
double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd);
double short_range_error = 0.0;
double table_error =
estimate_table_accuracy(q2_over_sqrt,short_range_error);
double estimated_total_error = sqrt(grid_error*grid_error +
short_range_error*short_range_error + table_error*table_error);
return estimated_total_error;
}
/* ----------------------------------------------------------------------
adjust MSM coeffs, called initially and whenever volume has changed
------------------------------------------------------------------------- */
void MSM::setup()
{
int i,j,k,l,m,n;
double *prd;
double a = cutoff;
// volume-dependent factors
prd = domain->prd;
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
volume = xprd * yprd * zprd;
// loop over grid levels
for (int n=0; n<levels; n++) {
delxinv[n] = nx_msm[n]/xprd;
delyinv[n] = ny_msm[n]/yprd;
delzinv[n] = nz_msm[n]/zprd;
delvolinv[n] = delxinv[n]*delyinv[n]*delzinv[n];
}
nxhi_direct = static_cast<int> (2.0*a*delxinv[0]);
nxlo_direct = -nxhi_direct;
nyhi_direct = static_cast<int> (2.0*a*delyinv[0]);
nylo_direct = -nyhi_direct;
nzhi_direct = static_cast<int> (2.0*a*delzinv[0]);
nzlo_direct = -nzhi_direct;
nmax_direct = 8*(nxhi_direct+1)*(nyhi_direct+1)*(nzhi_direct+1);
deallocate();
deallocate_peratom();
if (!peratom_allocate_flag) { // Timestep 0
get_g_direct();
get_virial_direct();
if (domain->nonperiodic) {
get_g_direct_top(levels-1);
get_virial_direct_top(levels-1);
}
} else {
get_g_direct();
if (domain->nonperiodic) get_g_direct_top(levels-1);
if (vflag_either) {
get_virial_direct();
if (domain->nonperiodic) get_virial_direct_top(levels-1);
}
}
boxlo = domain->boxlo;
set_grid_local();
// allocate K-space dependent memory
// don't invoke allocate_peratom(), compute() will allocate when needed
allocate();
peratom_allocate_flag = 0;
for (int n=0; n<levels; n++) {
cg[n]->ghost_notify();
cg[n]->setup();
}
}
/* ----------------------------------------------------------------------
compute the MSM long-range force, energy, virial
------------------------------------------------------------------------- */
void MSM::compute(int eflag, int vflag)
{
int i,j;
// set energy/virial flags
// invoke allocate_peratom() if needed for first time
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = eflag_either = vflag_either = 0;
if (vflag_atom && !peratom_allocate_flag) {
allocate_peratom();
for (int n=0; n<levels; n++) {
cg_peratom[n]->ghost_notify();
cg_peratom[n]->setup();
}
peratom_allocate_flag = 1;
}
// extend size of per-atom arrays if necessary
if (atom->nlocal > nmax) {
memory->destroy(part2grid);
nmax = atom->nmax;
memory->create(part2grid,nmax,3,"msm:part2grid");
}
// find grid points for all my particles
// map my particle charge onto my local 3d density grid (aninterpolation)
particle_map();
make_rho();
current_level = 0;
cg[0]->reverse_comm(this,REVERSE_RHO);
// all procs communicate density values from their ghost cells
// to fully sum contribution in their 3d bricks
for (int n=0; n<=levels-2; n++) {
current_level = n;
cg[n]->forward_comm(this,FORWARD_RHO);
direct(n);
if (vflag_atom) direct_peratom(n);
restriction(n);
}
// top grid level
if (domain->nonperiodic) {
direct_top(levels-1);
if (vflag_atom) direct_peratom_top(levels-1);
}
for (int n=levels-2; n>=0; n--) {
prolongation(n);
current_level = n;
cg[n]->reverse_comm(this,REVERSE_AD);
// extra per-atom virial communication
if (vflag_atom)
cg_peratom[n]->reverse_comm(this,REVERSE_AD_PERATOM);
}
// all procs communicate E-field values
// to fill ghost cells surrounding their 3d bricks
current_level = 0;
cg[0]->forward_comm(this,FORWARD_AD);
// extra per-atom energy/virial communication
if (vflag_atom)
cg_peratom[0]->forward_comm(this,FORWARD_AD_PERATOM);
// calculate the force on my particles (interpolation)
fieldforce();
// calculate the per-atom energy for my particles
if (evflag_atom) fieldforce_peratom();
const double qscale = force->qqrd2e * scale;
// Total long-range energy
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
double e_self = qsqsum*gamma(0.0)/cutoff; // Self-energy term
energy -= e_self;
energy *= 0.5*qscale;
}
// Total long-range virial
if (vflag_global) {
double virial_all[6];
MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*virial_all[i];
}
// per-atom energy/virial
// energy includes self-energy correction
if (evflag_atom) {
double *q = atom->q;
int nlocal = atom->nlocal;
if (eflag_atom) {
for (i = 0; i < nlocal; i++) {
eatom[i] -= q[i]*q[i]*gamma(0.0)/cutoff;
eatom[i] *= 0.5*qscale;
}
}
if (vflag_atom) {
for (i = 0; i < nlocal; i++)
for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*qscale;
}
}
}
/* ----------------------------------------------------------------------
allocate memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::allocate()
{
// summation coeffs
memory->create2d_offset(phi1d,3,-order,order,"msm:phi1d");
memory->create2d_offset(dphi1d,3,-order,order,"msm:dphi1d");
// allocate grid levels
for (int n=0; n<levels; n++) {
memory->create3d_offset(qgrid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:qgrid");
memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid");
// create ghost grid object for rho and electric field communication
int (*procneigh)[2] = comm->procneigh;
cg[n] = new CommGrid(lmp,world,1,1,
nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n],nzlo_in[n],nzhi_in[n],
nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n],nzlo_out[n],nzhi_out[n],
procneigh[0][0],procneigh[0][1],procneigh[1][0],
procneigh[1][1],procneigh[2][0],procneigh[2][1]);
}
}
/* ----------------------------------------------------------------------
allocate per-atom memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::allocate_peratom()
{
// allocate grid levels
for (int n=0; n<levels; n++) {
memory->create3d_offset(v0grid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:v0grid");
memory->create3d_offset(v1grid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:v1grid");
memory->create3d_offset(v2grid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:v2grid");
memory->create3d_offset(v3grid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:v3grid");
memory->create3d_offset(v4grid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:v4grid");
memory->create3d_offset(v5grid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:v5grid");
// create ghost grid object for per-atom energy/virial
int (*procneigh)[2] = comm->procneigh;
cg_peratom[n] =
new CommGrid(lmp,world,6,6,
nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n],nzlo_in[n],nzhi_in[n],
nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n],nzlo_out[n],nzhi_out[n],
procneigh[0][0],procneigh[0][1],procneigh[1][0],
procneigh[1][1],procneigh[2][0],procneigh[2][1]);
}
}
/* ----------------------------------------------------------------------
deallocate memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::deallocate()
{
memory->destroy2d_offset(phi1d,-order);
memory->destroy2d_offset(dphi1d,-order);
// deallocate grid levels
for (int n=0; n<levels; n++) {
if (qgrid[n])
memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (egrid[n])
memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (cg)
if (cg[n]) delete cg[n];
}
}
/* ----------------------------------------------------------------------
deallocate per-atom memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::deallocate_peratom()
{
for (int n=0; n<levels; n++) {
if (v0grid[n])
memory->destroy3d_offset(v0grid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (v1grid[n])
memory->destroy3d_offset(v1grid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (v2grid[n])
memory->destroy3d_offset(v2grid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (v3grid[n])
memory->destroy3d_offset(v3grid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (v4grid[n])
memory->destroy3d_offset(v4grid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (v5grid[n])
memory->destroy3d_offset(v5grid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (cg_peratom)
if (cg_peratom[n]) delete cg_peratom[n];
}
}
/* ----------------------------------------------------------------------
allocate memory that depends on # of grid levels
------------------------------------------------------------------------- */
void MSM::allocate_levels()
{
ngrid = new int[levels];
cg = new CommGrid*[levels];
cg_peratom = new CommGrid*[levels];
alpha = new int[levels];
betax = new int[levels];
betay = new int[levels];
betaz = new int[levels];
nx_msm = new int[levels];
ny_msm = new int[levels];
nz_msm = new int[levels];
nxlo_in = new int[levels];
nylo_in = new int[levels];
nzlo_in = new int[levels];
nxhi_in = new int[levels];
nyhi_in = new int[levels];
nzhi_in = new int[levels];
nxlo_out = new int[levels];
nylo_out = new int[levels];
nzlo_out = new int[levels];
nxhi_out = new int[levels];
nyhi_out = new int[levels];
nzhi_out = new int[levels];
delxinv = new double[levels];
delyinv = new double[levels];
delzinv = new double[levels];
delvolinv = new double[levels];
qgrid = new double***[levels];
egrid = new double***[levels];
v0grid = new double***[levels];
v1grid = new double***[levels];
v2grid = new double***[levels];
v3grid = new double***[levels];
v4grid = new double***[levels];
v5grid = new double***[levels];
for (int n=0; n<levels; n++) {
cg[n] = NULL;
cg_peratom[n] = NULL;
qgrid[n] = NULL;
egrid[n] = NULL;
v0grid[n] = NULL;
v1grid[n] = NULL;
v2grid[n] = NULL;
v3grid[n] = NULL;
v4grid[n] = NULL;
v5grid[n] = NULL;
}
}
/* ----------------------------------------------------------------------
deallocate memory that depends on # of grid levels
------------------------------------------------------------------------- */
void MSM::deallocate_levels()
{
delete [] ngrid;
delete [] cg;
delete [] cg_peratom;
delete [] alpha;
delete [] betax;
delete [] betay;
delete [] betaz;
delete [] nx_msm;
delete [] ny_msm;
delete [] nz_msm;
delete [] nxlo_in;
delete [] nylo_in;
delete [] nzlo_in;
delete [] nxhi_in;
delete [] nyhi_in;
delete [] nzhi_in;
delete [] nxlo_out;
delete [] nylo_out;
delete [] nzlo_out;
delete [] nxhi_out;
delete [] nyhi_out;
delete [] nzhi_out;
delete [] delxinv;
delete [] delyinv;
delete [] delzinv;
delete [] delvolinv;
delete [] qgrid;
delete [] egrid;
delete [] v0grid;
delete [] v1grid;
delete [] v2grid;
delete [] v3grid;
delete [] v4grid;
delete [] v5grid;
}
/* ----------------------------------------------------------------------
set size of MSM grids
------------------------------------------------------------------------- */
void MSM::set_grid_global()
{
if (accuracy_relative <= 0.0)
error->all(FLERR,"KSpace accuracy must be > 0");
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int nx_max,ny_max,nz_max;
double hx,hy,hz;
if (adjust_cutoff_flag && !gridflag) {
int p = order - 1;
double hmin = 3072.0*(p+1)/(p-1)/
(448.0*MY_PI + 56.0*MY_PI*order/2 + 1701.0);
hmin = pow(hmin,1.0/6.0)/pow(atom->natoms,1.0/3.0);
hx = hmin*xprd;
hy = hmin*yprd;
hz = hmin*zprd;
nx_max = static_cast<int>(xprd/hx);
ny_max = static_cast<int>(yprd/hy);
nz_max = static_cast<int>(zprd/hz);
} else if (!gridflag) {
nx_max = ny_max = nz_max = 2;
hx = xprd/nx_max;
hy = yprd/ny_max;
hz = zprd/nz_max;
double x_error = 2.0*accuracy;
double y_error = 2.0*accuracy;
double z_error = 2.0*accuracy;
while (x_error > accuracy) {
nx_max *= 2;
hx = xprd/nx_max;
x_error = estimate_1d_error(hx,xprd);
}
while (y_error > accuracy) {
ny_max *= 2;
hy = yprd/ny_max;
y_error = estimate_1d_error(hy,yprd);
}
while (z_error > accuracy) {
nz_max *= 2;
hz = zprd/nz_max;
z_error = estimate_1d_error(hz,zprd);
}
} else {
nx_max = nx_msm_max;
ny_max = ny_msm_max;
nz_max = nz_msm_max;
}
// boost grid size until it is factorable
int flag = 0;
int xlevels,ylevels,zlevels;
while (!factorable(nx_max,flag,xlevels)) nx_max++;
while (!factorable(ny_max,flag,ylevels)) ny_max++;
while (!factorable(nz_max,flag,zlevels)) nz_max++;
if (flag && gridflag && me == 0)
error->warning(FLERR,"Number of MSM mesh points increased to be a multiple of 2");
// Find maximum number of levels
levels = MAX(xlevels,ylevels);
levels = MAX(levels,zlevels);
if (levels > MAX_LEVELS)
error->all(FLERR,"Too many MSM grid levels");
// Need at least 2 MSM levels for periodic systems
if (levels <= 1) {
levels = xlevels = ylevels = zlevels = 2;
nx_max = ny_max = nz_max = 2;
if (gridflag)
error->warning(FLERR,"MSM mesh too small, increasing to 2 points in each direction)");
}
if (adjust_cutoff_flag) {
hx = xprd/nx_max;
hy = yprd/ny_max;
hz = zprd/nz_max;
double ax,ay,az;
ax = estimate_cutoff(hx,xprd);
ay = estimate_cutoff(hy,yprd);
az = estimate_cutoff(hz,zprd);
cutoff = sqrt(ax*ax + ay*ay + az*az)/sqrt(3.0);
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_msm",itmp);
*p_cutoff = cutoff;
char str[128];
sprintf(str,"Adjusting Coulombic cutoff for MSM, new cutoff = %g",cutoff);
if (me == 0) error->warning(FLERR,str);
}
allocate_levels();
for (int n = 0; n < levels; n++) {
if (xlevels-n-1 > 0)
nx_msm[n] = static_cast<int> (pow(2.0,xlevels-n-1));
else
nx_msm[n] = 1;
if (ylevels-n-1 > 0)
ny_msm[n] = static_cast<int> (pow(2.0,ylevels-n-1));
else
ny_msm[n] = 1;
if (zlevels-n-1 > 0)
nz_msm[n] = static_cast<int> (pow(2.0,zlevels-n-1));
else
nz_msm[n] = 1;
}
if (nx_msm[0] >= OFFSET || ny_msm[0] >= OFFSET || nz_msm[0] >= OFFSET)
error->all(FLERR,"MSM grid is too large");
if (domain->nonperiodic) {
alpha[0] = -(order/2 - 1);
betax[0] = nx_msm[0] + (order/2 - 1);
betay[0] = ny_msm[0] + (order/2 - 1);
betaz[0] = nz_msm[0] + (order/2 - 1);
for (int n = 1; n < levels; n++) {
alpha[n] = -((-alpha[n-1]+1)/2) - (order/2 - 1);
betax[n] = ((betax[n-1]+1)/2) + (order/2 - 1);
betay[n] = ((betay[n-1]+1)/2) + (order/2 - 1);
betaz[n] = ((betaz[n-1]+1)/2) + (order/2 - 1);
}
}
}
/* ----------------------------------------------------------------------
set local subset of MSM grid that I own
n xyz lo/hi in = 3d brick that I own (inclusive)
n xyz lo/hi out = 3d brick + ghost cells in 6 directions (inclusive)
------------------------------------------------------------------------- */
void MSM::set_grid_local()
{
// global indices of MSM grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global MSM grid that I own without ghost cells
// loop over grid levels
for (int n=0; n<levels; n++) {
// global indices of MSM grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global MSM grid that I own without ghost cells
nxlo_in[n] = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx_msm[n]);
nxhi_in[n] = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx_msm[n]) - 1;
nylo_in[n] = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny_msm[n]);
nyhi_in[n] = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny_msm[n]) - 1;
nzlo_in[n] = static_cast<int> (comm->zsplit[comm->myloc[2]] * nz_msm[n]);
nzhi_in[n] = static_cast<int> (comm->zsplit[comm->myloc[2]+1] * nz_msm[n]) - 1;
// nlower,nupper = stencil size for mapping particles to MSM grid
nlower = -(order-1)/2;
nupper = order/2;
double *prd,*sublo,*subhi,*boxhi;
prd = domain->prd;
boxlo = domain->boxlo;
boxhi = domain->boxhi;
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
// shift values for particle <-> grid mapping
// add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
// nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of
// global MSM grid that my particles can contribute charge to
// effectively nlo_in,nhi_in + ghost cells
// nlo,nhi = global coords of grid pt to "lower left" of smallest/largest
// position a particle in my box can be at
// dist[3] = particle position bound = subbox + skin/2.0
// nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping
sublo = domain->sublo;
subhi = domain->subhi;
double dist[3];
double cuthalf = 0.0;
if (n == 0) cuthalf = 0.5*neighbor->skin; // Only applies to finest grid
dist[0] = dist[1] = dist[2] = cuthalf;
int nlo,nhi;
nlo = static_cast<int> ((sublo[0]-dist[0]-boxlo[0]) *
nx_msm[n]/xprd + OFFSET) - OFFSET;
nhi = static_cast<int> ((subhi[0]+dist[0]-boxlo[0]) *
nx_msm[n]/xprd + OFFSET) - OFFSET;
nxlo_out[n] = nlo + MIN(-order,nxlo_direct);
nxhi_out[n] = nhi + MAX(order,nxhi_direct);
nlo = static_cast<int> ((sublo[1]-dist[1]-boxlo[1]) *
ny_msm[n]/yprd + OFFSET) - OFFSET;
nhi = static_cast<int> ((subhi[1]+dist[1]-boxlo[1]) *
ny_msm[n]/yprd + OFFSET) - OFFSET;
nylo_out[n] = nlo + MIN(-order,nylo_direct);
nyhi_out[n] = nhi + MAX(order,nyhi_direct);
nlo = static_cast<int> ((sublo[2]-dist[2]-boxlo[2]) *
nz_msm[n]/zprd + OFFSET) - OFFSET;
nhi = static_cast<int> ((subhi[2]+dist[2]-boxlo[2]) *
nz_msm[n]/zprd + OFFSET) - OFFSET;
nzlo_out[n] = nlo + MIN(-order,nzlo_direct);
nzhi_out[n] = nhi + MAX(order,nzhi_direct);
// Add extra grid points for nonperiodic boundary conditions
if (domain->nonperiodic) {
if (!domain->xperiodic) {
if (nxlo_in[n] == 0)
nxlo_in[n] = alpha[n];
nxlo_out[n] = MAX(nxlo_out[n],alpha[n]);
if (nxhi_in[n] == nx_msm[n] - 1)
nxhi_in[n] = betax[n];
nxhi_out[n] = MIN(nxhi_out[n],betax[n]);
if (nxhi_in[n] < 0)
nxhi_in[n] = alpha[n] - 1;
}
if (!domain->yperiodic) {
if (nylo_in[n] == 0)
nylo_in[n] = alpha[n];
nylo_out[n] = MAX(nylo_out[n],alpha[n]);
if (nyhi_in[n] == ny_msm[n] - 1)
nyhi_in[n] = betay[n];
nyhi_out[n] = MIN(nyhi_out[n],betay[n]);
if (nyhi_in[n] < 0)
nyhi_in[n] = alpha[n] - 1;
}
if (!domain->zperiodic) {
if (nzlo_in[n] == 0)
nzlo_in[n] = alpha[n];
nzlo_out[n] = MAX(nzlo_out[n],alpha[n]);
if (nzhi_in[n] == nz_msm[n] - 1)
nzhi_in[n] = betaz[n];
nzhi_out[n] = MIN(nzhi_out[n],betaz[n]);
if (nzhi_in[n] < 0)
nzhi_in[n] = alpha[n] - 1;
}
}
// MSM grids for this proc, including ghosts
ngrid[n] = (nxhi_out[n]-nxlo_out[n]+1) * (nyhi_out[n]-nylo_out[n]+1) *
(nzhi_out[n]-nzlo_out[n]+1);
}
}
/* ----------------------------------------------------------------------
reset local grid arrays and communication stencils
called by fix balance b/c it changed sizes of processor sub-domains
------------------------------------------------------------------------- */
void MSM::setup_grid()
{
// free all arrays previously allocated
// pre-compute volume-dependent coeffs
// reset portion of global grid that each proc owns
// reallocate MSM long-range dependent memory
// don't invoke allocate_peratom(), compute() will allocate when needed
setup();
}
/* ----------------------------------------------------------------------
check if all factors of n are in list of factors
return 1 if yes, 0 if no
------------------------------------------------------------------------- */
int MSM::factorable(int n, int &flag, int &levels)
{
int i,norig;
norig = n;
levels = 1;
while (n > 1) {
for (i = 0; i < nfactors; i++) {
if (n % factors[i] == 0) {
n /= factors[i];
levels++;
break;
}
}
if (i == nfactors) {
flag = 1;
return 0;
}
}
return 1;
}
/* ----------------------------------------------------------------------
find center grid pt for each of my particles
check that full stencil for the particle will fit in my 3d brick
store central grid pt indices in part2grid array
------------------------------------------------------------------------- */
void MSM::particle_map()
{
int nx,ny,nz;
double **x = atom->x;
int nlocal = atom->nlocal;
int flag = 0;
for (int i = 0; i < nlocal; i++) {
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// current particle coord can be outside global and local box
// add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
nx = static_cast<int> ((x[i][0]-boxlo[0])*delxinv[0]+OFFSET) - OFFSET;
ny = static_cast<int> ((x[i][1]-boxlo[1])*delyinv[0]+OFFSET) - OFFSET;
nz = static_cast<int> ((x[i][2]-boxlo[2])*delzinv[0]+OFFSET) - OFFSET;
part2grid[i][0] = nx;
part2grid[i][1] = ny;
part2grid[i][2] = nz;
// check that entire stencil around nx,ny,nz will fit in my 3d brick
if (nx+nlower < nxlo_out[0] || nx+nupper > nxhi_out[0] ||
ny+nlower < nylo_out[0] || ny+nupper > nyhi_out[0] ||
nz+nlower < nzlo_out[0] || nz+nupper > nzhi_out[0]) flag = 1;
}
if (flag) error->one(FLERR,"Out of range atoms - cannot compute MSM");
}
/* ----------------------------------------------------------------------
create discretized "density" on section of global grid due to my particles
density(x,y,z) = charge "density" at grid points of my 3d brick
(nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
in global grid
------------------------------------------------------------------------- */
void MSM::make_rho()
{
//fprintf(screen,"MSM aninterpolation\n\n");
int i,l,m,n,nn,nx,ny,nz,mx,my,mz;
double dx,dy,dz,x0,y0,z0;
// clear 3d density array
double ***qgridn = qgrid[0];
memset(&(qgridn[nzlo_out[0]][nylo_out[0]][nxlo_out[0]]),0,ngrid[0]*sizeof(double));
// loop over my charges, add their contribution to nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
double *q = atom->q;
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx - (x[i][0]-boxlo[0])*delxinv[0];
dy = ny - (x[i][1]-boxlo[1])*delyinv[0];
dz = nz - (x[i][2]-boxlo[2])*delzinv[0];
compute_phis_and_dphis(dx,dy,dz);
z0 = q[i];
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
y0 = z0*phi1d[2][n];
for (m = nlower; m <= nupper; m++) {
my = m+ny;
x0 = y0*phi1d[1][m];
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
qgridn[mz][my][mx] += x0*phi1d[0][l];
}
}
}
}
}
/* ----------------------------------------------------------------------
MSM direct part procedure for intermediate grid levels
------------------------------------------------------------------------- */
void MSM::direct(int n)
{
//fprintf(screen,"Direct contribution on level %i\n\n",n);
double ***egridn = egrid[n];
double ***qgridn = qgrid[n];
// zero out electric potential
memset(&(egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
int jj,kk;
int imin,imax,jmin,jmax,kmin,kmax;
double qtmp;
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
int nx = nxhi_direct - nxlo_direct + 1;
int ny = nyhi_direct - nylo_direct + 1;
int nz = nzhi_direct - nzlo_direct + 1;
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
if (domain->zperiodic) {
kmin = nzlo_direct;
kmax = nzhi_direct;
} else {
kmin = MAX(nzlo_direct,alpha[n] - icz);
kmax = MIN(nzhi_direct,betaz[n] - icz);
}
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) {
if (domain->yperiodic) {
jmin = nylo_direct;
jmax = nyhi_direct;
} else {
jmin = MAX(nylo_direct,alpha[n] - icy);
jmax = MIN(nyhi_direct,betay[n] - icy);
}
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) {
if (domain->xperiodic) {
imin = nxlo_direct;
imax = nxhi_direct;
} else {
imin = MAX(nxlo_direct,alpha[n] - icx);
imax = MIN(nxhi_direct,betax[n] - icx);
}
if (evflag) {
esum = 0.0;
v0sum = v1sum = v2sum = 0.0;
v3sum = v4sum = v5sum = 0.0;
}
for (iz = kmin; iz <= kmax; iz++) {
kk = icz+iz;
zk = (iz + nzhi_direct)*ny;
for (iy = jmin; iy <= jmax; iy++) {
jj = icy+iy;
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = imin; ix <= imax; ix++) {
qtmp = qgridn[kk][jj][icx+ix];
k = zyk + ix + nxhi_direct;
egridn[icz][icy][icx] += g_direct[n][k] * qtmp;
if (evflag) {
if (eflag_global) esum += g_direct[n][k] * qtmp;
if (vflag_global) {
v0sum += v0_direct[n][k] * qtmp;
v1sum += v1_direct[n][k] * qtmp;
v2sum += v2_direct[n][k] * qtmp;
v3sum += v3_direct[n][k] * qtmp;
v4sum += v4_direct[n][k] * qtmp;
v5sum += v5_direct[n][k] * qtmp;
}
}
}
}
}
if (evflag) {
qtmp = qgridn[icz][icy][icx];
if (eflag_global) energy += esum * qtmp;
if (vflag_global) {
virial[0] += v0sum * qtmp;
virial[1] += v1sum * qtmp;
virial[2] += v2sum * qtmp;
virial[3] += v3sum * qtmp;
virial[4] += v4sum * qtmp;
virial[5] += v5sum * qtmp;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
MSM direct part procedure for intermediate grid levels
------------------------------------------------------------------------- */
void MSM::direct_peratom(int n)
{
double ***qgridn = qgrid[n];
double ***v0gridn = v0grid[n];
double ***v1gridn = v1grid[n];
double ***v2gridn = v2grid[n];
double ***v3gridn = v3grid[n];
double ***v4gridn = v4grid[n];
double ***v5gridn = v5grid[n];
// zero out virial
if (vflag_atom) {
memset(&(v0gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v1gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v2gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v3gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v4gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v5gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
}
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
int jj,kk;
int imin,imax,jmin,jmax,kmin,kmax;
double qtmp;
int nx = nxhi_direct - nxlo_direct + 1;
int ny = nyhi_direct - nylo_direct + 1;
int nz = nzhi_direct - nzlo_direct + 1;
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
if (domain->zperiodic) {
kmin = nzlo_direct;
kmax = nzhi_direct;
} else {
kmin = MAX(nzlo_direct,alpha[n] - icz);
kmax = MIN(nzhi_direct,betaz[n] - icz);
}
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) {
if (domain->yperiodic) {
jmin = nylo_direct;
jmax = nyhi_direct;
} else {
jmin = MAX(nylo_direct,alpha[n] - icy);
jmax = MIN(nyhi_direct,betay[n] - icy);
}
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) {
if (domain->xperiodic) {
imin = nxlo_direct;
imax = nxhi_direct;
} else {
imin = MAX(nxlo_direct,alpha[n] - icx);
imax = MIN(nxhi_direct,betax[n] - icx);
}
for (iz = kmin; iz <= kmax; iz++) {
kk = icz+iz;
zk = (iz + nzhi_direct)*ny;
for (iy = jmin; iy <= jmax; iy++) {
jj = icy+iy;
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = imin; ix <= imax; ix++) {
qtmp = qgridn[kk][jj][icx+ix];
k = zyk + ix + nxhi_direct;
v0gridn[icz][icy][icx] += v0_direct[n][k] * qtmp;
v1gridn[icz][icy][icx] += v1_direct[n][k] * qtmp;
v2gridn[icz][icy][icx] += v2_direct[n][k] * qtmp;
v3gridn[icz][icy][icx] += v3_direct[n][k] * qtmp;
v4gridn[icz][icy][icx] += v4_direct[n][k] * qtmp;
v5gridn[icz][icy][icx] += v5_direct[n][k] * qtmp;
}
}
}
}
}
}
}
/* ----------------------------------------------------------------------
MSM direct part procedure for top grid level
------------------------------------------------------------------------- */
void MSM::direct_top(int n)
{
//fprintf(screen,"Direct contribution on level %i\n\n",n);
double ***egridn = egrid[n];
double ***qgridn = qgrid[n];
// zero out electric potential
memset(&(egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
int jj,kk;
int imin,imax,jmin,jmax,kmin,kmax;
double qtmp;
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
int nx_top = betax[n] - alpha[n];
int ny_top = betay[n] - alpha[n];
int nz_top = betaz[n] - alpha[n];
int nx = 2*nx_top + 1;
int ny = 2*ny_top + 1;
int nz = 2*nz_top + 1;
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
kmin = alpha[n] - icz;
kmax = betaz[n] - icz;
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) {
jmin = alpha[n] - icy;
jmax = betay[n] - icy;
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) {
imin = alpha[n] - icx;
imax = betax[n] - icx;
if (evflag) {
esum = 0.0;
v0sum = v1sum = v2sum = 0.0;
v3sum = v4sum = v5sum = 0.0;
}
for (iz = kmin; iz <= kmax; iz++) {
kk = icz+iz;
zk = (iz + nz_top)*ny;
for (iy = jmin; iy <= jmax; iy++) {
jj = icy+iy;
zyk = (zk + iy + ny_top)*nx;
for (ix = imin; ix <= imax; ix++) {
qtmp = qgridn[kk][jj][icx+ix];
k = zyk + ix + nx_top;
qtmp = qgridn[kk][jj][icx+ix];
egridn[icz][icy][icx] += g_direct_top[k] * qtmp;
if (evflag) {
if (eflag_global) esum += g_direct_top[k] * qtmp;
if (vflag_global) {
v0sum += v0_direct_top[k] * qtmp;
v1sum += v1_direct_top[k] * qtmp;
v2sum += v2_direct_top[k] * qtmp;
v3sum += v3_direct_top[k] * qtmp;
v4sum += v4_direct_top[k] * qtmp;
v5sum += v5_direct_top[k] * qtmp;
}
}
}
}
}
if (evflag) {
qtmp = qgridn[icz][icy][icx];
if (eflag_global) energy += esum * qtmp;
if (vflag_global) {
virial[0] += v0sum * qtmp;
virial[1] += v1sum * qtmp;
virial[2] += v2sum * qtmp;
virial[3] += v3sum * qtmp;
virial[4] += v4sum * qtmp;
virial[5] += v5sum * qtmp;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
MSM direct part procedure for top grid level
------------------------------------------------------------------------- */
void MSM::direct_peratom_top(int n)
{
//fprintf(screen,"Direct contribution on level %i\n\n",n);
double ***qgridn = qgrid[n];
double ***v0gridn = v0grid[n];
double ***v1gridn = v1grid[n];
double ***v2gridn = v2grid[n];
double ***v3gridn = v3grid[n];
double ***v4gridn = v4grid[n];
double ***v5gridn = v5grid[n];
// zero out virial
if (vflag_atom) {
memset(&(v0gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v1gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v2gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v3gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v4gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
memset(&(v5gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
}
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
int jj,kk;
int imin,imax,jmin,jmax,kmin,kmax;
double qtmp;
int nx_top = betax[n] - alpha[n];
int ny_top = betay[n] - alpha[n];
int nz_top = betaz[n] - alpha[n];
int nx = 2*nx_top + 1;
int ny = 2*ny_top + 1;
int nz = 2*nz_top + 1;
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
kmin = alpha[n] - icz;
kmax = betaz[n] - icz;
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) {
jmin = alpha[n] - icy;
jmax = betay[n] - icy;
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) {
imin = alpha[n] - icx;
imax = betax[n] - icx;
for (iz = kmin; iz <= kmax; iz++) {
kk = icz+iz;
zk = (iz + nz_top)*ny;
for (iy = jmin; iy <= jmax; iy++) {
jj = icy+iy;
zyk = (zk + iy + ny_top)*nx;
for (ix = imin; ix <= imax; ix++) {
qtmp = qgridn[kk][jj][icx+ix];
k = zyk + ix + nx_top;
qtmp = qgridn[kk][jj][icx+ix];
v0gridn[icz][icy][icx] += v0_direct_top[k] * qtmp;
v1gridn[icz][icy][icx] += v1_direct_top[k] * qtmp;
v2gridn[icz][icy][icx] += v2_direct_top[k] * qtmp;
v3gridn[icz][icy][icx] += v3_direct_top[k] * qtmp;
v4gridn[icz][icy][icx] += v4_direct_top[k] * qtmp;
v5gridn[icz][icy][icx] += v5_direct_top[k] * qtmp;
}
}
}
}
}
}
}
/* ----------------------------------------------------------------------
MSM restriction procedure for intermediate grid levels
------------------------------------------------------------------------- */
void MSM::restriction(int n)
{
//fprintf(screen,"Restricting from level %i to %i\n\n",n,n+1);
int p = order-1;
double ***qgrid1 = qgrid[n];
double ***qgrid2 = qgrid[n+1];
//restrict grid (going from grid n to grid n+1, i.e. to a coarser grid)
int k = 0;
int index[p+1];
for (int nu=-p; nu<=p; nu++) {
if (nu%2 == 0 && nu != 0) continue;
phi1d[0][k] = compute_phi(nu*delxinv[n+1]/delxinv[n]);
phi1d[1][k] = compute_phi(nu*delyinv[n+1]/delyinv[n]);
phi1d[2][k] = compute_phi(nu*delzinv[n+1]/delzinv[n]);
index[k] = nu;
k++;
}
int ip,jp,kp,ic,jc,kc,i,j;
int ii,jj,kk;
double phiz,phizy;
// zero out charge on coarser grid
memset(&(qgrid2[nzlo_out[n+1]][nylo_out[n+1]][nxlo_out[n+1]]),0,ngrid[n+1]*sizeof(double));
for (kp = nzlo_in[n+1]; kp <= nzhi_in[n+1]; kp++)
for (jp = nylo_in[n+1]; jp <= nyhi_in[n+1]; jp++)
for (ip = nxlo_in[n+1]; ip <= nxhi_in[n+1]; ip++) {
ic = ip * static_cast<int> (delxinv[n]/delxinv[n+1]);
jc = jp * static_cast<int> (delyinv[n]/delyinv[n+1]);
kc = kp * static_cast<int> (delzinv[n]/delzinv[n+1]);
for (k=0; k<=p+1; k++) {
kk = kc+index[k];
if (!domain->zperiodic) {
if (kk < alpha[n]) continue;
if (kk > betaz[n]) break;
}
phiz = phi1d[2][k];
for (j=0; j<=p+1; j++) {
jj = jc+index[j];
if (!domain->yperiodic) {
if (jj < alpha[n]) continue;
if (jj > betay[n]) break;
}
phizy = phi1d[1][j]*phiz;
for (i=0; i<=p+1; i++) {
ii = ic+index[i];
if (!domain->xperiodic) {
if (ii < alpha[n]) continue;
if (ii > betax[n]) break;
}
qgrid2[kp][jp][ip] += qgrid1[kk][jj][ii] *
phi1d[0][i]*phizy;
}
}
}
}
}
/* ----------------------------------------------------------------------
MSM prolongation procedure for intermediate grid levels
------------------------------------------------------------------------- */
void MSM::prolongation(int n)
{
//fprintf(screen,"Prolongating from level %i to %i\n\n",n+1,n);
int p = order-1;
double ***egrid1 = egrid[n];
double ***egrid2 = egrid[n+1];
double ***v0grid1 = v0grid[n];
double ***v0grid2 = v0grid[n+1];
double ***v1grid1 = v1grid[n];
double ***v1grid2 = v1grid[n+1];
double ***v2grid1 = v2grid[n];
double ***v2grid2 = v2grid[n+1];
double ***v3grid1 = v3grid[n];
double ***v3grid2 = v3grid[n+1];
double ***v4grid1 = v4grid[n];
double ***v4grid2 = v4grid[n+1];
double ***v5grid1 = v5grid[n];
double ***v5grid2 = v5grid[n+1];
//prolongate grid (going from grid n to grid n-1, i.e. to a finer grid)
int k = 0;
int index[p+1];
for (int nu=-p; nu<=p; nu++) {
if (nu%2 == 0 && nu != 0) continue;
phi1d[0][k] = compute_phi(nu*delxinv[n+1]/delxinv[n]);
phi1d[1][k] = compute_phi(nu*delyinv[n+1]/delyinv[n]);
phi1d[2][k] = compute_phi(nu*delzinv[n+1]/delzinv[n]);
index[k] = nu;
k++;
}
int ip,jp,kp,ic,jc,kc,i,j;
int ii,jj,kk;
double phiz,phizy,phi3d;
for (kp = nzlo_in[n+1]; kp <= nzhi_in[n+1]; kp++)
for (jp = nylo_in[n+1]; jp <= nyhi_in[n+1]; jp++)
for (ip = nxlo_in[n+1]; ip <= nxhi_in[n+1]; ip++) {
ic = ip * static_cast<int> (delxinv[n]/delxinv[n+1]);
jc = jp * static_cast<int> (delyinv[n]/delyinv[n+1]);
kc = kp * static_cast<int> (delzinv[n]/delzinv[n+1]);
for (k=0; k<=p+1; k++) {
kk = kc+index[k];
if (!domain->zperiodic) {
if (kk < alpha[n]) continue;
if (kk > betaz[n]) break;
}
phiz = phi1d[2][k];
for (j=0; j<=p+1; j++) {
jj = jc+index[j];
if (!domain->yperiodic) {
if (jj < alpha[n]) continue;
if (jj > betay[n]) break;
}
phizy = phi1d[1][j]*phiz;
for (i=0; i<=p+1; i++) {
ii = ic+index[i];
if (!domain->xperiodic) {
if (ii < alpha[n]) continue;
if (ii > betax[n]) break;
}
phi3d = phi1d[0][i]*phizy;
egrid1[kk][jj][ii] += egrid2[kp][jp][ip] * phi3d;
if (vflag_atom) {
v0grid1[kk][jj][ii] += v0grid2[kp][jp][ip] * phi3d;
v1grid1[kk][jj][ii] += v1grid2[kp][jp][ip] * phi3d;
v2grid1[kk][jj][ii] += v2grid2[kp][jp][ip] * phi3d;
v3grid1[kk][jj][ii] += v3grid2[kp][jp][ip] * phi3d;
v4grid1[kk][jj][ii] += v4grid2[kp][jp][ip] * phi3d;
v5grid1[kk][jj][ii] += v5grid2[kp][jp][ip] * phi3d;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
pack own values to buf to send to another proc
------------------------------------------------------------------------- */
void MSM::pack_forward(int flag, double *buf, int nlist, int *list)
{
int n = current_level;
double ***qgridn = qgrid[n];
double ***egridn = egrid[n];
double ***v0gridn = v0grid[n];
double ***v1gridn = v1grid[n];
double ***v2gridn = v2grid[n];
double ***v3gridn = v3grid[n];
double ***v4gridn = v4grid[n];
double ***v5gridn = v5grid[n];
int k = 0;
if (flag == FORWARD_RHO) {
double *qsrc = &qgridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
buf[k++] = qsrc[list[i]];
}
} else if (flag == FORWARD_AD) {
double *src = &egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++)
buf[i] = src[list[i]];
} else if (flag == FORWARD_AD_PERATOM) {
double *v0src = &v0gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v1src = &v1gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v2src = &v2gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v3src = &v3gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v4src = &v4gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v5src = &v5gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
buf[k++] = v0src[list[i]];
buf[k++] = v1src[list[i]];
buf[k++] = v2src[list[i]];
buf[k++] = v3src[list[i]];
buf[k++] = v4src[list[i]];
buf[k++] = v5src[list[i]];
}
}
}
/* ----------------------------------------------------------------------
unpack another proc's own values from buf and set own ghost values
------------------------------------------------------------------------- */
void MSM::unpack_forward(int flag, double *buf, int nlist, int *list)
{
int n = current_level;
double ***qgridn = qgrid[n];
double ***egridn = egrid[n];
double ***v0gridn = v0grid[n];
double ***v1gridn = v1grid[n];
double ***v2gridn = v2grid[n];
double ***v3gridn = v3grid[n];
double ***v4gridn = v4grid[n];
double ***v5gridn = v5grid[n];
int k = 0;
if (flag == FORWARD_RHO) {
double *dest = &qgridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
dest[list[i]] = buf[k++];
}
} else if (flag == FORWARD_AD) {
double *dest = &egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++)
dest[list[i]] = buf[k++];
} else if (flag == FORWARD_AD_PERATOM) {
double *v0src = &v0gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v1src = &v1gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v2src = &v2gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v3src = &v3gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v4src = &v4gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v5src = &v5gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
v0src[list[i]] = buf[k++];
v1src[list[i]] = buf[k++];
v2src[list[i]] = buf[k++];
v3src[list[i]] = buf[k++];
v4src[list[i]] = buf[k++];
v5src[list[i]] = buf[k++];
}
}
}
/* ----------------------------------------------------------------------
pack ghost values into buf to send to another proc
------------------------------------------------------------------------- */
void MSM::pack_reverse(int flag, double *buf, int nlist, int *list)
{
int n = current_level;
double ***qgridn = qgrid[n];
double ***egridn = egrid[n];
double ***v0gridn = v0grid[n];
double ***v1gridn = v1grid[n];
double ***v2gridn = v2grid[n];
double ***v3gridn = v3grid[n];
double ***v4gridn = v4grid[n];
double ***v5gridn = v5grid[n];
int k = 0;
if (flag == REVERSE_RHO) {
double *qsrc = &qgridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
buf[k++] = qsrc[list[i]];
}
} else if (flag == REVERSE_AD) {
double *src = &egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++)
buf[i] = src[list[i]];
} else if (flag == REVERSE_AD_PERATOM) {
double *v0src = &v0gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v1src = &v1gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v2src = &v2gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v3src = &v3gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v4src = &v4gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v5src = &v5gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
buf[k++] = v0src[list[i]];
buf[k++] = v1src[list[i]];
buf[k++] = v2src[list[i]];
buf[k++] = v3src[list[i]];
buf[k++] = v4src[list[i]];
buf[k++] = v5src[list[i]];
}
}
}
/* ----------------------------------------------------------------------
unpack another proc's ghost values from buf and add to own values
------------------------------------------------------------------------- */
void MSM::unpack_reverse(int flag, double *buf, int nlist, int *list)
{
int n = current_level;
double ***qgridn = qgrid[n];
double ***egridn = egrid[n];
double ***v0gridn = v0grid[n];
double ***v1gridn = v1grid[n];
double ***v2gridn = v2grid[n];
double ***v3gridn = v3grid[n];
double ***v4gridn = v4grid[n];
double ***v5gridn = v5grid[n];
int k = 0;
if (flag == REVERSE_RHO) {
double *dest = &qgridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
dest[list[i]] += buf[k++];
}
} else if (flag == REVERSE_AD) {
double *dest = &egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++)
dest[list[i]] += buf[k++];
} else if (flag == REVERSE_AD_PERATOM) {
double *v0src = &v0gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v1src = &v1gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v2src = &v2gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v3src = &v3gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v4src = &v4gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
double *v5src = &v5gridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]];
for (int i = 0; i < nlist; i++) {
v0src[list[i]] += buf[k++];
v1src[list[i]] += buf[k++];
v2src[list[i]] += buf[k++];
v3src[list[i]] += buf[k++];
v4src[list[i]] += buf[k++];
v5src[list[i]] += buf[k++];
}
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get force on my particles
------------------------------------------------------------------------- */
void MSM::fieldforce()
{
//fprintf(screen,"MSM interpolation\n\n");
double ***egridn = egrid[0];
int i,l,m,n,nx,ny,nz,mx,my,mz;
double dx,dy,dz;
double phi_x,phi_y,phi_z;
double dphi_x,dphi_y,dphi_z;
double ekx,eky,ekz;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx - (x[i][0]-boxlo[0])*delxinv[0];
dy = ny - (x[i][1]-boxlo[1])*delyinv[0];
dz = nz - (x[i][2]-boxlo[2])*delzinv[0];
compute_phis_and_dphis(dx,dy,dz);
ekx = eky = ekz = 0.0;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
phi_z = phi1d[2][n];
dphi_z = dphi1d[2][n];
for (m = nlower; m <= nupper; m++) {
my = m+ny;
phi_y = phi1d[1][m];
dphi_y = dphi1d[1][m];
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
phi_x = phi1d[0][l];
dphi_x = dphi1d[0][l];
ekx += dphi_x*phi_y*phi_z*egridn[mz][my][mx];
eky += phi_x*dphi_y*phi_z*egridn[mz][my][mx];
ekz += phi_x*phi_y*dphi_z*egridn[mz][my][mx];
}
}
}
ekx *= delxinv[0];
eky *= delyinv[0];
ekz *= delzinv[0];
// convert E-field to force
const double qfactor = force->qqrd2e*scale*q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
f[i][2] += qfactor*ekz;
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get per-atom energy/virial
------------------------------------------------------------------------- */
void MSM::fieldforce_peratom()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
double dx,dy,dz,x0,y0,z0;
double u,v0,v1,v2,v3,v4,v5;
double ***egridn = egrid[0];
double ***v0gridn = v0grid[0];
double ***v1gridn = v1grid[0];
double ***v2gridn = v2grid[0];
double ***v3gridn = v3grid[0];
double ***v4gridn = v4grid[0];
double ***v5gridn = v5grid[0];
// loop over my charges, interpolate from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
double *q = atom->q;
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx - (x[i][0]-boxlo[0])*delxinv[0];
dy = ny - (x[i][1]-boxlo[1])*delyinv[0];
dz = nz - (x[i][2]-boxlo[2])*delzinv[0];
compute_phis_and_dphis(dx,dy,dz);
u = v0 = v1 = v2 = v3 = v4 = v5 = 0.0;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
z0 = phi1d[2][n];
for (m = nlower; m <= nupper; m++) {
my = m+ny;
y0 = z0*phi1d[1][m];
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
x0 = y0*phi1d[0][l];
if (eflag_atom) u += x0*egridn[mz][my][mx];
if (vflag_atom) {
v0 += x0*v0gridn[mz][my][mx];
v1 += x0*v1gridn[mz][my][mx];
v2 += x0*v2gridn[mz][my][mx];
v3 += x0*v3gridn[mz][my][mx];
v4 += x0*v4gridn[mz][my][mx];
v5 += x0*v5gridn[mz][my][mx];
}
}
}
}
if (eflag_atom) eatom[i] += q[i]*u;
if (vflag_atom) {
vatom[i][0] += q[i]*v0;
vatom[i][1] += q[i]*v1;
vatom[i][2] += q[i]*v2;
vatom[i][3] += q[i]*v3;
vatom[i][4] += q[i]*v4;
vatom[i][5] += q[i]*v5;
}
}
}
/* ----------------------------------------------------------------------
charge assignment into phi1d
------------------------------------------------------------------------- */
void MSM::compute_phis_and_dphis(const double &dx, const double &dy,
const double &dz)
{
double delx,dely,delz;
for (int nu = nlower; nu <= nupper; nu++) {
delx = dx + double(nu);
dely = dy + double(nu);
delz = dz + double(nu);
phi1d[0][nu] = compute_phi(delx);
phi1d[1][nu] = compute_phi(dely);
phi1d[2][nu] = compute_phi(delz);
dphi1d[0][nu] = compute_dphi(delx);
dphi1d[1][nu] = compute_dphi(dely);
dphi1d[2][nu] = compute_dphi(delz);
}
}
/* ----------------------------------------------------------------------
compute phi using interpolating polynomial
see Eq 7 from Parallel Computing 35 (2009) 164177
and Hardy's thesis
------------------------------------------------------------------------- */
double MSM::compute_phi(const double &xi)
{
double phi;
double abs_xi = fabs(xi);
double xi2 = xi*xi;
if (order == 4) {
if (abs_xi <= 1) {
phi = (1.0 - abs_xi)*(1.0 + abs_xi - 1.5*xi2);
} else if (abs_xi <= 2) {
phi = -0.5*(abs_xi - 1.0)*(2.0 - abs_xi)*(2.0 - abs_xi);
} else {
phi = 0.0;
}
} else if (order == 6) {
if (abs_xi <= 1) {
phi = (1.0 - xi2)*(2.0 - abs_xi)*(6.0 + 3.0*abs_xi -
5.0*xi2)/12.0;
} else if (abs_xi <= 2) {
phi = -(abs_xi - 1.0)*(2.0 - abs_xi)*(3.0 - abs_xi)*
(4.0 + 9.0*abs_xi - 5.0*xi2)/24.0;
} else if (abs_xi <= 3) {
phi = (abs_xi - 1.0)*(abs_xi - 2.0)*(3.0 - abs_xi)*
(3.0 - abs_xi)*(4.0 - abs_xi)/24.0;
} else {
phi = 0.0;
}
} else if (order == 8) {
if (abs_xi <= 1) {
phi = (1.0 - xi2)*(4.0 - xi2)*(3.0 - abs_xi)*
(12.0 + 4.0*abs_xi - 7.0*xi2)/144.0;
} else if (abs_xi <= 2) {
phi = -(xi2 - 1.0)*(2.0 - abs_xi)*(3.0 - abs_xi)*
(4.0 - abs_xi)*(10.0 + 12.0*abs_xi - 7.0*xi2)/240.0;
} else if (abs_xi <= 3) {
phi = (abs_xi - 1.0)*(abs_xi - 2.0)*(3.0 - abs_xi)*(4.0 - abs_xi)*
(5.0 - abs_xi)*(6.0 + 20.0*abs_xi - 7.0*xi2)/720.0;
} else if (abs_xi <= 4) {
phi = -(abs_xi - 1.0)*(abs_xi - 2.0)*(abs_xi - 3.0)*(4.0 - abs_xi)*
(4.0 - abs_xi)*(5.0 - abs_xi)*(6.0 - abs_xi)/720.0;
} else {
phi = 0.0;
}
} else if (order == 10) {
if (abs_xi <= 1) {
phi = (1.0 - xi2)*(4.0 - xi2)*(9.0 - xi2)*
(4.0 - abs_xi)*(20.0 + 5.0*abs_xi - 9.0*xi2)/2880.0;
} else if (abs_xi <= 2) {
phi = -(xi2 - 1.0)*(4.0 - xi2)*(3.0 - abs_xi)*(4.0 - abs_xi)*
(5.0 - abs_xi)*(6.0 + 5.0*abs_xi - 3.0*xi2)/1440.0;
} else if (abs_xi <= 3) {
phi = (xi2 - 1.0)*(abs_xi - 2.0)*(3.0 - abs_xi)*(4.0 - abs_xi)*
(5.0 - abs_xi)*(6.0 - abs_xi)*(14.0 + 25.0*abs_xi - 9.0*xi2)/10080.0;
} else if (abs_xi <= 4) {
phi = -(abs_xi - 1.0)*(abs_xi - 2.0)*(abs_xi - 3.0)*(4.0 - abs_xi)*
(5.0 - abs_xi)*(6.0 - abs_xi)*(7.0 - abs_xi)*
(8.0 + 35.0*abs_xi - 9.0*xi2)/40320.0;
} else if (abs_xi <= 5) {
phi = (abs_xi - 1.0)*(abs_xi - 2.0)*(abs_xi - 3.0)*
(abs_xi - 4.0)*(5.0 - abs_xi)*(5.0 - abs_xi)*(6.0 - abs_xi)*
(7.0 - abs_xi)*(8.0 - abs_xi)/40320.0;
} else {
phi = 0.0;
}
}
return phi;
}
/* ----------------------------------------------------------------------
compute the derivative of phi
phi is an interpolating polynomial
see Eq 7 from Parallel Computing 35 (2009) 164177
and Hardy's thesis
------------------------------------------------------------------------- */
double MSM::compute_dphi(const double &xi)
{
double dphi;
double abs_xi = fabs(xi);
if (order == 4) {
double xi2 = xi*xi;
double abs_xi2 = abs_xi*abs_xi;
if (abs_xi == 0.0) {
dphi = 0.0;
} else if (abs_xi <= 1) {
dphi = xi*(3*xi2 + 6*abs_xi2 - 10*abs_xi)/2.0/abs_xi;
} else if (abs_xi <= 2) {
dphi = xi*(2 - abs_xi)*(3*abs_xi - 4)/2.0/abs_xi;
} else {
dphi = 0.0;
}
} else if (order == 6) {
double xi2 = xi*xi;
double xi4 = xi2*xi2;
double abs_xi2 = abs_xi*abs_xi;
double abs_xi3 = abs_xi2*abs_xi;
double abs_xi4 = abs_xi2*abs_xi2;
if (abs_xi == 0.0) {
dphi = 0.0;
} else if (abs_xi <= 1) {
dphi = xi*(46*xi2*abs_xi - 20*xi2*abs_xi2 - 5*xi4 + 5*xi2 +
6*abs_xi3 + 10*abs_xi2 - 50*abs_xi)/12.0/abs_xi;
} else if (abs_xi <= 2) {
dphi = xi*(15*xi2*abs_xi2 - 60*xi2*abs_xi + 55*xi2 +
10*abs_xi4 - 96*abs_xi3 + 260*abs_xi2 - 210*abs_xi + 10)/
24.0/abs_xi;
} else if (abs_xi <= 3) {
dphi = -xi*(abs_xi - 3)*(5*abs_xi3 - 37*abs_xi2 +
84*abs_xi - 58)/24.0/abs_xi;
} else {
dphi = 0.0;
}
} else if (order == 8) {
double xi2 = xi*xi;
double xi4 = xi2*xi2;
double xi6 = xi4*xi2;
double abs_xi2 = abs_xi*abs_xi;
double abs_xi3 = abs_xi2*abs_xi;
double abs_xi4 = abs_xi2*abs_xi2;
double abs_xi5 = abs_xi4*abs_xi;
double abs_xi6 = abs_xi5*abs_xi;
if (abs_xi == 0.0) {
dphi = 0.0;
} else if (abs_xi <= 1) {
dphi = xi*(7*xi6 + 42*xi4*abs_xi2 - 134*xi4*abs_xi - 35*xi4 -
16*xi2*abs_xi3 - 140*xi2*abs_xi2 + 604*xi2*abs_xi + 28*xi2 +
40*abs_xi3 + 56*abs_xi2 - 560*abs_xi)/144.0/abs_xi;
} else if (abs_xi <= 2) {
dphi = xi*(126*xi4*abs_xi - 21*xi4*abs_xi2 - 182*xi4 -
28*xi2*abs_xi4 + 300*xi2*abs_xi3 - 1001*xi2*abs_xi2 +
990*xi2*abs_xi + 154*xi2 + 24*abs_xi5 - 182*abs_xi4 +
270*abs_xi3 + 602*abs_xi2 - 1260*abs_xi + 28)/240.0/abs_xi;
} else if (abs_xi <= 3) {
dphi = xi*(35*xi2*abs_xi4 - 420*xi2*abs_xi3 +
1785*xi2*abs_xi2 - 3150*xi2*abs_xi + 1918*xi2 +
14*abs_xi6 - 330*abs_xi5 + 2660*abs_xi4 -
9590*abs_xi3 + 15806*abs_xi2 - 9940*abs_xi + 756)/720.0/abs_xi;
} else if (abs_xi <= 4) {
dphi = -xi*(abs_xi - 4)*(7*abs_xi5 - 122*abs_xi4 +
807*abs_xi3 - 2512*abs_xi2 + 3644*abs_xi - 1944)/720.0/abs_xi;
} else {
dphi = 0.0;
}
} else if (order == 10) {
double xi2 = xi*xi;
double xi4 = xi2*xi2;
double xi6 = xi4*xi2;
double xi8 = xi6*xi2;
double abs_xi2 = abs_xi*abs_xi;
double abs_xi3 = abs_xi2*abs_xi;
double abs_xi4 = abs_xi2*abs_xi2;
double abs_xi5 = abs_xi4*abs_xi;
double abs_xi6 = abs_xi5*abs_xi;
double abs_xi7 = abs_xi6*abs_xi;
double abs_xi8 = abs_xi7*abs_xi;
if (abs_xi == 0.0) {
dphi = 0.0;
} else if (abs_xi <= 1) {
dphi = xi*(298*xi6*abs_xi - 72*xi6*abs_xi2 - 9*xi8 +
126*xi6 + 30*xi4*abs_xi3 + 756*xi4*abs_xi2 - 3644*xi4*abs_xi -
441*xi4 - 280*xi2*abs_xi3 - 1764*xi2*abs_xi2 + 12026*xi2*abs_xi +
324*xi2 + 490*abs_xi3 + 648*abs_xi2 - 10792*abs_xi)/2880.0/abs_xi;
} else if (abs_xi <= 2) {
dphi = xi*(9*xi6*abs_xi2 - 72*xi6*abs_xi + 141*xi6 +
18*xi4*abs_xi4 - 236*xi4*abs_xi3 + 963*xi4*abs_xi2 -
1046*xi4*abs_xi - 687*xi4 - 20*xi2*abs_xi5 + 156*xi2*abs_xi4 +
168*xi2*abs_xi3 - 3522*xi2*abs_xi2 + 6382*xi2*abs_xi + 474*xi2 +
50*abs_xi5 - 516*abs_xi4 + 1262*abs_xi3 + 1596*abs_xi2 -
6344*abs_xi + 72)/1440.0/abs_xi;
} else if (abs_xi <= 3) {
dphi = xi*(720*xi4*abs_xi3 - 45*xi4*abs_xi4 - 4185*xi4*abs_xi2 +
10440*xi4*abs_xi - 9396*xi4 - 36*xi2*abs_xi6 + 870*xi2*abs_xi5 -
7965*xi2*abs_xi4 + 34540*xi2*abs_xi3 - 70389*xi2*abs_xi2 +
51440*xi2*abs_xi + 6012*xi2 + 50*abs_xi7 - 954*abs_xi6 +
6680*abs_xi5 - 19440*abs_xi4 + 11140*abs_xi3 + 49014*abs_xi2 -
69080*abs_xi + 3384)/10080.0/abs_xi;
} else if (abs_xi <= 4) {
dphi = xi*(63*xi2*abs_xi6 - 1512*xi2*abs_xi5 + 14490*xi2*abs_xi4 -
70560*xi2*abs_xi3 + 182763*xi2*abs_xi2 - 236376*xi2*abs_xi +
117612*xi2 + 18*abs_xi8 - 784*abs_xi7 + 12600*abs_xi6 -
101556*abs_xi5 + 451962*abs_xi4 - 1121316*abs_xi3 +
1451628*abs_xi2 - 795368*abs_xi + 71856)/40320.0/abs_xi;
} else if (abs_xi <= 5) {
dphi = -xi*(abs_xi - 5)*(9*abs_xi7 - 283*abs_xi6 +
3667*abs_xi5 - 25261*abs_xi4 + 99340*abs_xi3 -
221416*abs_xi2 + 256552*abs_xi - 117648)/40320.0/abs_xi;
} else {
dphi = 0.0;
}
}
return dphi;
}
/* ----------------------------------------------------------------------
Compute direct interaction for intermediate grid levels
------------------------------------------------------------------------- */
void MSM::get_g_direct()
{
if (g_direct) memory->destroy(g_direct);
memory->create(g_direct,levels,nmax_direct,"msm:g_direct");
double a = cutoff;
int n,zk,zyk,k,ix,iy,iz;
double xdiff,ydiff,zdiff;
double rsq,rho,two_n;
two_n = 1.0;
int nx = nxhi_direct - nxlo_direct + 1;
int ny = nyhi_direct - nylo_direct + 1;
int nz = nzhi_direct - nzlo_direct + 1;
for (n=0; n<levels; n++) {
for (iz = nzlo_direct; iz <= nzhi_direct; iz++) {
zdiff = iz/delzinv[n];
zk = (iz + nzhi_direct)*ny;
for (iy = nylo_direct; iy <= nyhi_direct; iy++) {
ydiff = iy/delyinv[n];
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = nxlo_direct; ix <= nxhi_direct; ix++) {
xdiff = ix/delxinv[n];
rsq = xdiff*xdiff + ydiff*ydiff + zdiff*zdiff;
rho = sqrt(rsq)/(two_n*a);
k = zyk + ix + nxhi_direct;
g_direct[n][k] = gamma(rho)/(two_n*a) - gamma(rho/2.0)/(2.0*two_n*a);
}
}
}
two_n *= 2.0;
}
}
/* ----------------------------------------------------------------------
Compute direct interaction for intermediate grid levels
------------------------------------------------------------------------- */
void MSM::get_virial_direct()
{
if (v0_direct) memory->destroy(v0_direct);
memory->create(v0_direct,levels,nmax_direct,"msm:v0_direct");
if (v1_direct) memory->destroy(v1_direct);
memory->create(v1_direct,levels,nmax_direct,"msm:v1_direct");
if (v2_direct) memory->destroy(v2_direct);
memory->create(v2_direct,levels,nmax_direct,"msm:v2_direct");
if (v3_direct) memory->destroy(v3_direct);
memory->create(v3_direct,levels,nmax_direct,"msm:v3_direct");
if (v4_direct) memory->destroy(v4_direct);
memory->create(v4_direct,levels,nmax_direct,"msm:v4_direct");
if (v5_direct) memory->destroy(v5_direct);
memory->create(v5_direct,levels,nmax_direct,"msm:v5_direct");
double a = cutoff;
double a_sq = cutoff*cutoff;
int n,zk,zyk,k,ix,iy,iz;
double xdiff,ydiff,zdiff;
double rsq,r,rho,two_n,two_nsq,dg;
two_n = 1.0;
int nx = nxhi_direct - nxlo_direct + 1;
int ny = nyhi_direct - nylo_direct + 1;
int nz = nzhi_direct - nzlo_direct + 1;
for (n=0; n<levels; n++) {
two_nsq = two_n * two_n;
for (iz = nzlo_direct; iz <= nzhi_direct; iz++) {
zdiff = iz/delzinv[n];
zk = (iz + nzhi_direct)*ny;
for (iy = nylo_direct; iy <= nyhi_direct; iy++) {
ydiff = iy/delyinv[n];
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = nxlo_direct; ix <= nxhi_direct; ix++) {
xdiff = ix/delxinv[n];
rsq = xdiff*xdiff + ydiff*ydiff + zdiff*zdiff;
k = zyk + ix + nxhi_direct;
r = sqrt(rsq);
if (r == 0) {
v0_direct[n][k] = 0.0;
v1_direct[n][k] = 0.0;
v2_direct[n][k] = 0.0;
v3_direct[n][k] = 0.0;
v4_direct[n][k] = 0.0;
v5_direct[n][k] = 0.0;
} else {
rho = r/(two_n*a);
dg = -(dgamma(rho)/(two_nsq*a_sq) -
dgamma(rho/2.0)/(4.0*two_nsq*a_sq))/r;
v0_direct[n][k] = dg * xdiff * xdiff;
v1_direct[n][k] = dg * ydiff * ydiff;
v2_direct[n][k] = dg * zdiff * zdiff;
v3_direct[n][k] = dg * xdiff * ydiff;
v4_direct[n][k] = dg * xdiff * zdiff;
v5_direct[n][k] = dg * ydiff * zdiff;
}
}
}
}
two_n *= 2.0;
}
}
/* ----------------------------------------------------------------------
Compute direct interaction for top grid level
------------------------------------------------------------------------- */
void MSM::get_g_direct_top(int n)
{
int nx_top = betax[n] - alpha[n];
int ny_top = betay[n] - alpha[n];
int nz_top = betaz[n] - alpha[n];
int nx = 2*nx_top + 1;
int ny = 2*ny_top + 1;
int nz = 2*nz_top + 1;
int nmax_top = 8*(nx+1)*(ny*1)*(nz+1);
if (g_direct_top) memory->destroy(g_direct_top);
memory->create(g_direct_top,nmax_top,"msm:g_direct_top");
double a = cutoff;
int zk,zyk,k,ix,iy,iz;
double xdiff,ydiff,zdiff;
double rsq,rho,two_n;
two_n = pow(2.0,n);
for (iz = -nz_top; iz <= nz_top; iz++) {
zdiff = iz/delzinv[n];
zk = (iz + nz_top)*ny;
for (iy = -ny_top; iy <= ny_top; iy++) {
ydiff = iy/delyinv[n];
zyk = (zk + iy + ny_top)*nx;
for (ix = -nx_top; ix <= nx_top; ix++) {
xdiff = ix/delxinv[n];
rsq = xdiff*xdiff + ydiff*ydiff + zdiff*zdiff;
rho = sqrt(rsq)/(two_n*a);
k = zyk + ix + nx_top;
g_direct_top[k] = gamma(rho)/(two_n*a);
}
}
}
}
/* ----------------------------------------------------------------------
Compute direct interaction for top grid level
------------------------------------------------------------------------- */
void MSM::get_virial_direct_top(int n)
{
int nx_top = betax[n] - alpha[n];
int ny_top = betay[n] - alpha[n];
int nz_top = betaz[n] - alpha[n];
int nx = 2*nx_top + 1;
int ny = 2*ny_top + 1;
int nz = 2*nz_top + 1;
int nmax_top = 8*(nx+1)*(ny*1)*(nz+1);
if (v0_direct_top) memory->destroy(v0_direct_top);
memory->create(v0_direct_top,nmax_top,"msm:v0_direct_top");
if (v1_direct_top) memory->destroy(v1_direct_top);
memory->create(v1_direct_top,nmax_top,"msm:v1_direct_top");
if (v2_direct_top) memory->destroy(v2_direct_top);
memory->create(v2_direct_top,nmax_top,"msm:v2_direct_top");
if (v3_direct_top) memory->destroy(v3_direct_top);
memory->create(v3_direct_top,nmax_top,"msm:v3_direct_top");
if (v4_direct_top) memory->destroy(v4_direct_top);
memory->create(v4_direct_top,nmax_top,"msm:v4_direct_top");
if (v5_direct_top) memory->destroy(v5_direct_top);
memory->create(v5_direct_top,nmax_top,"msm:v5_direct_top");
double a = cutoff;
double a_sq = cutoff*cutoff;
int zk,zyk,k,ix,iy,iz;
double xdiff,ydiff,zdiff;
double rsq,r,rho,two_n,two_nsq,dg;
two_n = pow(2.0,n);
two_nsq = two_n * two_n;
for (iz = -nz_top; iz <= nz_top; iz++) {
zdiff = iz/delzinv[n];
zk = (iz + nz_top)*ny;
for (iy = -ny_top; iy <= ny_top; iy++) {
ydiff = iy/delyinv[n];
zyk = (zk + iy + ny_top)*nx;
for (ix = -nx_top; ix <= nx_top; ix++) {
xdiff = ix/delxinv[n];
rsq = xdiff*xdiff + ydiff*ydiff + zdiff*zdiff;
k = zyk + ix + nx_top;
r = sqrt(rsq);
if (r == 0) {
v0_direct_top[k] = 0.0;
v1_direct_top[k] = 0.0;
v2_direct_top[k] = 0.0;
v3_direct_top[k] = 0.0;
v4_direct_top[k] = 0.0;
v5_direct_top[k] = 0.0;
} else {
rho = r/(two_n*a);
dg = -(dgamma(rho)/(two_nsq*a_sq))/r;
v0_direct_top[k] = dg * xdiff * xdiff;
v1_direct_top[k] = dg * ydiff * ydiff;
v2_direct_top[k] = dg * zdiff * zdiff;
v3_direct_top[k] = dg * xdiff * ydiff;
v4_direct_top[k] = dg * xdiff * zdiff;
v5_direct_top[k] = dg * ydiff * zdiff;
}
}
}
}
}
Event Timeline
Log In to Comment