Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F67086229
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
Wed, Jun 19, 20:15
Size
99 KB
Mime Type
text/x-c
Expires
Fri, Jun 21, 20:15 (2 d)
Engine
blob
Format
Raw Data
Handle
18336655
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 "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
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);
procneigh_levels = NULL;
world_levels = NULL;
active_flag = NULL;
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_all = cg_peratom_all = 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
triclinic_check();
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,"Slab correction not needed for 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)");
pair_check();
// extract short-range Coulombic cutoff from pair style
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_coul",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];
}
qqrd2e = force->qqrd2e;
scale = 1.0;
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);
}
triclinic = domain->triclinic;
// 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(h_x,xprd);
double error_y = estimate_1d_error(h_y,yprd);
double error_z = estimate_1d_error(h_z,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()
{
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 and compute grid spacing
for (int n=0; n<levels; n++) {
if (triclinic == 0) {
delxinv[n] = nx_msm[n]/xprd;
delyinv[n] = ny_msm[n]/yprd;
delzinv[n] = nz_msm[n]/zprd;
} else { // use lamda (0-1) coordinates
delxinv[n] = nx_msm[n];
delyinv[n] = ny_msm[n];
delzinv[n] = nz_msm[n];
}
}
double ax = a;
double ay = a;
double az = a;
// transform the interaction sphere in box coords to an ellipsoid in lamda (0-1) coords to
// get the direct sum interaction limits for a triclinic system
if (triclinic) {
double tmp[3];
kspacebbox(a,&tmp[0]);
ax = tmp[0];
ay = tmp[1];
az = tmp[2];
}
// direct sum interaction limits
nxhi_direct = static_cast<int> (2.0*ax*delxinv[0]);
nxlo_direct = -nxhi_direct;
nyhi_direct = static_cast<int> (2.0*ay*delyinv[0]);
nylo_direct = -nyhi_direct;
nzhi_direct = static_cast<int> (2.0*az*delzinv[0]);
nzlo_direct = -nzhi_direct;
nmax_direct = 8*(nxhi_direct+1)*(nyhi_direct+1)*(nzhi_direct+1);
deallocate();
if (peratom_allocate_flag) deallocate_peratom();
// compute direct sum interaction weights
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);
}
}
if (!triclinic)
boxlo = domain->boxlo;
else
boxlo = domain->boxlo_lamda;
// ghost grid points depend on direct sum interaction limits, so need to recompute local grid
set_grid_local();
// allocate K-space dependent memory
// don't invoke allocate_peratom(), compute() will allocate when needed
allocate();
peratom_allocate_flag = 0;
// setup commgrid
cg_all->ghost_notify();
cg_all->setup();
for (int n=0; n<levels; n++) {
if (!active_flag[n]) continue;
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
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = evflag_atom = eflag_global = vflag_global =
eflag_atom = vflag_atom = eflag_either = vflag_either = 0;
// invoke allocate_peratom() if needed for first time
if (vflag_atom && !peratom_allocate_flag) {
allocate_peratom();
cg_peratom_all->ghost_notify();
cg_peratom_all->setup();
for (int n=0; n<levels; n++) {
if (!active_flag[n]) continue;
cg_peratom[n]->ghost_notify();
cg_peratom[n]->setup();
}
peratom_allocate_flag = 1;
}
// convert atoms from box to lamda coords
if (triclinic)
domain->x2lamda(atom->nlocal);
// 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();
// all procs reverse communicate charge density values from their ghost grid points
// to fully sum contribution in their 3d grid
current_level = 0;
cg_all->reverse_comm(this,REVERSE_RHO);
// forward communicate charge density values to fill ghost grid points
// compute direct sum interaction and then restrict to coarser grid
for (int n=0; n<=levels-2; n++) {
if (!active_flag[n]) continue;
current_level = n;
cg[n]->forward_comm(this,FORWARD_RHO);
direct(n);
restriction(n);
}
// compute direct interation for top grid level for nonperiodic
// and for second from top grid level for periodic
if (active_flag[levels-1]) {
if (domain->nonperiodic) {
current_level = levels-1;
cg[levels-1]->forward_comm(this,FORWARD_RHO);
direct_top(levels-1);
cg[levels-1]->reverse_comm(this,REVERSE_AD);
if (vflag_atom)
cg_peratom[levels-1]->reverse_comm(this,REVERSE_AD_PERATOM);
} else {
// Here using MPI_Allreduce is cheaper than using commgrid
grid_swap_forward(levels-1,qgrid[levels-1]);
direct(levels-1);
grid_swap_reverse(levels-1,egrid[levels-1]);
current_level = levels-1;
if (vflag_atom)
cg_peratom[levels-1]->reverse_comm(this,REVERSE_AD_PERATOM);
}
}
// prolongate energy/virial from coarser grid to finer grid
// reverse communicate from ghost grid points to get full sum
for (int n=levels-2; n>=0; n--) {
if (!active_flag[n]) continue;
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_all->forward_comm(this,FORWARD_AD);
// extra per-atom energy/virial communication
if (vflag_atom)
cg_peratom_all->forward_comm(this,FORWARD_AD_PERATOM);
// calculate the force on my particles (interpolation)
fieldforce();
// calculate the per-atom energy/virial for my particles
if (evflag_atom) fieldforce_peratom();
// total long-range energy
const double qscale = force->qqrd2e * scale;
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;
}
}
// convert atoms back from lamda to box coords
if (triclinic)
domain->lamda2x(atom->nlocal);
}
/* ----------------------------------------------------------------------
allocate memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::allocate()
{
// interpolation coeffs
memory->create2d_offset(phi1d,3,-order,order,"msm:phi1d");
memory->create2d_offset(dphi1d,3,-order,order,"msm:dphi1d");
// commgrid using all processors for finest grid level
int (*procneigh_all)[2] = comm->procneigh;
cg_all = new CommGrid(lmp,world,1,1,
nxlo_in[0],nxhi_in[0],nylo_in[0],nyhi_in[0],nzlo_in[0],nzhi_in[0],
nxlo_out_all,nxhi_out_all,nylo_out_all,nyhi_out_all,nzlo_out_all,nzhi_out_all,
nxlo_out[0],nxhi_out[0],nylo_out[0],nyhi_out[0],nzlo_out[0],nzhi_out[0],
procneigh_all[0][0],procneigh_all[0][1],procneigh_all[1][0],
procneigh_all[1][1],procneigh_all[2][0],procneigh_all[2][1]);
// allocate memory for each grid level
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 commgrid object for rho and electric field communication
if (active_flag[n]) {
int **procneigh = procneigh_levels[n];
cg[n] = new CommGrid(lmp,world_levels[n],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 virial memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::allocate_peratom()
{
// create commgrid object for per-atom virial using all processors
int (*procneigh_all)[2] = comm->procneigh;
cg_peratom_all =
new CommGrid(lmp,world,6,6,
nxlo_in[0],nxhi_in[0],nylo_in[0],nyhi_in[0],nzlo_in[0],nzhi_in[0],
nxlo_out_all,nxhi_out_all,nylo_out_all,nyhi_out_all,nzlo_out_all,nzhi_out_all,
nxlo_out[0],nxhi_out[0],nylo_out[0],nyhi_out[0],nzlo_out[0],nzhi_out[0],
procneigh_all[0][0],procneigh_all[0][1],procneigh_all[1][0],
procneigh_all[1][1],procneigh_all[2][0],procneigh_all[2][1]);
// allocate memory for each grid level
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 commgrid object for per-atom virial
if (active_flag[n]) {
int **procneigh = procneigh_levels[n];
cg_peratom[n] =
new CommGrid(lmp,world_levels[n],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);
if (cg_all) delete cg_all;
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 (world_levels)
if (world_levels[n] != MPI_COMM_NULL)
MPI_Comm_free(&world_levels[n]);
if (cg)
if (cg[n]) delete cg[n];
}
}
/* ----------------------------------------------------------------------
deallocate per-atom virial memory that depends on # of grid points
------------------------------------------------------------------------- */
void MSM::deallocate_peratom()
{
if (cg_peratom_all) delete cg_peratom_all;
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];
memory->create(procneigh_levels,levels,3,2,"msm:procneigh_levels");
world_levels = new MPI_Comm[levels];
active_flag = new int[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];
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;
world_levels[n] = MPI_COMM_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;
memory->destroy(procneigh_levels);
delete [] world_levels;
delete [] active_flag;
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 [] qgrid;
delete [] egrid;
delete [] v0grid;
delete [] v1grid;
delete [] v2grid;
delete [] v3grid;
delete [] v4grid;
delete [] v5grid;
}
/* ----------------------------------------------------------------------
set total 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) {
// seek to choose optimal Coulombic cutoff and number of grid levels
// (based on a cost estimate in Hardy's thesis)
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(xprd*yprd*zprd/atom->natoms,1.0/3.0);
nx_max = static_cast<int>(xprd/hmin);
ny_max = static_cast<int>(yprd/hmin);
nz_max = static_cast<int>(zprd/hmin);
nx_max = MAX(nx_max,2);
ny_max = MAX(ny_max,2);
nz_max = MAX(nz_max,2);
} else if (!gridflag) {
// Coulombic cutoff is set by user, choose grid to give requested error
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 {
// cutoff and grid are set by user
nx_max = nx_msm_max;
ny_max = ny_msm_max;
nz_max = nz_msm_max;
}
// adjust Coulombic cutoff to give desired error (if requested)
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_coul",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);
}
// scale grid for triclinic skew
if (triclinic && !gridflag) {
double tmp[3];
tmp[0] = nx_max/xprd;
tmp[1] = ny_max/yprd;
tmp[2] = nz_max/zprd;
lamda2xT(&tmp[0],&tmp[0]);
nx_max = static_cast<int>(tmp[0]);
ny_max = static_cast<int>(tmp[1]);
nz_max = static_cast<int>(tmp[2]);
}
// boost grid size until it is factorable by 2
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");
if (triclinic == 0) {
h_x = xprd/nx_max;
h_y = yprd/ny_max;
h_z = zprd/nz_max;
} else {
double tmp[3];
tmp[0] = nx_max;
tmp[1] = ny_max;
tmp[2] = nz_max;
x2lamdaT(&tmp[0],&tmp[0]);
h_x = 1.0/tmp[0];
h_y = 1.0/tmp[1];
h_z = 1.0/tmp[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");
}
// omit top grid level for periodic systems
if (!domain->nonperiodic) levels -= 1;
allocate_levels();
// find number of grid levels in each direction
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");
// compute number of extra grid points needed for nonperiodic boundary conditions
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);
}
}
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 grid that I own (inclusive)
n xyz lo/hi out = 3d grid + ghost cells in 6 directions (inclusive)
------------------------------------------------------------------------- */
void MSM::set_grid_local()
{
// 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 (interpolating) particles to MSM grid
nlower = -(order-1)/2;
nupper = order/2;
// lengths of box and processor sub-domains
double *prd,*sublo,*subhi,*boxhi;
if (!triclinic) {
prd = domain->prd;
sublo = domain->sublo;
subhi = domain->subhi;
} else {
prd = domain->prd_lamda;
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
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
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;
if (triclinic) kspacebbox(cuthalf,&dist[0]);
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;
if (n == 0) {
// use a smaller ghost region for interpolation
nxlo_out_all = nlo + nlower;
nxhi_out_all = nhi + nupper;
}
// a larger ghost region is needed for the direct sum and for restriction/prolongation
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;
if (n == 0) {
nylo_out_all = nlo + nlower;
nyhi_out_all = nhi + nupper;
}
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;
if (n == 0) {
nzlo_out_all = nlo + nlower;
nzhi_out_all = nhi + nupper;
}
// a hemisphere is used for direct sum interactions,
// so no ghosting is needed for direct sum in the -z direction
nzlo_out[n] = nlo - order;
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 (n == 0) nxlo_out_all = MAX(nxlo_out_all,alpha[0]);
if (nxhi_in[n] == nx_msm[n] - 1)
nxhi_in[n] = betax[n];
nxhi_out[n] = MIN(nxhi_out[n],betax[n]);
if (n == 0) nxhi_out_all = MIN(nxhi_out_all,betax[0]);
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 (n == 0) nylo_out_all = MAX(nylo_out_all,alpha[0]);
if (nyhi_in[n] == ny_msm[n] - 1)
nyhi_in[n] = betay[n];
nyhi_out[n] = MIN(nyhi_out[n],betay[n]);
if (n == 0) nyhi_out_all = MIN(nyhi_out_all,betay[0]);
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 (n == 0) nzlo_out_all = MAX(nzlo_out_all,alpha[0]);
if (nzhi_in[n] == nz_msm[n] - 1)
nzhi_in[n] = betaz[n];
nzhi_out[n] = MIN(nzhi_out[n],betaz[n]);
if (n == 0) nzhi_out_all = MIN(nzhi_out_all,betaz[0]);
if (nzhi_in[n] < 0)
nzhi_in[n] = alpha[n] - 1;
}
}
// prevent inactive processors from participating in MPI communication routines
set_proc_grid(n);
// 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);
}
}
/* ----------------------------------------------------------------------
find active procs and neighbor procs for each grid level
------------------------------------------------------------------------- */
void MSM::set_proc_grid(int n)
{
for (int i=0; i<3; i++)
myloc[i] = comm->myloc[i];
// size of inner MSM grid owned by this proc
int nxgrid_in = nxhi_in[n]-nxlo_in[n]+1;
int nygrid_in = nyhi_in[n]-nylo_in[n]+1;
int nzgrid_in = nzhi_in[n]-nzlo_in[n]+1;
int ngrid_in = nxgrid_in * nygrid_in * nzgrid_in;
// check to see if this proc owns any inner grid points on this grid level
// if not, deactivate by setting active_flag = 0
int cnt[3];
cnt[0] = 0;
if (myloc[1] == 0 && myloc[2] == 0)
if (nxgrid_in > 0)
cnt[0] = 1;
cnt[1] = 0;
if (myloc[0] == 0 && myloc[2] == 0)
if (nygrid_in > 0)
cnt[1] = 1;
cnt[2] = 0;
if (myloc[0] == 0 && myloc[1] == 0)
if (nzgrid_in > 0)
cnt[2] = 1;
MPI_Allreduce(&cnt[0],&procgrid[0],3,MPI_INT,MPI_SUM,world);
int color;
if (ngrid_in > 0) {
active_flag[n] = 1;
color = 0;
} else {
active_flag[n] = 0;
color = MPI_UNDEFINED;
}
// define a new MPI communicator for this grid level that only includes active procs
MPI_Comm_split(world,color,me,&world_levels[n]);
if (!active_flag[n]) return;
int procneigh[3][2]; // my 6 neighboring procs, 0/1 = left/right
// map processor IDs to new 3d processor grid
// produces myloc, procneigh
int periods[3];
periods[0] = periods[1] = periods[2] = 1;
MPI_Comm cartesian;
MPI_Cart_create(world_levels[n],3,procgrid,periods,0,&cartesian);
MPI_Cart_get(cartesian,3,procgrid,periods,myloc);
MPI_Cart_shift(cartesian,0,1,&procneigh[0][0],&procneigh[0][1]);
MPI_Cart_shift(cartesian,1,1,&procneigh[1][0],&procneigh[1][1]);
MPI_Cart_shift(cartesian,2,1,&procneigh[2][0],&procneigh[2][1]);
MPI_Comm_free(&cartesian);
for (int i=0; i<3; i++)
for (int j=0; j<2; j++)
procneigh_levels[n][i][j] = procneigh[i][j];
}
/* ----------------------------------------------------------------------
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;
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");
}
/* ----------------------------------------------------------------------
aninterpolation: interpolate charges from particles to grid
------------------------------------------------------------------------- */
void MSM::make_rho()
{
//fprintf(screen,"MSM aninterpolation\n\n");
int i,l,m,n,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(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 sum procedure for intermediate grid levels, solve Poisson's
equation to get energy, virial, etc.
------------------------------------------------------------------------- */
void MSM::direct(int n)
{
//fprintf(screen,"Direct contribution on level %i\n\n",n);
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];
double *g_directn = g_direct[n];
double *v0_directn = v0_direct[n];
double *v1_directn = v1_direct[n];
double *v2_directn = v2_direct[n];
double *v3_directn = v3_direct[n];
double *v4_directn = v4_direct[n];
double *v5_directn = v5_direct[n];
// zero out electric potential
memset(&(egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
// 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 ii,jj,kk;
int imin,imax,jmin,jmax,kmin,kmax;
double qtmp,qtmp2,gtmp;
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
double **qk,**ek;
double *qkj,*ekj;
int nx = nxhi_direct - nxlo_direct + 1;
int ny = nyhi_direct - nylo_direct + 1;
// loop over inner grid points
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);
}
qtmp = qgridn[icz][icy][icx]; // charge on center grid point
esum = 0.0;
if (vflag_either)
v0sum = v1sum = v2sum = v3sum = v4sum = v5sum = 0.0;
// use hemisphere to avoid double computation of pair-wise
// interactions in direct sum (no computations in -z direction)
for (iz = 1; iz <= kmax; iz++) {
kk = icz+iz;
qk = qgridn[kk];
ek = egridn[kk];
zk = (iz + nzhi_direct)*ny;
for (iy = jmin; iy <= jmax; iy++) {
jj = icy+iy;
qkj = qk[jj];
ekj = ek[jj];
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = imin; ix <= imax; ix++) {
ii = icx+ix;
qtmp2 = qkj[ii]; // charge on outer grid point
k = zyk + ix + nxhi_direct;
gtmp = g_directn[k];
esum += gtmp * qtmp2;
ekj[ii] += gtmp * qtmp;
if (vflag_either) {
v0sum += v0_directn[k] * qtmp2;
v1sum += v1_directn[k] * qtmp2;
v2sum += v2_directn[k] * qtmp2;
v3sum += v3_directn[k] * qtmp2;
v4sum += v4_directn[k] * qtmp2;
v5sum += v5_directn[k] * qtmp2;
}
}
}
}
// iz=0
iz = 0;
kk = icz+iz;
qk = qgridn[kk];
ek = egridn[kk];
zk = (iz + nzhi_direct)*ny;
for (iy = 1; iy <= jmax; iy++) {
jj = icy+iy;
qkj = qk[jj];
ekj = ek[jj];
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = imin; ix <= imax; ix++) {
ii = icx+ix;
qtmp2 = qkj[ii];
k = zyk + ix + nxhi_direct;
gtmp = g_directn[k];
esum += gtmp * qtmp2;
ekj[ii] += gtmp * qtmp;
if (vflag_either) {
v0sum += v0_directn[k] * qtmp2;
v1sum += v1_directn[k] * qtmp2;
v2sum += v2_directn[k] * qtmp2;
v3sum += v3_directn[k] * qtmp2;
v4sum += v4_directn[k] * qtmp2;
v5sum += v5_directn[k] * qtmp2;
}
}
}
// iz=0, iy=0
iz = 0;
kk = icz+iz;
qk = qgridn[kk];
ek = egridn[kk];
zk = (iz + nzhi_direct)*ny;
iy = 0;
jj = icy+iy;
qkj = qk[jj];
ekj = ek[jj];
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = 1; ix <= imax; ix++) {
ii = icx+ix;
qtmp2 = qkj[ii];
k = zyk + ix + nxhi_direct;
gtmp = g_directn[k];
esum += gtmp * qtmp2;
ekj[ii] += gtmp * qtmp;
if (vflag_either) {
v0sum += v0_directn[k] * qtmp2;
v1sum += v1_directn[k] * qtmp2;
v2sum += v2_directn[k] * qtmp2;
v3sum += v3_directn[k] * qtmp2;
v4sum += v4_directn[k] * qtmp2;
v5sum += v5_directn[k] * qtmp2;
}
}
// iz=0, iy=0, ix=0
iz = 0;
zk = (iz + nzhi_direct)*ny;
iy = 0;
zyk = (zk + iy + nyhi_direct)*nx;
ix = 0;
k = zyk + ix + nxhi_direct;
gtmp = g_directn[k];
esum += 0.5 * gtmp * qtmp;
egridn[icz][icy][icx] += 0.5 * gtmp * qtmp;
// virial is zero for iz=0, iy=0, ix=0
// accumulate per-atom energy/virial
egridn[icz][icy][icx] += esum;
if (vflag_atom) {
v0gridn[icz][icy][icx] += v0sum;
v1gridn[icz][icy][icx] += v1sum;
v2gridn[icz][icy][icx] += v2sum;
v3gridn[icz][icy][icx] += v3sum;
v4gridn[icz][icy][icx] += v4sum;
v5gridn[icz][icy][icx] += v5sum;
}
// accumulate total energy/virial
if (evflag) {
qtmp = qgridn[icz][icy][icx];
if (eflag_global) energy += 2.0 * esum * qtmp;
if (vflag_global) {
virial[0] += 2.0 * v0sum * qtmp;
virial[1] += 2.0 * v1sum * qtmp;
virial[2] += 2.0 * v2sum * qtmp;
virial[3] += 2.0 * v3sum * qtmp;
virial[4] += 2.0 * v4sum * qtmp;
virial[5] += 2.0 * v5sum * qtmp;
}
}
}
}
}
// compute per-atom virial (if requested)
if (vflag_atom)
direct_peratom(n);
}
/* ----------------------------------------------------------------------
MSM direct sum procedure for intermediate grid levels, solve Poisson's
equation to get per-atom virial, separate method used for performance
reasons
------------------------------------------------------------------------- */
void MSM::direct_peratom(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];
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
int ii,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;
// loop over inner grid points
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);
}
qtmp = qgridn[icz][icy][icx]; // center grid point
// use hemisphere to avoid double computation of pair-wise
// interactions in direct sum (no computations in -z direction)
for (iz = 1; 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++) {
ii = icx+ix;
k = zyk + ix + nxhi_direct;
v0gridn[kk][jj][ii] += v0_direct[n][k] * qtmp;
v1gridn[kk][jj][ii] += v1_direct[n][k] * qtmp;
v2gridn[kk][jj][ii] += v2_direct[n][k] * qtmp;
v3gridn[kk][jj][ii] += v3_direct[n][k] * qtmp;
v4gridn[kk][jj][ii] += v4_direct[n][k] * qtmp;
v5gridn[kk][jj][ii] += v5_direct[n][k] * qtmp;
}
}
}
// iz=0
iz = 0;
kk = icz+iz;
zk = (iz + nzhi_direct)*ny;
for (iy = 1; iy <= jmax; iy++) {
jj = icy+iy;
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = imin; ix <= imax; ix++) {
ii = icx+ix;
k = zyk + ix + nxhi_direct;
v0gridn[kk][jj][ii] += v0_direct[n][k] * qtmp;
v1gridn[kk][jj][ii] += v1_direct[n][k] * qtmp;
v2gridn[kk][jj][ii] += v2_direct[n][k] * qtmp;
v3gridn[kk][jj][ii] += v3_direct[n][k] * qtmp;
v4gridn[kk][jj][ii] += v4_direct[n][k] * qtmp;
v5gridn[kk][jj][ii] += v5_direct[n][k] * qtmp;
}
}
// iz=0, iy=0
iz = 0;
kk = icz+iz;
zk = (iz + nzhi_direct)*ny;
iy = 0;
jj = icy+iy;
zyk = (zk + iy + nyhi_direct)*nx;
for (ix = 1; ix <= imax; ix++) {
ii = icx+ix;
k = zyk + ix + nxhi_direct;
v0gridn[kk][jj][ii] += v0_direct[n][k] * qtmp;
v1gridn[kk][jj][ii] += v1_direct[n][k] * qtmp;
v2gridn[kk][jj][ii] += v2_direct[n][k] * qtmp;
v3gridn[kk][jj][ii] += v3_direct[n][k] * qtmp;
v4gridn[kk][jj][ii] += v4_direct[n][k] * qtmp;
v5gridn[kk][jj][ii] += v5_direct[n][k] * qtmp;
}
// virial is zero for iz=0, iy=0, ix=0
}
}
}
}
/* ----------------------------------------------------------------------
MSM direct sum procedure for top grid level (nonperiodic systems only)
------------------------------------------------------------------------- */
void MSM::direct_top(int n)
{
//fprintf(screen,"Direct contribution on level %i\n\n",n);
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];
// zero out electric potential
memset(&(egridn[nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
// 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 ii,jj,kk;
int imin,imax,jmin,jmax,kmin,kmax;
double qtmp,qtmp2,gtmp;
double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
double **qk,**ek;
double *qkj,*ekj;
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;
// loop over inner grid points
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
if (domain->zperiodic) {
kmin = 0;
kmax = nz_msm[n]-1;
} else {
kmin = alpha[n] - icz;
kmax = betaz[n] - icz;
}
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) {
if (domain->yperiodic) {
jmin = 0;
jmax = ny_msm[n]-1;
} else {
jmin = alpha[n] - icy;
jmax = betay[n] - icy;
}
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) {
if (domain->xperiodic) {
imin = 0;
imax = nx_msm[n]-1;
} else {
imin = alpha[n] - icx;
imax = betax[n] - icx;
}
qtmp = qgridn[icz][icy][icx];
esum = 0.0;
if (vflag_either)
v0sum = v1sum = v2sum = v3sum = v4sum = v5sum = 0.0;
// use hemisphere to avoid double computation of pair-wise
// interactions in direct sum (no computations in -z direction)
for (iz = 1; iz <= kmax; iz++) {
kk = icz+iz;
qk = qgridn[kk];
ek = egridn[kk];
zk = (iz + nz_top)*ny;
for (iy = jmin; iy <= jmax; iy++) {
jj = icy+iy;
qkj = qk[jj];
ekj = ek[jj];
zyk = (zk + iy + ny_top)*nx;
for (ix = imin; ix <= imax; ix++) {
ii = icx+ix;
qtmp2 = qkj[ii];
k = zyk + ix + nx_top;
gtmp = g_direct_top[k];
esum += gtmp * qtmp2;
ekj[ii] += gtmp * qtmp;
if (vflag_either) {
v0sum += v0_direct_top[k] * qtmp2;
v1sum += v1_direct_top[k] * qtmp2;
v2sum += v2_direct_top[k] * qtmp2;
v3sum += v3_direct_top[k] * qtmp2;
v4sum += v4_direct_top[k] * qtmp2;
v5sum += v5_direct_top[k] * qtmp2;
}
}
}
}
// iz=0
iz = 0;
kk = icz+iz;
qk = qgridn[kk];
ek = egridn[kk];
zk = (iz + nz_top)*ny;
for (iy = 1; iy <= jmax; iy++) {
jj = icy+iy;
qkj = qk[jj];
ekj = ek[jj];
zyk = (zk + iy + ny_top)*nx;
for (ix = imin; ix <= imax; ix++) {
ii = icx+ix;
qtmp2 = qkj[ii];
k = zyk + ix + nx_top;
gtmp = g_direct_top[k];
esum += gtmp * qtmp2;
ekj[ii] += gtmp * qtmp;
if (vflag_either) {
v0sum += v0_direct_top[k] * qtmp2;
v1sum += v1_direct_top[k] * qtmp2;
v2sum += v2_direct_top[k] * qtmp2;
v3sum += v3_direct_top[k] * qtmp2;
v4sum += v4_direct_top[k] * qtmp2;
v5sum += v5_direct_top[k] * qtmp2;
}
}
}
// iz=0, iy=0
iz = 0;
kk = icz+iz;
qk = qgridn[kk];
ek = egridn[kk];
zk = (iz + nz_top)*ny;
iy = 0;
jj = icy+iy;
qkj = qk[jj];
ekj = ek[jj];
zyk = (zk + iy + ny_top)*nx;
for (ix = 1; ix <= imax; ix++) {
ii = icx+ix;
qtmp2 = qkj[ii];
k = zyk + ix + nx_top;
gtmp = g_direct_top[k];
esum += gtmp * qtmp2;
ekj[ii] += gtmp * qtmp;
if (vflag_either) {
v0sum += v0_direct_top[k] * qtmp2;
v1sum += v1_direct_top[k] * qtmp2;
v2sum += v2_direct_top[k] * qtmp2;
v3sum += v3_direct_top[k] * qtmp2;
v4sum += v4_direct_top[k] * qtmp2;
v5sum += v5_direct_top[k] * qtmp2;
}
}
// iz=0, iy=0, ix=0
iz = 0;
zk = (iz + nz_top)*ny;
iy = 0;
zyk = (zk + iy + ny_top)*nx;
ix = 0;
ii = icx+ix;
k = zyk + ix + nx_top;
gtmp = g_direct_top[k];
esum += 0.5 * gtmp * qtmp;
egridn[icz][icy][icx] += 0.5 * gtmp * qtmp;
if (vflag_either) {
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;
}
// accumulate per-atom energy/virial
egridn[icz][icy][icx] += esum;
if (vflag_atom) {
v0gridn[icz][icy][icx] += v0sum;
v1gridn[icz][icy][icx] += v1sum;
v2gridn[icz][icy][icx] += v2sum;
v3gridn[icz][icy][icx] += v3sum;
v4gridn[icz][icy][icx] += v4sum;
v5gridn[icz][icy][icx] += v5sum;
}
// accumulate total energy/virial
if (evflag) {
qtmp = qgridn[icz][icy][icx];
if (eflag_global) energy += 2.0 * esum * qtmp;
if (vflag_global) {
virial[0] += 2.0 * v0sum * qtmp;
virial[1] += 2.0 * v1sum * qtmp;
virial[2] += 2.0 * v2sum * qtmp;
virial[3] += 2.0 * v3sum * qtmp;
virial[4] += 2.0 * v4sum * qtmp;
virial[5] += 2.0 * v5sum * qtmp;
}
}
}
}
}
// compute per-atom virial (if requested)
if (vflag_atom)
direct_peratom_top(n);
}
/* ----------------------------------------------------------------------
MSM direct sum procedure for top grid level, solve Poisson's
equation to get per-atom virial, separate method used for performance
reasons
------------------------------------------------------------------------- */
void MSM::direct_peratom_top(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];
int icx,icy,icz,ix,iy,iz,zk,zyk,k;
int ii,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;
// loop over inner grid points
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++) {
if (domain->zperiodic) {
kmin = 0;
kmax = nz_msm[n]-1;
} else {
kmin = alpha[n] - icz;
kmax = betaz[n] - icz;
}
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++) {
if (domain->yperiodic) {
jmin = 0;
jmax = ny_msm[n]-1;
} else {
jmin = alpha[n] - icy;
jmax = betay[n] - icy;
}
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++) {
if (domain->xperiodic) {
imin = 0;
imax = nx_msm[n]-1;
} else {
imin = alpha[n] - icx;
imax = betax[n] - icx;
}
qtmp = qgridn[icz][icy][icx]; // center grid point
// use hemisphere to avoid double computation of pair-wise
// interactions in direct sum (no computations in -z direction)
for (iz = 1; 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++) {
ii = icx+ix;
k = zyk + ix + nx_top;
v0gridn[kk][jj][ii] += v0_direct_top[k] * qtmp;
v1gridn[kk][jj][ii] += v1_direct_top[k] * qtmp;
v2gridn[kk][jj][ii] += v2_direct_top[k] * qtmp;
v3gridn[kk][jj][ii] += v3_direct_top[k] * qtmp;
v4gridn[kk][jj][ii] += v4_direct_top[k] * qtmp;
v5gridn[kk][jj][ii] += v5_direct_top[k] * qtmp;
}
}
}
// iz=0
iz = 0;
kk = icz+iz;
zk = (iz + nz_top)*ny;
for (iy = 1; iy <= jmax; iy++) {
jj = icy+iy;
zyk = (zk + iy + ny_top)*nx;
for (ix = imin; ix <= imax; ix++) {
ii = icx+ix;
k = zyk + ix + nx_top;
v0gridn[kk][jj][ii] += v0_direct_top[k] * qtmp;
v1gridn[kk][jj][ii] += v1_direct_top[k] * qtmp;
v2gridn[kk][jj][ii] += v2_direct_top[k] * qtmp;
v3gridn[kk][jj][ii] += v3_direct_top[k] * qtmp;
v4gridn[kk][jj][ii] += v4_direct_top[k] * qtmp;
v5gridn[kk][jj][ii] += v5_direct_top[k] * qtmp;
}
}
// iz=0, iy=0
iz = 0;
kk = icz+iz;
zk = (iz + nz_top)*ny;
iy = 0;
jj = icy+iy;
zyk = (zk + iy + ny_top)*nx;
for (ix = 1; ix <= imax; ix++) {
ii = icx+ix;
k = zyk + ix + nx_top;
v0gridn[kk][jj][ii] += v0_direct_top[k] * qtmp;
v1gridn[kk][jj][ii] += v1_direct_top[k] * qtmp;
v2gridn[kk][jj][ii] += v2_direct_top[k] * qtmp;
v3gridn[kk][jj][ii] += v3_direct_top[k] * qtmp;
v4gridn[kk][jj][ii] += v4_direct_top[k] * qtmp;
v5gridn[kk][jj][ii] += v5_direct_top[k] * qtmp;
}
// virial is zero for iz=0, iy=0, ix=0
}
}
}
}
/* ----------------------------------------------------------------------
MSM restriction procedure for intermediate grid levels, interpolate
charges from finer grid to coarser grid
------------------------------------------------------------------------- */
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];
int k = 0;
int index[p+2];
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,q2sum;
// 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]);
q2sum = 0.0;
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;
}
q2sum += qgrid1[kk][jj][ii] *
phi1d[0][i]*phizy;
}
}
}
qgrid2[kp][jp][ip] += q2sum;
}
}
/* ----------------------------------------------------------------------
MSM prolongation procedure for intermediate grid levels, interpolate
per-atom energy/virial from coarser grid to finer grid
------------------------------------------------------------------------- */
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];
int k = 0;
int index[p+2];
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;
double etmp2,v0tmp2,v1tmp2,v2tmp2,v3tmp2,v4tmp2,v5tmp2;
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]);
etmp2 = egrid2[kp][jp][ip];
if (vflag_atom) {
v0tmp2 = v0grid2[kp][jp][ip];
v1tmp2 = v1grid2[kp][jp][ip];
v2tmp2 = v2grid2[kp][jp][ip];
v3tmp2 = v3grid2[kp][jp][ip];
v4tmp2 = v4grid2[kp][jp][ip];
v5tmp2 = v5grid2[kp][jp][ip];
}
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] += etmp2 * phi3d;
if (vflag_atom) {
v0grid1[kk][jj][ii] += v0tmp2 * phi3d;
v1grid1[kk][jj][ii] += v1tmp2 * phi3d;
v2grid1[kk][jj][ii] += v2tmp2 * phi3d;
v3grid1[kk][jj][ii] += v3tmp2 * phi3d;
v4grid1[kk][jj][ii] += v4tmp2 * phi3d;
v5grid1[kk][jj][ii] += v5tmp2 * phi3d;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
Use MPI_Allreduce to fill ghost grid values, for coarse grids this may
be cheaper than using nearest-neighbor communication (commgrid), right
now only works for periodic boundary conditions
------------------------------------------------------------------------- */
void MSM::grid_swap_forward(int n, double*** &gridn)
{
double ***gridn_tmp;
memory->create(gridn_tmp,nz_msm[n],ny_msm[n],nx_msm[n],"msm:grid_tmp");
double ***gridn_all;
memory->create(gridn_all,nz_msm[n],ny_msm[n],nx_msm[n],"msm:grid_all");
int ngrid_in = nx_msm[n] * ny_msm[n] * nz_msm[n];
memset(&(gridn_tmp[0][0][0]),0,ngrid_in*sizeof(double));
memset(&(gridn_all[0][0][0]),0,ngrid_in*sizeof(double));
// copy inner grid cell values from gridn into gridn_tmp
int icx,icy,icz;
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++)
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++)
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++)
gridn_tmp[icz][icy][icx] = gridn[icz][icy][icx];
MPI_Allreduce(&(gridn_tmp[0][0][0]),
&(gridn_all[0][0][0]),
ngrid_in,MPI_DOUBLE,MPI_SUM,world_levels[n]);
// bitmask for PBCs (only works for power of 2 numbers)
int PBCx,PBCy,PBCz;
PBCx = nx_msm[n]-1;
PBCy = ny_msm[n]-1;
PBCz = nz_msm[n]-1;
// copy from gridn_all into gridn to fill ghost grid cell values
for (icz = nzlo_out[n]; icz <= nzhi_out[n]; icz++)
for (icy = nylo_out[n]; icy <= nyhi_out[n]; icy++)
for (icx = nxlo_out[n]; icx <= nxhi_out[n]; icx++)
gridn[icz][icy][icx] = gridn_all[icz&PBCz][icy&PBCy][icx&PBCx];
memory->destroy(gridn_tmp);
memory->destroy(gridn_all);
}
/* ----------------------------------------------------------------------
Use MPI_Allreduce to get contribution from ghost grid cells, for coarse
grids this may be cheaper than using nearest-neighbor communication
(commgrid), right now only works for periodic boundary conditions
------------------------------------------------------------------------- */
void MSM::grid_swap_reverse(int n, double*** &gridn)
{
double ***gridn_tmp;
memory->create(gridn_tmp,nz_msm[n],ny_msm[n],nx_msm[n],"msm:grid_tmp");
double ***gridn_all;
memory->create(gridn_all,nz_msm[n],ny_msm[n],nx_msm[n],"msm:grid_all");
int ngrid_in = nx_msm[n] * ny_msm[n] * nz_msm[n];
memset(&(gridn_tmp[0][0][0]),0,ngrid_in*sizeof(double));
memset(&(gridn_all[0][0][0]),0,ngrid_in*sizeof(double));
// bitmask for PBCs (only works for power of 2 numbers)
int icx,icy,icz;
int PBCx,PBCy,PBCz;
PBCx = nx_msm[n]-1;
PBCy = ny_msm[n]-1;
PBCz = nz_msm[n]-1;
// copy ghost grid cell values from gridn into inner portion of gridn_tmp
for (icz = nzlo_out[n]; icz <= nzhi_out[n]; icz++)
for (icy = nylo_out[n]; icy <= nyhi_out[n]; icy++)
for (icx = nxlo_out[n]; icx <= nxhi_out[n]; icx++)
gridn_tmp[icz&PBCz][icy&PBCy][icx&PBCx] += gridn[icz][icy][icx];
MPI_Allreduce(&(gridn_tmp[0][0][0]),
&(gridn_all[0][0][0]),
ngrid_in,MPI_DOUBLE,MPI_SUM,world_levels[n]);
// copy inner grid cell values from gridn_all into gridn
for (icz = nzlo_in[n]; icz <= nzhi_in[n]; icz++)
for (icy = nylo_in[n]; icy <= nyhi_in[n]; icy++)
for (icx = nxlo_in[n]; icx <= nxhi_in[n]; icx++)
gridn[icz][icy][icx] = gridn_all[icz][icy][icx];
memory->destroy(gridn_tmp);
memory->destroy(gridn_all);
}
/* ----------------------------------------------------------------------
pack own values to buf to send to another proc (used by commgrid)
------------------------------------------------------------------------- */
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,etmp;
// 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];
etmp = egridn[mz][my][mx];
ekx += dphi_x*phi_y*phi_z*etmp;
eky += phi_x*dphi_y*phi_z*etmp;
ekz += phi_x*phi_y*dphi_z*etmp;
}
}
}
ekx *= delxinv[0];
eky *= delyinv[0];
ekz *= delzinv[0];
// effectively divide by length for a triclinic system
if (triclinic) {
double tmp[3];
tmp[0] = ekx;
tmp[1] = eky;
tmp[2] = ekz;
x2lamdaT(&tmp[0],&tmp[0]);
ekx = tmp[0];
eky = tmp[1];
ekz = tmp[2];
}
// 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 (interpolation coefficients)
------------------------------------------------------------------------- */
void MSM::compute_phis(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);
}
}
/* ----------------------------------------------------------------------
charge assignment into phi1d and dphi1d (interpolation coefficients)
------------------------------------------------------------------------- */
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
------------------------------------------------------------------------- */
inline 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
------------------------------------------------------------------------- */
inline 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_xi3 = xi2*abs_xi;
double abs_xi5 = xi4*abs_xi;
if (abs_xi == 0.0) {
dphi = 0.0;
} else if (abs_xi <= 1) {
dphi = xi*(49*xi6 - 175*xi4 + 84*xi2 - 150*abs_xi5 +
644*abs_xi3 - 560*abs_xi)/144.0/abs_xi;
} else if (abs_xi <= 2) {
dphi = xi*(-49*xi6 - 1365*xi4 + 756*xi2 +
450*abs_xi5 + 1260*abs_xi3 - 1260*abs_xi + 28)/240.0/abs_xi;
} else if (abs_xi <= 3) {
dphi = xi*(49*xi6 + 4445*xi4 + 17724*xi2 -
750*abs_xi5 - 12740*abs_xi3 - 9940*abs_xi + 756)/720.0/abs_xi;
} else if (abs_xi <= 4) {
dphi = -xi*(abs_xi - 4)*(7*abs_xi5 - 122*xi4 +
807*abs_xi3 - 2512*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 (energy) weights 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 dx,dy,dz;
double tmp[3];
double rsq,rho,two_n;
two_n = 1.0;
int nx = nxhi_direct - nxlo_direct + 1;
int ny = nyhi_direct - nylo_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];
// transform grid point pair-wise distance from lamda (0-1) coords to box coords
if (triclinic) {
tmp[0] = xdiff;
tmp[1] = ydiff;
tmp[2] = zdiff;
lamda2xvector(&tmp[0],&tmp[0]);
dx = tmp[0];
dy = tmp[1];
dz = tmp[2];
} else {
dx = xdiff;
dy = ydiff;
dz = zdiff;
}
rsq = dx*dx + dy*dy + dz*dz;
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 (virial) weights 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 dx,dy,dz;
double tmp[3];
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;
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];
if (triclinic) {
tmp[0] = xdiff;
tmp[1] = ydiff;
tmp[2] = zdiff;
lamda2xvector(&tmp[0],&tmp[0]);
dx = tmp[0];
dy = tmp[1];
dz = tmp[2];
} else {
dx = xdiff;
dy = ydiff;
dz = zdiff;
}
rsq = dx*dx + dy*dy + dz*dz;
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 * dx * dx;
v1_direct[n][k] = dg * dy * dy;
v2_direct[n][k] = dg * dz * dz;
v3_direct[n][k] = dg * dx * dy;
v4_direct[n][k] = dg * dx * dz;
v5_direct[n][k] = dg * dy * dz;
}
}
}
}
two_n *= 2.0;
}
}
/* ----------------------------------------------------------------------
Compute direct interaction (energy) weights for top grid level
(nonperiodic systems only)
------------------------------------------------------------------------- */
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 dx,dy,dz;
double tmp[3];
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];
if (triclinic) {
tmp[0] = xdiff;
tmp[1] = ydiff;
tmp[2] = zdiff;
lamda2xvector(&tmp[0],&tmp[0]);
dx = tmp[0];
dy = tmp[1];
dz = tmp[2];
} else {
dx = xdiff;
dy = ydiff;
dz = zdiff;
}
rsq = dx*dx + dy*dy + dz*dz;
rho = sqrt(rsq)/(two_n*a);
k = zyk + ix + nx_top;
g_direct_top[k] = gamma(rho)/(two_n*a);
}
}
}
}
/* ----------------------------------------------------------------------
Compute direct interaction (virial) weights for top grid level
(nonperiodic systems only)
------------------------------------------------------------------------- */
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 dx,dy,dz;
double tmp[3];
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];
if (triclinic) {
tmp[0] = xdiff;
tmp[1] = ydiff;
tmp[2] = zdiff;
lamda2xvector(&tmp[0],&tmp[0]);
dx = tmp[0];
dy = tmp[1];
dz = tmp[2];
} else {
dx = xdiff;
dy = ydiff;
dz = zdiff;
}
rsq = dx*dx + dy*dy + dz*dz;
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 * dx * dx;
v1_direct_top[k] = dg * dy * dy;
v2_direct_top[k] = dg * dz * dz;
v3_direct_top[k] = dg * dx * dy;
v4_direct_top[k] = dg * dx * dz;
v5_direct_top[k] = dg * dy * dz;
}
}
}
}
}
Event Timeline
Log In to Comment