forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9401 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
a480021fa6
commit
0beb95626f
|
|
@ -6,6 +6,7 @@ if (test $1 = 1) then
|
|||
cp ewald.cpp ..
|
||||
cp ewald_disp.cpp ..
|
||||
cp msm.cpp ..
|
||||
cp msm_cg.cpp ..
|
||||
cp pppm.cpp ..
|
||||
cp pppm_old.cpp ..
|
||||
cp pppm_cg.cpp ..
|
||||
|
|
@ -35,6 +36,7 @@ if (test $1 = 1) then
|
|||
cp ewald.h ..
|
||||
cp ewald_disp.h ..
|
||||
cp msm.h ..
|
||||
cp msm_cg.h ..
|
||||
cp pppm.h ..
|
||||
cp pppm_old.h ..
|
||||
cp pppm_cg.h ..
|
||||
|
|
@ -69,6 +71,7 @@ elif (test $1 = 0) then
|
|||
rm -f ../ewald.cpp
|
||||
rm -f ../ewald_disp.cpp
|
||||
rm -f ../msm.cpp
|
||||
rm -f ../msm_cg.cpp
|
||||
rm -f ../pppm.cpp
|
||||
rm -f ../pppm_old.cpp
|
||||
rm -f ../pppm_cg.cpp
|
||||
|
|
@ -98,6 +101,7 @@ elif (test $1 = 0) then
|
|||
rm -f ../ewald.h
|
||||
rm -f ../ewald_disp.h
|
||||
rm -f ../msm.h
|
||||
rm -f ../msm_cg.h
|
||||
rm -f ../pppm.h
|
||||
rm -f ../pppm_old.h
|
||||
rm -f ../pppm_cg.h
|
||||
|
|
|
|||
|
|
@ -23,6 +23,7 @@
|
|||
#include "ewald_disp.h"
|
||||
#include "math_vector.h"
|
||||
#include "math_const.h"
|
||||
#include "math_special.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
|
|
@ -33,6 +34,7 @@
|
|||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using namespace MathSpecial;
|
||||
|
||||
#define SMALL 0.00001
|
||||
|
||||
|
|
@ -418,7 +420,7 @@ void EwaldDisp::coefficients()
|
|||
if (func12) { // -Bij/r^6 coeffs
|
||||
b1 = sqrt(b2); // minus sign folded
|
||||
h1 = sqrt(h2); // into constants
|
||||
*(ke++) = c1 = -h1*h2*((c2=sqrt(MY_PI)*erfc(b1))+(0.5/b2-1.0)*expb2/b1);
|
||||
*(ke++) = c1 = -h1*h2*((c2=MY_PIS*erfc(b1))+(0.5/b2-1.0)*expb2/b1);
|
||||
*(kv++) = c1-(c2 = 3.0*h1*(c2-expb2/b1))*h[0]*h[0];
|
||||
*(kv++) = c1-c2*h[1]*h[1]; // lammps convention
|
||||
*(kv++) = c1-c2*h[2]*h[2]; // instead of voigt
|
||||
|
|
@ -517,20 +519,20 @@ void EwaldDisp::init_self()
|
|||
|
||||
if (function[0]) { // 1/r
|
||||
virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*sum[0].x*sum[0].x;
|
||||
energy_self[0] = sum[0].x2*qscale*g1/sqrt(MY_PI)-virial_self[0];
|
||||
energy_self[0] = sum[0].x2*qscale*g1/MY_PIS-virial_self[0];
|
||||
}
|
||||
if (function[1]) { // geometric 1/r^6
|
||||
virial_self[1] = MY_PI*sqrt(MY_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x;
|
||||
virial_self[1] = MY_PI*MY_PIS*g3/(6.0*volume)*sum[1].x*sum[1].x;
|
||||
energy_self[1] = -sum[1].x2*g3*g3/12.0+virial_self[1];
|
||||
}
|
||||
if (function[2]) { // arithmetic 1/r^6
|
||||
virial_self[2] = MY_PI*sqrt(MY_PI)*g3/(48.0*volume)*(sum[2].x*sum[8].x+
|
||||
virial_self[2] = MY_PI*MY_PIS*g3/(48.0*volume)*(sum[2].x*sum[8].x+
|
||||
sum[3].x*sum[7].x+sum[4].x*sum[6].x+0.5*sum[5].x*sum[5].x);
|
||||
energy_self[2] = -sum[2].x2*g3*g3/3.0+virial_self[2];
|
||||
}
|
||||
if (function[3]) { // dipole
|
||||
virial_self[3] = 0; // in surface
|
||||
energy_self[3] = sum[9].x2*mumurd2e*2.0*g3/3.0/sqrt(MY_PI)-virial_self[3];
|
||||
energy_self[3] = sum[9].x2*mumurd2e*2.0*g3/3.0/MY_PIS-virial_self[3];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -551,7 +553,7 @@ void EwaldDisp::init_self_peratom()
|
|||
if (function[0]) { // 1/r
|
||||
double *ei = energy;
|
||||
double *vi = virial;
|
||||
double ce = qscale*g1/sqrt(MY_PI);
|
||||
double ce = qscale*g1/MY_PIS;
|
||||
double cv = -0.5*MY_PI*qscale/(g2*volume);
|
||||
double *qi = atom->q, *qn = qi + nlocal;
|
||||
for (; qi < qn; qi++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) {
|
||||
|
|
@ -564,7 +566,7 @@ void EwaldDisp::init_self_peratom()
|
|||
double *ei = energy+1;
|
||||
double *vi = virial+1;
|
||||
double ce = -g3*g3/12.0;
|
||||
double cv = MY_PI*sqrt(MY_PI)*g3/(6.0*volume);
|
||||
double cv = MY_PI*MY_PIS*g3/(6.0*volume);
|
||||
int *typei = atom->type, *typen = typei + atom->nlocal;
|
||||
for (; typei < typen; typei++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) {
|
||||
double b = B[*typei];
|
||||
|
|
@ -577,7 +579,7 @@ void EwaldDisp::init_self_peratom()
|
|||
double *ei = energy+2;
|
||||
double *vi = virial+2;
|
||||
double ce = -g3*g3/3.0;
|
||||
double cv = 0.5*MY_PI*sqrt(MY_PI)*g3/(48.0*volume);
|
||||
double cv = 0.5*MY_PI*MY_PIS*g3/(48.0*volume);
|
||||
int *typei = atom->type, *typen = typei + atom->nlocal;
|
||||
for (; typei < typen; typei++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) {
|
||||
bi = B+7*typei[0]+7;
|
||||
|
|
@ -601,7 +603,7 @@ void EwaldDisp::init_self_peratom()
|
|||
double *ei = energy+3;
|
||||
double *vi = virial+3;
|
||||
double *imu = atom->mu[0], *nmu = imu+4*atom->nlocal;
|
||||
double ce = mumurd2e*2.0*g3/3.0/sqrt(MY_PI);
|
||||
double ce = mumurd2e*2.0*g3/3.0/MY_PIS;
|
||||
for (; imu < nmu; imu += 4, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) {
|
||||
*vi = 0; // in surface
|
||||
*ei = ce*imu[3]*imu[3]-vi[0];
|
||||
|
|
@ -729,9 +731,9 @@ void EwaldDisp::compute_force()
|
|||
double *mu = atom->mu ? atom->mu[0] : NULL;
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
double *ke, c[EWALD_NFUNCS] = {
|
||||
8.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(12.0*volume),
|
||||
2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 8.0*MY_PI*mumurd2e/volume};
|
||||
double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(MY_PI)/c[3];
|
||||
8.0*MY_PI*qscale/volume, 2.0*MY_PI*MY_PIS/(12.0*volume),
|
||||
2.0*MY_PI*MY_PIS/(192.0*volume), 8.0*MY_PI*mumurd2e/volume};
|
||||
double kt = 4.0*cube(g_ewald)/3.0/MY_PIS/c[3];
|
||||
int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type;
|
||||
int func[EWALD_NFUNCS];
|
||||
|
||||
|
|
@ -848,8 +850,8 @@ void EwaldDisp::compute_energy()
|
|||
double *ke = kenergy;
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
double c[EWALD_NFUNCS] = {
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume),
|
||||
2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*MY_PIS/(24.0*volume),
|
||||
2.0*MY_PI*MY_PIS/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
double sum[EWALD_NFUNCS];
|
||||
int func[EWALD_NFUNCS];
|
||||
|
||||
|
|
@ -892,8 +894,8 @@ void EwaldDisp::compute_energy_peratom()
|
|||
const double qscale = force->qqrd2e * scale;
|
||||
double *ke = kenergy;
|
||||
double c[EWALD_NFUNCS] = {
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume),
|
||||
2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*MY_PIS/(24.0*volume),
|
||||
2.0*MY_PI*MY_PIS/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type;
|
||||
int func[EWALD_NFUNCS];
|
||||
|
||||
|
|
@ -968,8 +970,8 @@ void EwaldDisp::compute_virial()
|
|||
double *kv = kvirial;
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
double c[EWALD_NFUNCS] = {
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume),
|
||||
2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*MY_PIS/(24.0*volume),
|
||||
2.0*MY_PI*MY_PIS/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
shape sum[EWALD_NFUNCS];
|
||||
int func[EWALD_NFUNCS];
|
||||
|
||||
|
|
@ -1025,8 +1027,8 @@ void EwaldDisp::compute_virial_peratom()
|
|||
double *mu = atom->mu ? atom->mu[0] : NULL;
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
double c[EWALD_NFUNCS] = {
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume),
|
||||
2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
4.0*MY_PI*qscale/volume, 2.0*MY_PI*MY_PIS/(24.0*volume),
|
||||
2.0*MY_PI*MY_PIS/(192.0*volume), 4.0*MY_PI*mumurd2e/volume};
|
||||
shape sum[EWALD_MAX_NSUMS];
|
||||
int func[EWALD_NFUNCS];
|
||||
|
||||
|
|
@ -1200,8 +1202,8 @@ double EwaldDisp::NewtonSolve(double x, double Rc,
|
|||
double EwaldDisp::f(double x, double Rc, bigint natoms, double vol, double b2)
|
||||
{
|
||||
double a = Rc*x;
|
||||
double f = (4.0*MY_PI*b2*pow(x,4.0)/vol/sqrt((double)natoms)*erfc(a) *
|
||||
(6.0*pow(a,-5.0) + 6.0*pow(a,-3.0) + 3.0/a + a) - accuracy);
|
||||
double f = (4.0*MY_PI*b2*powint(x,4)/vol/sqrt((double)natoms)*erfc(a) *
|
||||
(6.0*powint(a,-5) + 6.0*powint(a,-3) + 3.0/a + a) - accuracy);
|
||||
return f;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,508 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "commgrid.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "msm_cg.h"
|
||||
|
||||
#include "math_const.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define OFFSET 16384
|
||||
#define SMALLQ 0.00001
|
||||
|
||||
enum{REVERSE_RHO,REVERSE_AD,REVERSE_AD_PERATOM};
|
||||
enum{FORWARD_RHO,FORWARD_AD,FORWARD_AD_PERATOM};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MSMCG::MSMCG(LAMMPS *lmp, int narg, char **arg) : MSM(lmp, narg, arg)
|
||||
{
|
||||
if ((narg < 1) || (narg > 2))
|
||||
error->all(FLERR,"Illegal kspace_style msm/cg command");
|
||||
|
||||
if (narg == 2)
|
||||
smallq = atof(arg[1]);
|
||||
else
|
||||
smallq = SMALLQ;
|
||||
|
||||
num_charged = -1;
|
||||
is_charged = NULL;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
free all memory
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
MSMCG::~MSMCG()
|
||||
{
|
||||
memory->destroy(is_charged);
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute the MSM long-range force, energy, virial
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MSMCG::compute(int eflag, int vflag)
|
||||
{
|
||||
const double * const q = atom->q;
|
||||
const int nlocal = atom->nlocal;
|
||||
int i,j,n;
|
||||
|
||||
// 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 (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 (nlocal > nmax) {
|
||||
memory->destroy(part2grid);
|
||||
memory->destroy(is_charged);
|
||||
nmax = atom->nmax;
|
||||
memory->create(part2grid,nmax,3,"msm:part2grid");
|
||||
memory->create(is_charged,nmax,"msm/cg:is_charged");
|
||||
}
|
||||
|
||||
// one time setup message
|
||||
|
||||
if (num_charged < 0) {
|
||||
bigint charged_all, charged_num;
|
||||
double charged_frac, charged_fmax, charged_fmin;
|
||||
|
||||
num_charged=0;
|
||||
for (i=0; i < nlocal; ++i)
|
||||
if (fabs(q[i]) > smallq)
|
||||
++num_charged;
|
||||
|
||||
// get fraction of charged particles per domain
|
||||
|
||||
if (nlocal > 0)
|
||||
charged_frac = static_cast<double>(num_charged) * 100.0
|
||||
/ static_cast<double>(nlocal);
|
||||
else
|
||||
charged_frac = 0.0;
|
||||
|
||||
MPI_Reduce(&charged_frac,&charged_fmax,1,MPI_DOUBLE,MPI_MAX,0,world);
|
||||
MPI_Reduce(&charged_frac,&charged_fmin,1,MPI_DOUBLE,MPI_MIN,0,world);
|
||||
|
||||
// get fraction of charged particles overall
|
||||
|
||||
charged_num = num_charged;
|
||||
MPI_Reduce(&charged_num,&charged_all,1,MPI_LMP_BIGINT,MPI_SUM,0,world);
|
||||
charged_frac = static_cast<double>(charged_all) * 100.0
|
||||
/ static_cast<double>(atom->natoms);
|
||||
|
||||
if (me == 0) {
|
||||
if (screen)
|
||||
fprintf(screen,
|
||||
" MSM/cg optimization cutoff: %g\n"
|
||||
" Total charged atoms: %.1f%%\n"
|
||||
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
|
||||
smallq,charged_frac,charged_fmin,charged_fmax);
|
||||
if (logfile)
|
||||
fprintf(logfile,
|
||||
" MSM/cg optimization cutoff: %g\n"
|
||||
" Total charged atoms: %.1f%%\n"
|
||||
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
|
||||
smallq,charged_frac,charged_fmin,charged_fmax);
|
||||
}
|
||||
}
|
||||
|
||||
num_charged = 0;
|
||||
for (i = 0; i < nlocal; ++i)
|
||||
if (fabs(q[i]) > smallq) {
|
||||
is_charged[num_charged] = i;
|
||||
++num_charged;
|
||||
}
|
||||
|
||||
// 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 (n=0; n<=levels-2; n++) {
|
||||
current_level = n;
|
||||
cg[n]->forward_comm(this,FORWARD_RHO);
|
||||
|
||||
direct(n);
|
||||
restriction(n);
|
||||
}
|
||||
|
||||
// top grid level
|
||||
|
||||
current_level = levels-1;
|
||||
cg[levels-1]->forward_comm(this,FORWARD_RHO);
|
||||
direct_top(levels-1);
|
||||
|
||||
for (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) {
|
||||
const double qs = 0.5*qscale;
|
||||
|
||||
if (eflag_atom) {
|
||||
const double sf = gamma(0.0)/cutoff;
|
||||
for (j = 0; j < num_charged; j++) {
|
||||
i = is_charged[j];
|
||||
eatom[i] -= q[i]*q[i]*sf;
|
||||
eatom[i] *= qs;
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_atom) {
|
||||
for (n = 0; n < num_charged; n++) {
|
||||
i = is_charged[n];
|
||||
for (j = 0; j < 6; j++)
|
||||
vatom[i][j] *= qs;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
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 MSMCG::particle_map()
|
||||
{
|
||||
const double * const * const x = atom->x;
|
||||
|
||||
int flag = 0;
|
||||
int i;
|
||||
|
||||
for (int j = 0; j < num_charged; j++) {
|
||||
i = is_charged[j];
|
||||
|
||||
// (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
|
||||
|
||||
const int nx=static_cast<int>((x[i][0]-boxlo[0])*delxinv[0]+OFFSET)-OFFSET;
|
||||
const int ny=static_cast<int>((x[i][1]-boxlo[1])*delyinv[0]+OFFSET)-OFFSET;
|
||||
const int 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 MSMCG::make_rho()
|
||||
{
|
||||
const double * const q = atom->q;
|
||||
const double * const * const x = atom->x;
|
||||
|
||||
// clear 3d density array
|
||||
|
||||
double * const * const * const qgridn = qgrid[0];
|
||||
|
||||
memset(&(qgridn[nzlo_out[0]][nylo_out[0]][nxlo_out[0]]),0,ngrid[0]*sizeof(double));
|
||||
|
||||
double dx,dy,dz,x0,y0,z0;
|
||||
int i,j,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
|
||||
// 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
|
||||
|
||||
for (j = 0; j < num_charged; j++) {
|
||||
i = is_charged[j];
|
||||
|
||||
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];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
interpolate from grid to get force on my particles
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MSMCG::fieldforce()
|
||||
{
|
||||
|
||||
const double * const * const * const egridn = egrid[0];
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = atom->f;
|
||||
const double * const q = atom->q;
|
||||
|
||||
int i,j,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
|
||||
|
||||
for (j = 0; j < num_charged; j++) {
|
||||
i = is_charged[j];
|
||||
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 MSMCG::fieldforce_peratom()
|
||||
{
|
||||
const double * const q = atom->q;
|
||||
const double * const * const x = atom->x;
|
||||
|
||||
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];
|
||||
|
||||
int i,j,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
double dx,dy,dz,x0,y0,z0;
|
||||
double u,v0,v1,v2,v3,v4,v5;
|
||||
|
||||
// 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
|
||||
|
||||
for (j = 0; j < num_charged; j++) {
|
||||
i = is_charged[j];
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double MSMCG::memory_usage()
|
||||
{
|
||||
double bytes = MSM::memory_usage();
|
||||
bytes += nmax * sizeof(int);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,141 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef KSPACE_CLASS
|
||||
|
||||
KSpaceStyle(msm/cg,MSMCG)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_MSM_CG_H
|
||||
#define LMP_MSM_CG_H
|
||||
|
||||
#include "msm.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class MSMCG : public MSM {
|
||||
public:
|
||||
MSMCG(class LAMMPS *, int, char **);
|
||||
virtual ~MSMCG();
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
protected:
|
||||
int num_charged;
|
||||
int *is_charged;
|
||||
double smallq;
|
||||
|
||||
protected:
|
||||
virtual void particle_map();
|
||||
virtual void make_rho();
|
||||
virtual void fieldforce();
|
||||
virtual void fieldforce_peratom();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Cannot (yet) use MSM with triclinic box
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot (yet) use MSM with 2d simulation
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Kspace style requires atom attribute q
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
E: Cannot use slab correction with MSM
|
||||
|
||||
Slab correction can only be used with Ewald and PPPM, not MSM.
|
||||
|
||||
E: MSM order must be 4, 6, 8, or 10
|
||||
|
||||
This is a limitation of the MSM implementation in LAMMPS:
|
||||
the MSM order can only be 4, 6, 8, or 10.
|
||||
|
||||
E: Cannot (yet) use single precision with MSM (remove -DFFT_SINGLE from Makefile and recompile)
|
||||
|
||||
Single precision cannot be used with MSM.
|
||||
|
||||
E: KSpace style is incompatible with Pair style
|
||||
|
||||
Setting a kspace style requires that a pair style with a long-range
|
||||
Coulombic component be selected that is compatible with MSM. Note
|
||||
that TIP4P is not (yet) supported by MSM.
|
||||
|
||||
E: Cannot use kspace solver on system with no charge
|
||||
|
||||
No atoms in system have a non-zero charge.
|
||||
|
||||
E: System is not charge neutral, net charge = %g
|
||||
|
||||
The total charge on all atoms on the system is not 0.0, which
|
||||
is not valid for MSM.
|
||||
|
||||
E: MSM grid is too large
|
||||
|
||||
The global MSM grid is larger than OFFSET in one or more dimensions.
|
||||
OFFSET is currently set to 16384. You likely need to decrease the
|
||||
requested accuracy.
|
||||
|
||||
W: MSM mesh too small, increasing to 2 points in each direction
|
||||
|
||||
The global MSM grid is too small, so the number of grid points has been
|
||||
increased
|
||||
|
||||
E: KSpace accuracy must be > 0
|
||||
|
||||
The kspace accuracy designated in the input must be greater than zero.
|
||||
|
||||
W: Number of MSM mesh points increased to be a multiple of 2
|
||||
|
||||
MSM requires that the number of grid points in each direction be a multiple
|
||||
of two and the number of grid points in one or more directions have been
|
||||
adjusted to meet this requirement.
|
||||
|
||||
W: Adjusting Coulombic cutoff for MSM, new cutoff = %g
|
||||
|
||||
The adjust/cutoff command is turned on and the Coulombic cutoff has been
|
||||
adjusted to match the user-specified accuracy.
|
||||
|
||||
E: Out of range atoms - cannot compute MSM
|
||||
|
||||
One or more atoms are attempting to map their charge to a MSM grid point
|
||||
that is not owned by a processor. This is likely for one of two
|
||||
reasons, both of them bad. First, it may mean that an atom near the
|
||||
boundary of a processor's sub-domain has moved more than 1/2 the
|
||||
"neighbor skin distance"_neighbor.html without neighbor lists being
|
||||
rebuilt and atoms being migrated to new processors. This also means
|
||||
you may be missing pairwise interactions that need to be computed.
|
||||
The solution is to change the re-neighboring criteria via the
|
||||
"neigh_modify"_neigh_modify command. The safest settings are "delay 0
|
||||
every 1 check yes". Second, it may mean that an atom has moved far
|
||||
outside a processor's sub-domain or even the entire simulation box.
|
||||
This indicates bad physics, e.g. due to highly overlapping atoms, too
|
||||
large a timestep, etc.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue