ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

modernize and reformat code

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-06-05 21:43:38 -04:00
parent f7ca10b070
commit 0bc86a7eea
3 changed files with 461 additions and 476 deletions
Download Patch File Download Diff File Expand all files Collapse all files

378
src/USER-DIELECTRIC/fix_polarize_bem_gmres.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
@ -47,13 +46,13 @@
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "msm_dielectric.h"
#include "pair_coul_cut_dielectric.h"
#include "pair_coul_long_dielectric.h"
#include "pair_lj_cut_coul_cut_dielectric.h"
#include "pair_lj_cut_coul_long_dielectric.h"
#include "pair_lj_cut_coul_msm_dielectric.h"
#include "pppm_dielectric.h"
#include "msm_dielectric.h"
#include "random_park.h"
#include "timer.h"
#include "update.h"
@ -70,19 +69,19 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
FixPolarizeBEMGMRES::FixPolarizeBEMGMRES(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
q_backup(NULL), c(NULL), g(NULL), h(NULL), r(NULL), s(NULL), v(NULL), y(NULL)
Fix(lmp, narg, arg), q_backup(NULL), c(NULL), g(NULL), h(NULL), r(NULL), s(NULL), v(NULL),
y(NULL)
{
if (narg < 5) error->all(FLERR,"Illegal fix polarize/bem/gmres command");
if (narg < 5) error->all(FLERR, "Illegal fix polarize/bem/gmres command");
avec = (AtomVecDielectric *) atom->style_match("dielectric");
if (!avec) error->all(FLERR,"Fix polarize requires atom style dielectric");
if (!avec) error->all(FLERR, "Fix polarize requires atom style dielectric");
// parse required arguments
nevery = utils::numeric(FLERR,arg[3],false,lmp);
if (nevery < 0) error->all(FLERR,"Illegal fix polarize/bem/gmres command");
double tol = utils::numeric(FLERR,arg[4],false,lmp);
nevery = utils::numeric(FLERR, arg[3], false, lmp);
if (nevery < 0) error->all(FLERR, "Illegal fix polarize/bem/gmres command");
double tol = utils::numeric(FLERR, arg[4], false, lmp);
tol_abs = tol_rel = tol;
itr_max = 20;
@ -112,7 +111,7 @@ FixPolarizeBEMGMRES::FixPolarizeBEMGMRES(LAMMPS *lmp, int narg, char **arg) :
if (atom->avec->forceclearflag) extraflag = 1;
grow_arrays(atom->nmax);
atom->add_callback(0); // to ensure to work with atom->sort()
atom->add_callback(0); // to ensure to work with atom->sort()
// output the residual and actual number of iterations
@ -135,7 +134,7 @@ FixPolarizeBEMGMRES::~FixPolarizeBEMGMRES()
memory->destroy(tag2mat);
if (allocated) deallocate();
atom->delete_callback(id,0);
atom->delete_callback(id, 0);
}
/* ---------------------------------------------------------------------- */
@ -153,7 +152,7 @@ void FixPolarizeBEMGMRES::init()
{
// mapping induced charge matrix/vector to atom tags and vice versa
int i,maxtag;
int i, maxtag;
double *q = atom->q;
int *mask = atom->mask;
tagint *tag = atom->tag;
@ -161,28 +160,30 @@ void FixPolarizeBEMGMRES::init()
tagint max_tag = -1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) max_tag = MAX(max_tag,tag[i]);
if (mask[i] & groupbit) max_tag = MAX(max_tag, tag[i]);
tagint itmp;
MPI_Allreduce(&max_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
MPI_Allreduce(&max_tag, &itmp, 1, MPI_LMP_TAGINT, MPI_MAX, world);
maxtag = (int) itmp;
int *ncount;
memory->create(ncount,maxtag+1,"polarize:ncount");
memory->create(ncount, maxtag + 1, "polarize:ncount");
for (i = 0; i <= maxtag; i++) ncount[i] = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) ncount[tag[i]]++;
memory->create(tag2mat,maxtag+1,"polarize:tag2mat");
MPI_Allreduce(ncount,tag2mat,maxtag+1,MPI_INT,MPI_SUM,world);
memory->create(tag2mat, maxtag + 1, "polarize:tag2mat");
MPI_Allreduce(ncount, tag2mat, maxtag + 1, MPI_INT, MPI_SUM, world);
num_induced_charges = 0;
for (i = 0; i <= maxtag; i++)
if (tag2mat[i]) tag2mat[i] = num_induced_charges++;
else tag2mat[i] = -1;
if (tag2mat[i])
tag2mat[i] = num_induced_charges++;
else
tag2mat[i] = -1;
memory->create(mat2tag,num_induced_charges,"polarize:mat2tag");
memory->create(mat2tag, num_induced_charges, "polarize:mat2tag");
num_induced_charges = 0;
for (i = 0; i <= maxtag; i++)
@ -197,9 +198,9 @@ void FixPolarizeBEMGMRES::init()
// allocate memory for the solver
memory->create(induced_charges,num_induced_charges,"polarize:induced_charges");
memory->create(rhs,num_induced_charges,"polarize:rhs");
memory->create(buffer,num_induced_charges,"polarize:buffer");
memory->create(induced_charges, num_induced_charges, "polarize:induced_charges");
memory->create(rhs, num_induced_charges, "polarize:rhs");
memory->create(buffer, num_induced_charges, "polarize:buffer");
mat_dim = num_induced_charges;
if (mr > mat_dim - 1 || mr <= 0) mr = mat_dim - 1;
@ -213,16 +214,16 @@ void FixPolarizeBEMGMRES::init()
if (randomized) {
RanPark *random = new RanPark(lmp,seed_charge + comm->me);
RanPark *random = new RanPark(lmp, seed_charge + comm->me);
for (i = 0; i < 100; i++) random->uniform();
double sum,tmp = 0;
double sum, tmp = 0;
for (i = 0; i < nlocal; i++) {
if (induced_charge_idx[i] < 0) continue;
q[i] = ave_charge*(random->uniform() - 0.5);
q[i] = ave_charge * (random->uniform() - 0.5);
tmp += q[i];
}
MPI_Allreduce(&tmp,&sum,1,MPI_DOUBLE,MPI_SUM,world);
sum /= (double)num_induced_charges;
MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
sum /= (double) num_induced_charges;
tmp = 0;
for (i = 0; i < nlocal; i++) {
@ -230,20 +231,17 @@ void FixPolarizeBEMGMRES::init()
q[i] -= sum;
tmp += q[i];
}
MPI_Allreduce(&tmp,&sum,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
if (comm->me == 0) {
if (screen) fprintf(screen, "ave induced charge q = %g\n", sum);
}
if (comm->me == 0) utils::logmesg(lmp, "ave induced charge q = {:.8}\n", sum);
delete random;
}
if (comm->me == 0) {
if (screen) fprintf(screen,"GMRES solver for %d induced charges "
"using maximum %d q-vectors\n",num_induced_charges,mr);
if (logfile) fprintf(logfile,"GMRES solver for %d induced charges "
"using maximum %d q-vectors\n",num_induced_charges,mr);
}
if (comm->me == 0)
utils::logmesg(lmp,
"GMRES solver for {} induced charges "
"using maximum {} q-vectors\n",
num_induced_charges, mr);
}
/* ---------------------------------------------------------------------- */
@ -252,30 +250,33 @@ void FixPolarizeBEMGMRES::setup(int vflag)
{
// check if the pair styles in use are compatible
if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric") == 0)
efield_pair = ((PairLJCutCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/msm/dielectric") == 0)
efield_pair = ((PairLJCutCoulMSMDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric") == 0)
efield_pair = ((PairLJCutCoulCutDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulCutDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"coul/long/dielectric") == 0)
efield_pair = ((PairCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"coul/cut/dielectric") == 0)
efield_pair = ((PairCoulCutDielectric*)force->pair)->efield;
else error->all(FLERR,"Pair style not compatible with fix polarize");
if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0)
efield_pair = ((PairLJCutCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = ((PairLJCutCoulMSMDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = ((PairLJCutCoulCutDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulCutDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = ((PairCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = ((PairCoulCutDielectric *) force->pair)->efield;
else
error->all(FLERR, "Pair style not compatible with fix polarize");
if (kspaceflag) {
if (force->kspace) {
if (strcmp(force->kspace_style,"pppm/dielectric") == 0)
efield_kspace = ((PPPMDielectric*)force->kspace)->efield;
else if (strcmp(force->kspace_style,"msm/dielectric") == 0)
efield_kspace = ((MSMDielectric*)force->kspace)->efield;
else error->all(FLERR,"Kspace style not compatible with fix polarize/bem/gmres");
} else error->all(FLERR,"No Kspace style available for fix polarize/bem/gmres");
if (strcmp(force->kspace_style, "pppm/dielectric") == 0)
efield_kspace = ((PPPMDielectric *) force->kspace)->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = ((MSMDielectric *) force->kspace)->efield;
else
error->all(FLERR, "Kspace style not compatible with fix polarize/bem/gmres");
} else
error->all(FLERR, "No Kspace style available for fix polarize/bem/gmres");
}
first = 1;
@ -336,8 +337,8 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
// Note: the right-hand side (b) is in the unit of charge density
force_clear();
force->pair->compute(eflag,vflag);
if (kspaceflag) force->kspace->compute(eflag,vflag);
force->pair->compute(eflag, vflag);
if (kspaceflag) force->kspace->compute(eflag, vflag);
if (force->newton) comm->reverse_comm();
for (int i = 0; i < num_induced_charges; i++) buffer[i] = 0;
@ -357,12 +358,12 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
Ey += efield_kspace[i][1];
Ez += efield_kspace[i][2];
}
double dot = (Ex*norm[i][0] + Ey*norm[i][1] + Ez*norm[i][2]) / epsilon[i];
double dot = (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / epsilon[i];
double sigma_f = q_real[i] / area[i];
buffer[idx] = (1 - em[i]) * sigma_f - epsilon0 * ed[i] * dot / (4*MY_PI);
buffer[idx] = (1 - em[i]) * sigma_f - epsilon0 * ed[i] * dot / (4 * MY_PI);
}
MPI_Allreduce(buffer,rhs,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer, rhs, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
// compute the initial residual r before iteration
// while it seems that assigning induced charges to the last values
@ -391,7 +392,7 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
for (int i = 0; i < nlocal; i++) {
if (induced_charge_idx[i] >= 0) {
int idx = induced_charge_idx[i];
q[i] = induced_charges[idx]*area[i] + q_real[i];
q[i] = induced_charges[idx] * area[i] + q_real[i];
} else {
q[i] = q_backup[i];
}
@ -400,22 +401,21 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
comm->forward_comm_fix(this);
if (first) first = 0;
}
/* ---------------------------------------------------------------------- */
void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
void FixPolarizeBEMGMRES::gmres_solve(double *x, double *r)
{
int i,j,k,k_copy,n,itr;
double av,htmp,mu,rho_tol;
int i, j, k, k_copy, n, itr;
double av, htmp, mu, rho_tol;
double delta = 1.0e-03;
n = mat_dim;
// compute the relative tolerance
// rho = norm(r)
rho = sqrt( vec_dot(r, r, n) );
rho = sqrt(vec_dot(r, r, n));
rho_tol = rho * tol_rel;
// the outer loop to itr_max
@ -425,15 +425,14 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
// the first vector v (i.e. v[0]) is the updated residual normalized
for (i = 0; i < n; i++)
v[i+0*n] = r[i] / rho;
for (i = 0; i < n; i++) v[i + 0 * n] = r[i] / rho;
g[0] = rho;
for (i = 1; i <= mr; i++) g[i] = 0.0;
// fill up h with zero
memset(h, 0, (mr+1)*mr*sizeof(double));
memset(h, 0, (mr + 1) * mr * sizeof(double));
// the inner loop k = 1..(n-1)
// build up the k-th Krylov space,
@ -448,43 +447,40 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
// here is the tricky part: v(k-1) plays a role as "charges"
// matvec(a, v+(k-1)*n, v+k*n, n);
apply_operator(v+(k-1)*n, v+k*n, n);
apply_operator(v + (k - 1) * n, v + k * n, n);
// compute the norm of the vector v(k)
av = sqrt(vec_dot(v+k*n, v+k*n, n));
av = sqrt(vec_dot(v + k * n, v + k * n, n));
// Arnoldi iteration to find v's
// orthogonalize the k vectors v(1) . . . v(k)
for (j = 1; j <= k; j++) {
h[(j-1)+(k-1)*(mr+1)] = vec_dot(v+k*n, v+(j-1)*n, n);
h[(j - 1) + (k - 1) * (mr + 1)] = vec_dot(v + k * n, v + (j - 1) * n, n);
for (i = 0; i < n; i++)
v[i+k*n] = v[i+k*n] - h[(j-1)+(k-1)*(mr+1)] * v[i+(j-1)*n];
v[i + k * n] = v[i + k * n] - h[(j - 1) + (k - 1) * (mr + 1)] * v[i + (j - 1) * n];
}
// compute the norm of the newly created vector v(k)
h[k+(k-1)*(mr+1)] = sqrt(vec_dot(v+k*n, v+k*n, n));
h[k + (k - 1) * (mr + 1)] = sqrt(vec_dot(v + k * n, v + k * n, n));
// if the norm is close to zero, repeat the above orthogonalization
if ((av + delta * h[k+(k-1)*(mr+1)]) == av) {
if ((av + delta * h[k + (k - 1) * (mr + 1)]) == av) {
for (j = 1; j <= k; j++) {
htmp = vec_dot(v+k*n, v+(j-1)*n, n);
h[(j-1)+(k-1)*(mr+1)] = h[(j-1)+(k-1)*(mr+1)] + htmp;
for (i = 0; i < n; i++)
v[i+k*n] = v[i+k*n] - htmp * v[i+(j-1)*n];
htmp = vec_dot(v + k * n, v + (j - 1) * n, n);
h[(j - 1) + (k - 1) * (mr + 1)] = h[(j - 1) + (k - 1) * (mr + 1)] + htmp;
for (i = 0; i < n; i++) v[i + k * n] = v[i + k * n] - htmp * v[i + (j - 1) * n];
}
h[k+(k-1)*(mr+1)] = sqrt( vec_dot(v+k*n, v+k*n, n) );
h[k + (k - 1) * (mr + 1)] = sqrt(vec_dot(v + k * n, v + k * n, n));
}
// if the norm of v(k) is nonzero, normalize v(k)
if (h[k+(k-1)*(mr+1)] != 0.0) {
for (i = 0; i < n; i++) {
v[i+k*n] = v[i+k*n] / h[k+(k-1)*(mr+1)];
}
if (h[k + (k - 1) * (mr + 1)] != 0.0) {
for (i = 0; i < n; i++) { v[i + k * n] = v[i + k * n] / h[k + (k - 1) * (mr + 1)]; }
}
// if k is not the first iteration,
@ -495,71 +491,65 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
// update y(i-1) <- h(k-1, i-1) for i = 1...(k+1)
for (i = 1; i <= k + 1; i++)
y[i-1] = h[(i-1)+(k-1)*(mr+1)];
for (i = 1; i <= k + 1; i++) y[i - 1] = h[(i - 1) + (k - 1) * (mr + 1)];
// apply the Given rotation to y[j-1] and y[j] for j = 1..(k-1)
for (j = 1; j <= k - 1; j++)
mult_givens(c[j-1], s[j-1], j-1, y);
for (j = 1; j <= k - 1; j++) mult_givens(c[j - 1], s[j - 1], j - 1, y);
// update h(k-1, i-1) <- y(i-1) for i = 1..(k_1)
for (i = 1; i <= k + 1; i++)
h[i-1+(k-1)*(mr+1)] = y[i-1];
for (i = 1; i <= k + 1; i++) h[i - 1 + (k - 1) * (mr + 1)] = y[i - 1];
}
// compute cosine and sine terms of the Given rotations
mu = sqrt(h[(k-1)+(k-1)*(mr+1)]*h[(k-1)+(k-1)*(mr+1)]
+ h[ k +(k-1)*(mr+1)]*h[ k +(k-1)*(mr+1)]);
c[k-1] = h[(k-1)+(k-1)*(mr+1)] / mu;
s[k-1] = -h[ k +(k-1)*(mr+1)] / mu;
mu = sqrt(h[(k - 1) + (k - 1) * (mr + 1)] * h[(k - 1) + (k - 1) * (mr + 1)] +
h[k + (k - 1) * (mr + 1)] * h[k + (k - 1) * (mr + 1)]);
c[k - 1] = h[(k - 1) + (k - 1) * (mr + 1)] / mu;
s[k - 1] = -h[k + (k - 1) * (mr + 1)] / mu;
// update h(k-1,k-1) and set h(k-1,k) to zero
h[(k-1)+(k-1)*(mr+1)] = c[k-1] * h[(k-1)+(k-1)*(mr+1)]
- s[k-1] * h[ k +(k-1)*(mr+1)];
h[k +(k-1)*(mr+1)] = 0;
h[(k - 1) + (k - 1) * (mr + 1)] =
c[k - 1] * h[(k - 1) + (k - 1) * (mr + 1)] - s[k - 1] * h[k + (k - 1) * (mr + 1)];
h[k + (k - 1) * (mr + 1)] = 0;
// apply the Givens rotation to g[k-1] and g[k]
mult_givens(c[k-1], s[k-1], k-1, g);
mult_givens(c[k - 1], s[k - 1], k - 1, g);
// compute the norm of the residual
rho = fabs(g[k]);
#ifdef _POLARIZE_DEBUG
#ifdef _POLARIZE_DEBUG
if (comm->me == 0) {
char message[256];
sprintf(message, "itr = %d: k = %d, norm(r) = %g norm(b) = %g",
itr, k, rho, normb);
sprintf(message, "itr = %d: k = %d, norm(r) = %g norm(b) = %g", itr, k, rho, normb);
error->warning(FLERR, message);
}
#endif
#endif
if (rho <= rho_tol && rho <= tol_abs)
break;
if (rho <= rho_tol && rho <= tol_abs) break;
}
k = k_copy - 1;
// compute the estimate y from h
y[k] = g[k] / h[k + k*(mr+1)];
y[k] = g[k] / h[k + k * (mr + 1)];
for (i = k; i >= 1; i--) {
y[i-1] = g[i-1];
y[i - 1] = g[i - 1];
for (j = i + 1; j <= k + 1; j++)
y[i-1] = y[i-1] - h[(i-1)+(j-1)*(mr+1)] * y[j-1];
y[i-1] = y[i-1] / h[(i-1)+(i-1)*(mr+1)];
y[i - 1] = y[i - 1] - h[(i - 1) + (j - 1) * (mr + 1)] * y[j - 1];
y[i - 1] = y[i - 1] / h[(i - 1) + (i - 1) * (mr + 1)];
}
// update x at the current iteration: x <- Q(n by k) * y (k by 1)
for (i = 1; i <= n; i++) {
for (j = 1; j <= k + 1; j++)
x[i-1] = x[i-1] + v[(i-1)+(j-1)*n] * y[j-1];
for (j = 1; j <= k + 1; j++) x[i - 1] = x[i - 1] + v[(i - 1) + (j - 1) * n] * y[j - 1];
}
// update the residual with the updated induced charges (x)
@ -568,16 +558,15 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
// rho = norm(r)
rho = sqrt( vec_dot(r, r, n) );
rho = sqrt(vec_dot(r, r, n));
#ifdef _POLARIZE_DEBUG
#ifdef _POLARIZE_DEBUG
if (comm->me == 0) {
char message[256];
sprintf(message, "itr = %d: norm(r) = %g norm(b) = %g",
itr, rho, normb);
sprintf(message, "itr = %d: norm(r) = %g norm(b) = %g", itr, rho, normb);
error->warning(FLERR, message);
}
#endif
#endif
// Barros et al. suggested the condition: norm(r) < EPSILON norm(b)
@ -596,7 +585,7 @@ void FixPolarizeBEMGMRES::gmres_solve(double* x, double* r)
matvec(A, v(k-1), v(k), n);
------------------------------------------------------------------------- */
void FixPolarizeBEMGMRES::apply_operator(double* w, double* Aw, int n)
void FixPolarizeBEMGMRES::apply_operator(double *w, double *Aw, int n)
{
int i;
double *q = atom->q;
@ -621,7 +610,7 @@ void FixPolarizeBEMGMRES::apply_operator(double* w, double* Aw, int n)
q[i] = 0;
} else {
int idx = induced_charge_idx[i];
q[i] = w[idx]*area[i];
q[i] = w[idx] * area[i];
}
}
@ -630,8 +619,8 @@ void FixPolarizeBEMGMRES::apply_operator(double* w, double* Aw, int n)
// compute the electrical field due to w*area: y = A (w*area)
force_clear();
force->pair->compute(eflag,vflag);
if (kspaceflag) force->kspace->compute(eflag,vflag);
force->pair->compute(eflag, vflag);
if (kspaceflag) force->kspace->compute(eflag, vflag);
if (force->newton) comm->reverse_comm();
// now efield is the electrical field due to induced charges only
@ -652,11 +641,11 @@ void FixPolarizeBEMGMRES::apply_operator(double* w, double* Aw, int n)
Ey += efield_kspace[i][1];
Ez += efield_kspace[i][2];
}
double dot = (Ex*norm[i][0] + Ey*norm[i][1] + Ez*norm[i][2]) / epsilon[i];
buffer[idx] = em[i] * w[idx] + epsilon0 * ed[i] * dot / (4*MY_PI);
double dot = (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / epsilon[i];
buffer[idx] = em[i] * w[idx] + epsilon0 * ed[i] * dot / (4 * MY_PI);
}
MPI_Allreduce(buffer,Aw,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer, Aw, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
}
/* ----------------------------------------------------------------------
@ -666,7 +655,7 @@ void FixPolarizeBEMGMRES::apply_operator(double* w, double* Aw, int n)
using Eq. (60) in Barros et al.
------------------------------------------------------------------------ */
void FixPolarizeBEMGMRES::update_residual(double* w, double* r, int n)
void FixPolarizeBEMGMRES::update_residual(double *w, double *r, int n)
{
int i;
double *q = atom->q;
@ -692,15 +681,15 @@ void FixPolarizeBEMGMRES::update_residual(double* w, double* r, int n)
q[i] = q_backup[i];
} else {
int idx = induced_charge_idx[i];
q[i] = w[idx]*area[i] + q_real[i];
q[i] = w[idx] * area[i] + q_real[i];
}
}
comm->forward_comm_fix(this);
force_clear();
force->pair->compute(eflag,vflag);
if (kspaceflag) force->kspace->compute(eflag,vflag);
force->pair->compute(eflag, vflag);
if (kspaceflag) force->kspace->compute(eflag, vflag);
if (force->newton) comm->reverse_comm();
// compute the residual according to Eq. (60) in Barros et al.
@ -725,13 +714,12 @@ void FixPolarizeBEMGMRES::update_residual(double* w, double* r, int n)
Ey += efield_kspace[i][1];
Ez += efield_kspace[i][2];
}
double dot = (Ex*norm[i][0] + Ey*norm[i][1] + Ez*norm[i][2]) / epsilon[i];
double dot = (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / epsilon[i];
double sigma_f = q_real[i] / area[i];
buffer[idx] = (1 - em[i]) * sigma_f - em[i] * w[idx] -
epsilon0 * ed[i] * dot / (4*MY_PI);
buffer[idx] = (1 - em[i]) * sigma_f - em[i] * w[idx] - epsilon0 * ed[i] * dot / (4 * MY_PI);
}
MPI_Allreduce(buffer,r,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer, r, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
}
/* ---------------------------------------------------------------------- */
@ -742,15 +730,15 @@ void FixPolarizeBEMGMRES::force_clear()
if (force->newton) nbytes += sizeof(double) * atom->nghost;
if (nbytes) {
memset(&atom->f[0][0],0,3*nbytes);
if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
if (extraflag) atom->avec->force_clear(0,nbytes);
memset(&atom->f[0][0], 0, 3 * nbytes);
if (torqueflag) memset(&atom->torque[0][0], 0, 3 * nbytes);
if (extraflag) atom->avec->force_clear(0, nbytes);
}
}
/* ---------------------------------------------------------------------- */
double FixPolarizeBEMGMRES::vec_dot(const double* a1, const double* a2, int n)
double FixPolarizeBEMGMRES::vec_dot(const double *a1, const double *a2, int n)
{
double value = 0.0;
for (int i = 0; i < n; i++) value += (a1[i] * a2[i]);
@ -764,18 +752,18 @@ double FixPolarizeBEMGMRES::vec_dot(const double* a1, const double* a2, int n)
double FixPolarizeBEMGMRES::memory_usage()
{
double bytes = 0;
bytes += mat_dim*sizeof(double); // induced_charges
bytes += mat_dim*sizeof(double); // buffer
bytes += mat_dim*sizeof(double); // rhs
bytes += atom->nmax*sizeof(double); // induced_charge_idx
bytes += atom->nmax*sizeof(double); // q_backup
bytes += mr*sizeof(double); // c
bytes += (mr+1)*sizeof(double); // g
bytes += (mr+1)*mr*sizeof(double); // h
bytes += mat_dim*sizeof(double); // r
bytes += mr*(mr+1)*sizeof(double); // s
bytes += mat_dim*sizeof(double); // v
bytes += (mr+1)*mr*sizeof(double); // y
bytes += mat_dim * sizeof(double); // induced_charges
bytes += mat_dim * sizeof(double); // buffer
bytes += mat_dim * sizeof(double); // rhs
bytes += atom->nmax * sizeof(double); // induced_charge_idx
bytes += atom->nmax * sizeof(double); // q_backup
bytes += mr * sizeof(double); // c
bytes += (mr + 1) * sizeof(double); // g
bytes += (mr + 1) * mr * sizeof(double); // h
bytes += mat_dim * sizeof(double); // r
bytes += mr * (mr + 1) * sizeof(double); // s
bytes += mat_dim * sizeof(double); // v
bytes += (mr + 1) * mr * sizeof(double); // y
return bytes;
}
@ -784,12 +772,12 @@ double FixPolarizeBEMGMRES::memory_usage()
void FixPolarizeBEMGMRES::allocate()
{
memory->create(c, mr, "polarize:c");
memory->create(g, mr+1, "polarize:g");
memory->create(h, (mr+1)*mr, "polarize:h");
memory->create(g, mr + 1, "polarize:g");
memory->create(h, (mr + 1) * mr, "polarize:h");
memory->create(r, mat_dim, "polarize:r");
memory->create(s, mr, "polarize:s");
memory->create(v, mat_dim*(mr+1), "polarize:v");
memory->create(y, mr+1, "polarize:y");
memory->create(v, mat_dim * (mr + 1), "polarize:v");
memory->create(y, mr + 1, "polarize:y");
}
/* ---------------------------------------------------------------------- */
@ -811,45 +799,49 @@ int FixPolarizeBEMGMRES::modify_param(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"itr_max") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
itr_max = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (strcmp(arg[iarg], "itr_max") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
itr_max = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"mr") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
mr = utils::numeric(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "mr") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
mr = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) kspaceflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) kspaceflag = 0;
else error->all(FLERR,"Illegal fix_modify command for fix polarize");
} else if (strcmp(arg[iarg], "kspace") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
if (strcmp(arg[iarg + 1], "yes") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg + 1], "no") == 0)
kspaceflag = 0;
else
error->all(FLERR, "Illegal fix_modify command for fix polarize");
iarg += 2;
} else if (strcmp(arg[iarg],"dielectrics") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix_modify command");
double epsiloni=-1, areai=-1;
double qreali=0;
int set_charge=0;
double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (strcmp(arg[iarg+3],"NULL") != 0)
epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (strcmp(arg[iarg+4],"NULL") != 0)
areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
if (strcmp(arg[iarg+5],"NULL") != 0) {
qreali = utils::numeric(FLERR,arg[iarg+5],false,lmp);
} else if (strcmp(arg[iarg], "dielectrics") == 0) {
if (iarg + 6 > narg) error->all(FLERR, "Illegal fix_modify command");
double epsiloni = -1, areai = -1;
double qreali = 0;
int set_charge = 0;
double ediff = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
double emean = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
if (strcmp(arg[iarg + 3], "NULL") != 0)
epsiloni = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (strcmp(arg[iarg + 4], "NULL") != 0)
areai = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
if (strcmp(arg[iarg + 5], "NULL") != 0) {
qreali = utils::numeric(FLERR, arg[iarg + 5], false, lmp);
set_charge = 1;
}
set_dielectric_params(ediff, emean, epsiloni, areai, set_charge, qreali);
iarg += 6;
} else if (strcmp(arg[iarg],"rand") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix_modify command");
ave_charge = utils::numeric(FLERR,arg[iarg+1],false,lmp);
seed_charge = utils::numeric(FLERR,arg[iarg+2],false,lmp);
} else if (strcmp(arg[iarg], "rand") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal fix_modify command");
ave_charge = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
seed_charge = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
randomized = 1;
iarg += 3;
} else error->all(FLERR,"Illegal fix_modify command");
} else
error->all(FLERR, "Illegal fix_modify command");
}
return iarg;
@ -862,8 +854,8 @@ int FixPolarizeBEMGMRES::modify_param(int narg, char **arg)
void FixPolarizeBEMGMRES::grow_arrays(int n)
{
if (n > nmax) nmax = n;
memory->grow(induced_charge_idx,nmax,"polarize:induced_charge_idx");
memory->grow(q_backup,nmax,"polarize:q_backup");
memory->grow(induced_charge_idx, nmax, "polarize:induced_charge_idx");
memory->grow(q_backup, nmax, "polarize:q_backup");
}
/* ----------------------------------------------------------------------
@ -886,8 +878,7 @@ void FixPolarizeBEMGMRES::set_arrays(int i)
/* ---------------------------------------------------------------------- */
int FixPolarizeBEMGMRES::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int FixPolarizeBEMGMRES::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int m;
for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];
@ -929,17 +920,20 @@ int FixPolarizeBEMGMRES::unpack_exchange(int nlocal, double *buf)
double FixPolarizeBEMGMRES::compute_vector(int n)
{
if (n == 0) return iterations;
else if (n == 1) return rho;
else return 0;
if (n == 0)
return iterations;
else if (n == 1)
return rho;
else
return 0;
}
/* ----------------------------------------------------------------------
set dielectric params for the atoms in the group
------------------------------------------------------------------------- */
void FixPolarizeBEMGMRES::set_dielectric_params(double ediff, double emean,
double epsiloni, double areai, int set_charge, double qvalue)
void FixPolarizeBEMGMRES::set_dielectric_params(double ediff, double emean, double epsiloni,
double areai, int set_charge, double qvalue)
{
double *area = atom->area;
double *ed = atom->ed;

200
src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
@ -43,7 +42,6 @@
#include "pair_lj_cut_coul_msm_dielectric.h"
#include "pppm_dielectric.h"
#include "random_park.h"
//#include "timer.h"
#include "update.h"
#include <cmath>
@ -57,19 +55,18 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 5) error->all(FLERR,"Illegal fix polarize/bem/icc command");
if (narg < 5) error->all(FLERR, "Illegal fix polarize/bem/icc command");
avec = (AtomVecDielectric *) atom->style_match("dielectric");
if (!avec) error->all(FLERR,"Fix polarize requires atom style dielectric");
if (!avec) error->all(FLERR, "Fix polarize requires atom style dielectric");
// parse required arguments
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery < 0) error->all(FLERR,"Illegal fix polarize/bem/icc command");
double tol = utils::numeric(FLERR,arg[4],false,lmp);
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
if (nevery < 0) error->all(FLERR, "Illegal fix polarize/bem/icc command");
double tol = utils::numeric(FLERR, arg[4], false, lmp);
tol_abs = tol_rel = tol;
itr_max = 20;
@ -97,9 +94,7 @@ FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *lmp, int narg, char **arg) :
/* ---------------------------------------------------------------------- */
FixPolarizeBEMICC::~FixPolarizeBEMICC()
{
}
FixPolarizeBEMICC::~FixPolarizeBEMICC() {}
/* ---------------------------------------------------------------------- */
@ -115,10 +110,7 @@ int FixPolarizeBEMICC::setmask()
void FixPolarizeBEMICC::init()
{
int ncount = group->count(igroup);
if (comm->me == 0) {
if (screen) fprintf(screen,"BEM/ICC solver for %d induced charges\n", ncount);
if (logfile) fprintf(logfile,"BEM/ICC solver for %d induced charges\n", ncount);
}
if (comm->me == 0) utils::logmesg(lmp, "BEM/ICC solver for {} induced charges\n", ncount);
// initialize random induced charges with zero sum
@ -129,16 +121,16 @@ void FixPolarizeBEMICC::init()
int *mask = atom->mask;
int nlocal = atom->nlocal;
RanPark *random = new RanPark(lmp,seed_charge + comm->me);
RanPark *random = new RanPark(lmp, seed_charge + comm->me);
for (i = 0; i < 100; i++) random->uniform();
double sum,tmp = 0;
double sum, tmp = 0;
for (i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
q[i] = ave_charge*(random->uniform() - 0.5);
q[i] = ave_charge * (random->uniform() - 0.5);
tmp += q[i];
}
MPI_Allreduce(&tmp,&sum,1,MPI_DOUBLE,MPI_SUM,world);
sum /= (double)ncount;
MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
sum /= (double) ncount;
tmp = 0;
for (i = 0; i < nlocal; i++) {
@ -146,11 +138,10 @@ void FixPolarizeBEMICC::init()
q[i] -= sum;
tmp += q[i];
}
MPI_Allreduce(&tmp,&sum,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
delete random;
}
}
/* ---------------------------------------------------------------------- */
@ -159,37 +150,39 @@ void FixPolarizeBEMICC::setup(int vflag)
{
// check if the pair styles in use are compatible
if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric") == 0)
efield_pair = ((PairLJCutCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/msm/dielectric") == 0)
efield_pair = ((PairLJCutCoulMSMDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric") == 0)
efield_pair = ((PairLJCutCoulCutDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulCutDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"coul/long/dielectric") == 0)
efield_pair = ((PairCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"coul/cut/dielectric") == 0)
efield_pair = ((PairCoulCutDielectric*)force->pair)->efield;
else error->all(FLERR,"Pair style not compatible with fix polarize/bem/icc");
if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0)
efield_pair = ((PairLJCutCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = ((PairLJCutCoulMSMDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = ((PairLJCutCoulCutDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulCutDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = ((PairCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = ((PairCoulCutDielectric *) force->pair)->efield;
else
error->all(FLERR, "Pair style not compatible with fix polarize/bem/icc");
// check if kspace is used for force computation
if (force->kspace) {
kspaceflag = 1;
if (strcmp(force->kspace_style,"pppm/dielectric") == 0)
efield_kspace = ((PPPMDielectric*)force->kspace)->efield;
else if (strcmp(force->kspace_style,"msm/dielectric") == 0)
efield_kspace = ((MSMDielectric*)force->kspace)->efield;
else error->all(FLERR,"Kspace style not compatible with fix polarize/bem/icc");
if (strcmp(force->kspace_style, "pppm/dielectric") == 0)
efield_kspace = ((PPPMDielectric *) force->kspace)->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = ((MSMDielectric *) force->kspace)->efield;
else
error->all(FLERR, "Kspace style not compatible with fix polarize/bem/icc");
} else {
if (kspaceflag == 1) { // users specified kspace yes
error->warning(FLERR,"No Kspace style available for fix polarize/bem/icc");
if (kspaceflag == 1) { // users specified kspace yes
error->warning(FLERR, "No Kspace style available for fix polarize/bem/icc");
kspaceflag = 0;
}
}
@ -239,8 +232,8 @@ void FixPolarizeBEMICC::compute_induced_charges()
// Let's choose that epsilon[i] = em[i] for the interface particles
force_clear();
force->pair->compute(eflag,vflag);
if (kspaceflag) force->kspace->compute(eflag,vflag);
force->pair->compute(eflag, vflag);
if (kspaceflag) force->kspace->compute(eflag, vflag);
if (force->newton) comm->reverse_comm();
int i10 = 0;
@ -257,7 +250,7 @@ void FixPolarizeBEMICC::compute_induced_charges()
}
// divide (Ex,Ey,Ez) by epsilon[i] here
double dot = (Ex*norm[i][0] + Ey*norm[i][1] + Ez*norm[i][2]) / (2*MY_PI) / epsilon[i];
double dot = (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / (2 * MY_PI) / epsilon[i];
double q_free = q_real[i];
double q_bound = 0;
q_bound = (1.0 / em[i] - 1) * q_free - epsilon0 * (ed[i] / (2 * em[i])) * dot * area[i];
@ -271,8 +264,8 @@ void FixPolarizeBEMICC::compute_induced_charges()
for (itr = 0; itr < itr_max; itr++) {
force_clear();
force->pair->compute(eflag,vflag);
if (kspaceflag) force->kspace->compute(eflag,vflag);
force->pair->compute(eflag, vflag);
if (kspaceflag) force->kspace->compute(eflag, vflag);
if (force->newton) comm->reverse_comm();
double tol = 0;
@ -294,10 +287,10 @@ void FixPolarizeBEMICC::compute_induced_charges()
// note the area[i] is included here to ensure correct charge unit
// for direct use in force/efield compute
double dot = (Ex*norm[i][0] + Ey*norm[i][1] + Ez*norm[i][2]) / (4*MY_PI) / epsilon[i];
double dot = (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / (4 * MY_PI) / epsilon[i];
double q_bound = q[i] - q_free;
q_bound = (1 - omega) * q_bound + omega * ((1.0 / em[i] - 1) * q_free -
epsilon0 * (ed[i] / em[i]) * dot * area[i]);
q_bound = (1 - omega) * q_bound +
omega * ((1.0 / em[i] - 1) * q_free - epsilon0 * (ed[i] / em[i]) * dot * area[i]);
q[i] = q_free + q_bound;
// Eq. (11) in Tyagi et al., with f from Eq. (6)
@ -314,20 +307,20 @@ void FixPolarizeBEMICC::compute_induced_charges()
//q[i] = (1 - omega) * q[i] - omega * epsilon0 * f * dot * area[i];
double delta = fabs(qtmp - q_bound);
double r = (fabs(qtmp) > 0) ? delta/fabs(qtmp) : 0;
double r = (fabs(qtmp) > 0) ? delta / fabs(qtmp) : 0;
if (tol < r) tol = r;
#ifdef _POLARIZE_DEBUG
//printf("i = %d: q_bound = %f \n", i, q_bound);
#endif
#ifdef _POLARIZE_DEBUG
//printf("i = %d: q_bound = %f \n", i, q_bound);
#endif
}
comm->forward_comm_fix(this);
MPI_Allreduce(&tol,&rho,1,MPI_DOUBLE,MPI_MAX,world);
#ifdef _POLARIZE_DEBUG
MPI_Allreduce(&tol, &rho, 1, MPI_DOUBLE, MPI_MAX, world);
#ifdef _POLARIZE_DEBUG
printf("itr = %d: rho = %f\n", itr, rho);
#endif
#endif
if (itr > 0 && rho < tol_rel) break;
}
@ -342,9 +335,9 @@ void FixPolarizeBEMICC::force_clear()
if (force->newton) nbytes += sizeof(double) * atom->nghost;
if (nbytes) {
memset(&atom->f[0][0],0,3*nbytes);
if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
if (extraflag) atom->avec->force_clear(0,nbytes);
memset(&atom->f[0][0], 0, 3 * nbytes);
if (torqueflag) memset(&atom->torque[0][0], 0, 3 * nbytes);
if (extraflag) atom->avec->force_clear(0, nbytes);
}
}
@ -354,46 +347,49 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"itr_max") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
itr_max = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (strcmp(arg[iarg], "itr_max") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
itr_max = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"omega") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
omega = utils::numeric(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "omega") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
omega = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) kspaceflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) kspaceflag = 0;
else error->all(FLERR,"Illegal fix_modify command for fix polarize");
} else if (strcmp(arg[iarg], "kspace") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
if (strcmp(arg[iarg + 1], "yes") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg + 1], "no") == 0)
kspaceflag = 0;
else
error->all(FLERR, "Illegal fix_modify command for fix polarize");
iarg += 2;
} else if (strcmp(arg[iarg],"dielectrics") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix_modify command");
double epsiloni=-1, areai=-1;
double qunscaledi=0;
int set_charge=0;
double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (strcmp(arg[iarg+3],"NULL") != 0)
epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (strcmp(arg[iarg+4],"NULL") != 0)
areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
if (strcmp(arg[iarg+5],"NULL") != 0) {
qunscaledi = utils::numeric(FLERR,arg[iarg+5],false,lmp);
} else if (strcmp(arg[iarg], "dielectrics") == 0) {
if (iarg + 6 > narg) error->all(FLERR, "Illegal fix_modify command");
double epsiloni = -1, areai = -1;
double qunscaledi = 0;
int set_charge = 0;
double ediff = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
double emean = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
if (strcmp(arg[iarg + 3], "NULL") != 0)
epsiloni = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (strcmp(arg[iarg + 4], "NULL") != 0)
areai = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
if (strcmp(arg[iarg + 5], "NULL") != 0) {
qunscaledi = utils::numeric(FLERR, arg[iarg + 5], false, lmp);
set_charge = 1;
}
set_dielectric_params(ediff, emean, epsiloni, areai, set_charge,
qunscaledi);
set_dielectric_params(ediff, emean, epsiloni, areai, set_charge, qunscaledi);
iarg += 6;
} else if (strcmp(arg[iarg],"rand") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix_modify command");
ave_charge = utils::numeric(FLERR,arg[iarg+1],false,lmp);
seed_charge = utils::numeric(FLERR,arg[iarg+2],false,lmp);
} else if (strcmp(arg[iarg], "rand") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal fix_modify command");
ave_charge = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
seed_charge = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
randomized = 1;
iarg += 3;
} else error->all(FLERR,"Illegal fix_modify command");
} else
error->all(FLERR, "Illegal fix_modify command");
}
return iarg;
@ -401,8 +397,7 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
/* ---------------------------------------------------------------------- */
int FixPolarizeBEMICC::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int FixPolarizeBEMICC::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int m;
for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];
@ -421,8 +416,8 @@ void FixPolarizeBEMICC::unpack_forward_comm(int n, int first, double *buf)
set dielectric params for the atoms in the group
------------------------------------------------------------------------- */
void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean,
double epsiloni, double areai, int set_charge, double qvalue)
void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean, double epsiloni,
double areai, int set_charge, double qvalue)
{
double *area = atom->area;
double *ed = atom->ed;
@ -450,7 +445,10 @@ void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean,
double FixPolarizeBEMICC::compute_vector(int n)
{
if (n == 0) return iterations;
else if (n == 1) return rho;
else return 0;
if (n == 0)
return iterations;
else if (n == 1)
return rho;
else
return 0;
}

359
src/USER-DIELECTRIC/fix_polarize_functional.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
@ -62,7 +61,7 @@ using namespace FixConst;
using namespace MathExtra;
using namespace MathConst;
enum {REAL2SCALED=0,SCALED2REAL=1};
enum { REAL2SCALED = 0, SCALED2REAL = 1 };
#define EPSILON 1e-6
@ -71,18 +70,18 @@ enum {REAL2SCALED=0,SCALED2REAL=1};
/* ---------------------------------------------------------------------- */
FixPolarizeFunctional::FixPolarizeFunctional(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
Fix(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal fix polarize/functional command");
if (narg < 4) error->all(FLERR, "Illegal fix polarize/functional command");
avec = (AtomVecDielectric *) atom->style_match("dielectric");
if (!avec) error->all(FLERR,"Fix polarize/functional requires atom style dielectric");
if (!avec) error->all(FLERR, "Fix polarize/functional requires atom style dielectric");
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery < 0) error->all(FLERR,"Illegal fix polarize/functional command");
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
if (nevery < 0) error->all(FLERR, "Illegal fix polarize/functional command");
tolerance = EPSILON;
if (narg == 5) tolerance = utils::numeric(FLERR,arg[4],false,lmp);
if (narg == 5) tolerance = utils::numeric(FLERR, arg[4], false, lmp);
comm_forward = 1;
nmax = 0;
@ -131,7 +130,7 @@ FixPolarizeFunctional::FixPolarizeFunctional(LAMMPS *lmp, int narg, char **arg)
cg_A = nullptr;
grow_arrays(atom->nmax);
atom->add_callback(0); // to ensure to work with atom->sort()
atom->add_callback(0); // to ensure to work with atom->sort()
}
/* ---------------------------------------------------------------------- */
@ -151,7 +150,7 @@ FixPolarizeFunctional::~FixPolarizeFunctional()
memory->destroy(buffer2);
if (allocated) deallocate();
atom->delete_callback(id,0);
atom->delete_callback(id, 0);
}
/* ---------------------------------------------------------------------- */
@ -169,7 +168,7 @@ void FixPolarizeFunctional::init()
{
// mapping induced charge matrix/vector to atom tags and vice versa
int i,maxtag;
int i, maxtag;
double *q = atom->q;
int *mask = atom->mask;
tagint *tag = atom->tag;
@ -182,26 +181,28 @@ void FixPolarizeFunctional::init()
max_tag = -1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) max_tag = MAX(max_tag,tag[i]);
if (mask[i] & groupbit) max_tag = MAX(max_tag, tag[i]);
MPI_Allreduce(&max_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
MPI_Allreduce(&max_tag, &itmp, 1, MPI_LMP_TAGINT, MPI_MAX, world);
maxtag = (int) itmp;
memory->create(ncount,maxtag+1,"polarize:ncount");
memory->create(ncount, maxtag + 1, "polarize:ncount");
for (i = 0; i <= maxtag; i++) ncount[i] = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) ncount[tag[i]]++;
memory->create(tag2mat,maxtag+1,"polarize:tag2mat");
MPI_Allreduce(ncount,tag2mat,maxtag+1,MPI_INT,MPI_SUM,world);
memory->create(tag2mat, maxtag + 1, "polarize:tag2mat");
MPI_Allreduce(ncount, tag2mat, maxtag + 1, MPI_INT, MPI_SUM, world);
num_induced_charges = 0;
for (i = 0; i <= maxtag; i++)
if (tag2mat[i]) tag2mat[i] = num_induced_charges++;
else tag2mat[i] = -1;
if (tag2mat[i])
tag2mat[i] = num_induced_charges++;
else
tag2mat[i] = -1;
memory->create(mat2tag,num_induced_charges,"polarize:mat2tag");
memory->create(mat2tag, num_induced_charges, "polarize:mat2tag");
num_induced_charges = 0;
for (i = 0; i <= maxtag; i++)
@ -218,34 +219,34 @@ void FixPolarizeFunctional::init()
max_tag = -1;
for (i = 0; i < nlocal; i++)
if (!(mask[i] & groupbit)) max_tag = MAX(max_tag,tag[i]);
if (!(mask[i] & groupbit)) max_tag = MAX(max_tag, tag[i]);
MPI_Allreduce(&max_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
MPI_Allreduce(&max_tag, &itmp, 1, MPI_LMP_TAGINT, MPI_MAX, world);
maxtag = (int) itmp;
memory->create(ncount,maxtag+1,"polarize:ncount");
memory->create(ncount, maxtag + 1, "polarize:ncount");
for (i = 0; i <= maxtag; i++) ncount[i] = 0;
for (i = 0; i < nlocal; i++)
if (!(mask[i] & groupbit)) ncount[tag[i]]++;
memory->create(tag2mat_ions,maxtag+1,"polarize:tag2mat_ions");
MPI_Allreduce(ncount,tag2mat_ions,maxtag+1,MPI_INT,MPI_SUM,world);
memory->create(tag2mat_ions, maxtag + 1, "polarize:tag2mat_ions");
MPI_Allreduce(ncount, tag2mat_ions, maxtag + 1, MPI_INT, MPI_SUM, world);
num_ions = 0;
for (i = 0; i <= maxtag; i++)
if (tag2mat_ions[i]) tag2mat_ions[i] = num_ions++;
else tag2mat_ions[i] = -1;
if (tag2mat_ions[i])
tag2mat_ions[i] = num_ions++;
else
tag2mat_ions[i] = -1;
memory->create(mat2tag_ions,num_ions,"polarize:mat2tag_ions");
memory->create(rhs1,num_induced_charges,"polarize:rhs1");
memory->create(rhs2,num_induced_charges,"polarize:rhs2");
memory->create(mat2tag_ions, num_ions, "polarize:mat2tag_ions");
memory->create(rhs1, num_induced_charges, "polarize:rhs1");
memory->create(rhs2, num_induced_charges, "polarize:rhs2");
int buffer_size = (num_induced_charges > num_ions) ? num_induced_charges : num_ions;
memory->create(buffer1,buffer_size,num_induced_charges,"polarize:buffer1");
memory->create(buffer2,num_induced_charges,num_induced_charges,"polarize:buffer2");
memory->create(induced_charges,num_induced_charges,"polarize:induced_charges");
memory->create(buffer1, buffer_size, num_induced_charges, "polarize:buffer1");
memory->create(buffer2, num_induced_charges, num_induced_charges, "polarize:buffer2");
memory->create(induced_charges, num_induced_charges, "polarize:induced_charges");
num_ions = 0;
for (i = 0; i <= maxtag; i++)
@ -267,22 +268,21 @@ void FixPolarizeFunctional::init()
// need a full neighbor list w/ Newton off and ghost neighbors
// built whenever re-neighboring occurs
int irequest = neighbor->request(this,instance_me);
int irequest = neighbor->request(this, instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->occasional = 0;
if (force->kspace) g_ewald = force->kspace->g_ewald;
else g_ewald = 0.01;
if (force->kspace)
g_ewald = force->kspace->g_ewald;
else
g_ewald = 0.01;
if (comm->me == 0) {
if (screen) fprintf(screen,"Direct solver using the variational approach "
"for %d induced charges\n", num_induced_charges);
if (logfile) fprintf(logfile,"Direct solver using the variational approach "
"for %d induced charges\n", num_induced_charges);
}
if (comm->me == 0)
utils::logmesg(lmp, "Direct solver using a variational approach for {} induced charges\n",
num_induced_charges);
}
/* ---------------------------------------------------------------------- */
@ -298,35 +298,37 @@ void FixPolarizeFunctional::setup(int vflag)
{
// check if the pair styles in use are compatible
if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric") == 0)
efield_pair = ((PairLJCutCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/msm/dielectric") == 0)
efield_pair = ((PairLJCutCoulMSMDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric") == 0)
efield_pair = ((PairLJCutCoulCutDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulCutDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"coul/long/dielectric") == 0)
efield_pair = ((PairCoulLongDielectric*)force->pair)->efield;
else if (strcmp(force->pair_style,"coul/cut/dielectric") == 0)
efield_pair = ((PairCoulCutDielectric*)force->pair)->efield;
else error->all(FLERR,"Pair style not compatible with fix polarize/functional");
if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0)
efield_pair = ((PairLJCutCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = ((PairLJCutCoulMSMDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = ((PairLJCutCoulCutDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = ((PairLJCutCoulCutDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = ((PairCoulLongDielectric *) force->pair)->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = ((PairCoulCutDielectric *) force->pair)->efield;
else
error->all(FLERR, "Pair style not compatible with fix polarize/functional");
if (force->kspace) {
kspaceflag = 1;
if (strcmp(force->kspace_style,"pppm/dielectric") == 0)
efield_kspace = ((PPPMDielectric*)force->kspace)->efield;
else if (strcmp(force->kspace_style,"msm/dielectric") == 0)
efield_kspace = ((MSMDielectric*)force->kspace)->efield;
else error->all(FLERR,"Kspace style not compatible with fix polarize/functional");
if (strcmp(force->kspace_style, "pppm/dielectric") == 0)
efield_kspace = ((PPPMDielectric *) force->kspace)->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = ((MSMDielectric *) force->kspace)->efield;
else
error->all(FLERR, "Kspace style not compatible with fix polarize/functional");
} else {
if (kspaceflag == 1) { // users specified kspace yes
error->warning(FLERR,"No Kspace style available for fix polarize/functional");
if (kspaceflag == 1) { // users specified kspace yes
error->warning(FLERR, "No Kspace style available for fix polarize/functional");
kspaceflag = 0;
}
}
@ -371,8 +373,7 @@ void FixPolarizeFunctional::update_induced_charges()
// conjugate gradient solver for w from Rww * w = -qRqw
for (int i = 0; i < num_induced_charges; i++)
for (int j = 0; j < num_induced_charges; j++)
cg_A[i][j] = Rww[i][j] + Rww[j][i];
for (int j = 0; j < num_induced_charges; j++) cg_A[i][j] = Rww[i][j] + Rww[j][i];
for (int i = 0; i < num_induced_charges; i++) induced_charges[i] = 0;
@ -386,7 +387,7 @@ void FixPolarizeFunctional::update_induced_charges()
for (int i = 0; i < nlocal; i++) {
if (induced_charge_idx[i] < 0) continue;
int idx = induced_charge_idx[i];
q[i] = -induced_charges[idx] / (4*MY_PI);
q[i] = -induced_charges[idx] / (4 * MY_PI);
}
// revert to scaled charges to calculate forces
@ -401,9 +402,9 @@ void FixPolarizeFunctional::update_induced_charges()
void FixPolarizeFunctional::charge_rescaled(int scaled2real)
{
double* q = atom->q;
double* q_real = atom->q_unscaled;
double* epsilon = atom->epsilon;
double *q = atom->q;
double *q_real = atom->q_unscaled;
double *epsilon = atom->epsilon;
int nlocal = atom->nlocal;
if (scaled2real) {
@ -478,32 +479,35 @@ int FixPolarizeFunctional::modify_param(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) kspaceflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) kspaceflag = 0;
else error->all(FLERR,"Illegal fix_modify command for fix polarize/functional");
if (strcmp(arg[iarg], "kspace") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
if (strcmp(arg[iarg + 1], "yes") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg + 1], "no") == 0)
kspaceflag = 0;
else
error->all(FLERR, "Illegal fix_modify command for fix polarize/functional");
iarg += 2;
} else if (strcmp(arg[iarg],"dielectrics") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix_modify command");
double epsiloni=-1, areai=-1;
double q_unscaled=0;
int set_charge=0;
double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (strcmp(arg[iarg+3],"nullptr") != 0)
epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (strcmp(arg[iarg+4],"nullptr") != 0)
areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
if (strcmp(arg[iarg+5],"nullptr") != 0) {
q_unscaled = utils::numeric(FLERR,arg[iarg+5],false,lmp);
} else if (strcmp(arg[iarg], "dielectrics") == 0) {
if (iarg + 6 > narg) error->all(FLERR, "Illegal fix_modify command");
double epsiloni = -1, areai = -1;
double q_unscaled = 0;
int set_charge = 0;
double ediff = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
double emean = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
if (strcmp(arg[iarg + 3], "nullptr") != 0)
epsiloni = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (strcmp(arg[iarg + 4], "nullptr") != 0)
areai = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
if (strcmp(arg[iarg + 5], "nullptr") != 0) {
q_unscaled = utils::numeric(FLERR, arg[iarg + 5], false, lmp);
set_charge = 1;
}
set_dielectric_params(ediff, emean, epsiloni, areai, set_charge,
q_unscaled);
set_dielectric_params(ediff, emean, epsiloni, areai, set_charge, q_unscaled);
iarg += 6;
} else error->all(FLERR,"Illegal fix_modify command");
} else
error->all(FLERR, "Illegal fix_modify command");
}
return iarg;
@ -533,8 +537,7 @@ int FixPolarizeFunctional::unpack_exchange(int nlocal, double *buf)
/* ---------------------------------------------------------------------- */
int FixPolarizeFunctional::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int FixPolarizeFunctional::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int m;
for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];
@ -556,8 +559,8 @@ void FixPolarizeFunctional::unpack_forward_comm(int n, int first, double *buf)
void FixPolarizeFunctional::grow_arrays(int n)
{
if (n > nmax) nmax = n;
memory->grow(induced_charge_idx,nmax,"fix:induced_charge_idx");
memory->grow(ion_idx,nmax,"fix:ion_idx");
memory->grow(induced_charge_idx, nmax, "fix:induced_charge_idx");
memory->grow(ion_idx, nmax, "fix:ion_idx");
}
/* ----------------------------------------------------------------------
@ -587,21 +590,21 @@ void FixPolarizeFunctional::set_arrays(int i)
double FixPolarizeFunctional::memory_usage()
{
double bytes = 0;
bytes += num_induced_charges*num_induced_charges*sizeof(double); // inverse_matrix
bytes += num_induced_charges*num_induced_charges*sizeof(double); // Rww
bytes += num_induced_charges*num_induced_charges*sizeof(double); // G1ww
bytes += num_induced_charges*num_induced_charges*sizeof(double); // ndotGww
bytes += num_induced_charges*num_induced_charges*sizeof(double); // G2ww
bytes += num_induced_charges*num_induced_charges*sizeof(double); // G3ww
bytes += num_induced_charges*sizeof(double); // qiRqwVector
bytes += num_induced_charges*sizeof(double); // sum2G2wq
bytes += num_induced_charges*sizeof(double); // sum1G2qw
bytes += num_induced_charges*sizeof(double); // sum1G1qw_epsilon
bytes += num_induced_charges*sizeof(double); // sum2ndotGwq_epsilon
bytes += num_ions*num_induced_charges*sizeof(double); // G1qw_real
bytes += nmax*sizeof(int); // induced_charge_idx
bytes += nmax*sizeof(int); // ion_idx
bytes += num_induced_charges*sizeof(double); // induced_charges
bytes += num_induced_charges * num_induced_charges * sizeof(double); // inverse_matrix
bytes += num_induced_charges * num_induced_charges * sizeof(double); // Rww
bytes += num_induced_charges * num_induced_charges * sizeof(double); // G1ww
bytes += num_induced_charges * num_induced_charges * sizeof(double); // ndotGww
bytes += num_induced_charges * num_induced_charges * sizeof(double); // G2ww
bytes += num_induced_charges * num_induced_charges * sizeof(double); // G3ww
bytes += num_induced_charges * sizeof(double); // qiRqwVector
bytes += num_induced_charges * sizeof(double); // sum2G2wq
bytes += num_induced_charges * sizeof(double); // sum1G2qw
bytes += num_induced_charges * sizeof(double); // sum1G1qw_epsilon
bytes += num_induced_charges * sizeof(double); // sum2ndotGwq_epsilon
bytes += num_ions * num_induced_charges * sizeof(double); // G1qw_real
bytes += nmax * sizeof(int); // induced_charge_idx
bytes += nmax * sizeof(int); // ion_idx
bytes += num_induced_charges * sizeof(double); // induced_charges
return bytes;
}
@ -611,7 +614,7 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
{
int *mask = atom->mask;
int *type = atom->type;
tagint* tag = atom->tag;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
double **x = atom->x;
double *area = atom->area;
@ -622,7 +625,7 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
// invoke full neighbor list
int inum,jnum,*ilist,*jlist,*numneigh,**firstneigh;
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -657,7 +660,7 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
double delx = xtmp - x[k][0];
double dely = ytmp - x[k][1];
double delz = ztmp - x[k][2];
domain->minimum_image(delx,dely,delz);
domain->minimum_image(delx, dely, delz);
int mk = tag2mat[tag[k]];
// G1ww[mi][mk] = calculate_greens_ewald(delx, dely, delz);
@ -668,7 +671,7 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
calculate_grad_greens_ewald(gradG1ww, delx, dely, delz);
// use mu to store the normal vector of interface vertex
buffer2[mi][mk] = MathExtra::dot3(norm[i], gradG1ww) / (4*MY_PI);
buffer2[mi][mk] = MathExtra::dot3(norm[i], gradG1ww) / (4 * MY_PI);
}
}
@ -676,15 +679,14 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
// even though in the above loop there is mk == mi
buffer1[mi][mi] = calculate_greens_ewald_self_vertex(area[i]);
buffer2[mi][mi] = calculate_ndotgreens_ewald_self_vertex(area[i], curvature[i]) /
(4*MY_PI);
buffer2[mi][mi] = calculate_ndotgreens_ewald_self_vertex(area[i], curvature[i]) / (4 * MY_PI);
}
}
MPI_Allreduce(buffer1[0],G1ww[0],num_induced_charges*num_induced_charges,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer2[0],ndotGww[0],num_induced_charges*num_induced_charges,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer1[0], G1ww[0], num_induced_charges * num_induced_charges, MPI_DOUBLE, MPI_SUM,
world);
MPI_Allreduce(buffer2[0], ndotGww[0], num_induced_charges * num_induced_charges, MPI_DOUBLE,
MPI_SUM, world);
// calculate G2ww
// fill up buffer1 with local G2ww
@ -734,8 +736,8 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
}
}
MPI_Allreduce(buffer1[0],G2ww[0],num_induced_charges*num_induced_charges,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer1[0], G2ww[0], num_induced_charges * num_induced_charges, MPI_DOUBLE, MPI_SUM,
world);
// calculate G3ww and Rww
// G3ww is implemented as in _exact(), but can be optionally excluded
@ -815,18 +817,18 @@ void FixPolarizeFunctional::calculate_Rww_cutoff()
}
}
MPI_Allreduce(buffer1[0],Rww[0],num_induced_charges*num_induced_charges,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(buffer1[0], Rww[0], num_induced_charges * num_induced_charges, MPI_DOUBLE, MPI_SUM,
world);
#ifdef _POLARIZE_DEBUG
#ifdef _POLARIZE_DEBUG
if (comm->me == 0) {
FILE* fp = fopen("Rww-functional.txt", "w");
FILE *fp = fopen("Rww-functional.txt", "w");
for (int i = 0; i < num_induced_charges; i++)
fprintf(fp, "%d %g %g %g\n", i, Rww[i][i], Rww[i][num_induced_charges/2],
Rww[num_induced_charges/2][i]);
fprintf(fp, "%d %g %g %g\n", i, Rww[i][i], Rww[i][num_induced_charges / 2],
Rww[num_induced_charges / 2][i]);
fclose(fp);
}
#endif
#endif
}
/* ---------------------------------------------------------------------- */
@ -849,7 +851,7 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
// invoke full neighbor list
int inum,*ilist,*jlist,*numneigh,**firstneigh;
int inum, *ilist, *jlist, *numneigh, **firstneigh;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -879,7 +881,7 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
delx = xtmp - x[k][0];
dely = ytmp - x[k][1];
delz = ztmp - x[k][2];
domain->minimum_image(delx,dely,delz);
domain->minimum_image(delx, dely, delz);
r = sqrt(delx * delx + dely * dely + delz * delz);
int mk = tag2mat[tag[k]];
@ -890,18 +892,18 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
}
}
MPI_Allreduce(&buffer1[0][0],&G1qw_real[0][0],num_ions*num_induced_charges,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&buffer1[0][0], &G1qw_real[0][0], num_ions * num_induced_charges, MPI_DOUBLE,
MPI_SUM, world);
// the following loop need the above results: gradG1wq_real
// calculate sum1G1qw_epsilon and sum2ndotGwq_epsilon
// fill up rhs1 with local sum1G1qw_epsilon and rhs2 with local sum2ndotGwq_epsilon
memset(rhs1, 0, num_induced_charges*sizeof(double));
memset(rhs2, 0, num_induced_charges*sizeof(double));
memset(rhs1, 0, num_induced_charges * sizeof(double));
memset(rhs2, 0, num_induced_charges * sizeof(double));
for (kk = 0; kk < inum; kk++) {
k = ilist[kk]; // k is local index
k = ilist[kk]; // k is local index
if (mask[k] & groupbit) {
// interface particles
int mk = induced_charge_idx[k];
@ -920,9 +922,9 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
delx = x[i][0] - xtmp;
dely = x[i][1] - ytmp;
delz = x[i][2] - ztmp;
domain->minimum_image(delx,dely,delz);
domain->minimum_image(delx, dely, delz);
int mi = tag2mat_ions[tag[i]]; //ion_idx[i];
int mi = tag2mat_ions[tag[i]]; //ion_idx[i];
//calculate_grad_greens_real(gradG1wq_real[mk][mi], delx, dely, delz);
double gradG1wq[3];
@ -944,18 +946,18 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
}
}
MPI_Allreduce(rhs1,sum1G1qw_epsilon,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(rhs2,sum2ndotGwq_epsilon,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(rhs1, sum1G1qw_epsilon, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(rhs2, sum2ndotGwq_epsilon, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
// calculate G2, gradient G2
// sum2G2wq and sum1G2qw
// for (int i = 0; i < num_induced_charges; i++) rhs1[i] = rhs2[i] = 0;
memset(rhs1, 0, num_induced_charges*sizeof(double));
memset(rhs2, 0, num_induced_charges*sizeof(double));
// for (int i = 0; i < num_induced_charges; i++) rhs1[i] = rhs2[i] = 0;
memset(rhs1, 0, num_induced_charges * sizeof(double));
memset(rhs2, 0, num_induced_charges * sizeof(double));
for (kk = 0; kk < inum; kk++) {
k = ilist[kk]; // k is local index
k = ilist[kk]; // k is local index
if (mask[k] & groupbit) {
// interface particles
int mk = induced_charge_idx[k];
@ -985,16 +987,16 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
}
}
MPI_Allreduce(rhs1,sum2G2wq,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(rhs2,sum1G2qw,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(rhs1, sum2G2wq, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(rhs2, sum1G2qw, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
// calculate G3, gradient G3
// fill up rhs1 with local qiRqwVector
memset(rhs1, 0, num_induced_charges*sizeof(double));
memset(rhs1, 0, num_induced_charges * sizeof(double));
for (kk = 0; kk < inum; kk++) {
k = ilist[kk]; // k is local index
k = ilist[kk]; // k is local index
if (mask[k] & groupbit) {
// interface particles
int mk = induced_charge_idx[k];
@ -1011,7 +1013,7 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
sum1G3qw += sum2ndotGwq_epsilon[mi] * G2ww[mi][mk] * area[i] * ed[i];
} else {
// ions particles: i can be ghost atoms
int mi = tag2mat_ions[tag[i]]; //ion_idx[i];
int mi = tag2mat_ions[tag[i]]; //ion_idx[i];
qiRwwVectorTemp1 += q[i] * (1.0 - em[k] / epsilon[i]) * G1qw_real[mi][mk];
}
}
@ -1021,32 +1023,30 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
sum1G3qw += sum2ndotGwq_epsilon[mk] * G2ww[mk][mk] * area[k] * ed[k];
// qiRwwVectorTemp2 is a significant contribution, of which sum2G2wq is significant
double qiRwwVectorTemp2 = (1.0 - 2.0 * em[k]) * sum2G2wq[mk] +
sum1G2qw[mk] + 2.0 * sum1G3qw;
double qiRwwVectorTemp2 = (1.0 - 2.0 * em[k]) * sum2G2wq[mk] + sum1G2qw[mk] + 2.0 * sum1G3qw;
// qiRqwVector[mk] = qiRwwVectorTemp1 + qiRwwVectorTemp2;
rhs1[mk] = qiRwwVectorTemp1 + qiRwwVectorTemp2;
}
}
MPI_Allreduce(rhs1,qiRqwVector,num_induced_charges,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(rhs1, qiRqwVector, num_induced_charges, MPI_DOUBLE, MPI_SUM, world);
#ifdef _POLARIZE_DEBUG
#ifdef _POLARIZE_DEBUG
if (comm->me == 0) {
FILE* fp = fopen("qRqw-functional.txt", "w");
for (int i = 0; i < num_induced_charges; i++)
fprintf(fp, "%d %g\n", i, qiRqwVector[i]);
FILE *fp = fopen("qRqw-functional.txt", "w");
for (int i = 0; i < num_induced_charges; i++) fprintf(fp, "%d %g\n", i, qiRqwVector[i]);
fclose(fp);
}
#endif
#endif
}
/* ----------------------------------------------------------------------
set dielectric params for the atom in the group
------------------------------------------------------------------------- */
void FixPolarizeFunctional::set_dielectric_params(double ediff, double emean,
double epsiloni, double areai, int set_charge, double qvalue)
void FixPolarizeFunctional::set_dielectric_params(double ediff, double emean, double epsiloni,
double areai, int set_charge, double qvalue)
{
double *area = atom->area;
double *ed = atom->ed;
@ -1081,16 +1081,14 @@ double FixPolarizeFunctional::greens_real(double r)
double FixPolarizeFunctional::grad_greens_real_factor(double r)
{
double alpharij = g_ewald * r;
double factor = erfc(alpharij) + 2.0 * alpharij / MY_PIS *
exp(-(alpharij * alpharij));
double r3 = r*r*r;
double factor = erfc(alpharij) + 2.0 * alpharij / MY_PIS * exp(-(alpharij * alpharij));
double r3 = r * r * r;
return (factor * (-1.0 / r3));
}
/* ---------------------------------------------------------------------- */
void FixPolarizeFunctional::calculate_grad_greens_real(double *vec,
double dx, double dy, double dz)
void FixPolarizeFunctional::calculate_grad_greens_real(double *vec, double dx, double dy, double dz)
{
double r = sqrt(dx * dx + dy * dy + dz * dz);
double real = grad_greens_real_factor(r);
@ -1101,8 +1099,7 @@ void FixPolarizeFunctional::calculate_grad_greens_real(double *vec,
/* ---------------------------------------------------------------------- */
double FixPolarizeFunctional::calculate_greens_ewald(double dx, double dy,
double dz)
double FixPolarizeFunctional::calculate_greens_ewald(double dx, double dy, double dz)
{
// excluding the reciprocal term
double r = sqrt(dx * dx + dy * dy + dz * dz);
@ -1111,8 +1108,8 @@ double FixPolarizeFunctional::calculate_greens_ewald(double dx, double dy,
/* ---------------------------------------------------------------------- */
void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec,
double dx, double dy, double dz)
void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec, double dx, double dy,
double dz)
{
// real part of grad greens, excluding the reciprocal term
calculate_grad_greens_real(vec, dx, dy, dz);
@ -1120,17 +1117,15 @@ void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec,
/* ---------------------------------------------------------------------- */
void FixPolarizeFunctional::calculate_matrix_multiply_vector(double **matrix,
double *in_vec, double *out_vec, int M)
void FixPolarizeFunctional::calculate_matrix_multiply_vector(double **matrix, double *in_vec,
double *out_vec, int M)
{
#if defined(OPENMP)
#pragma parallel omp for
#endif
for (int k = 0; k < M; ++k) {
double temp = 0.0;
for (int l = 0; l < M; ++l) {
temp += matrix[k][l] * in_vec[l];
}
for (int l = 0; l < M; ++l) { temp += matrix[k][l] * in_vec[l]; }
out_vec[k] = temp;
}
}
@ -1147,8 +1142,7 @@ double FixPolarizeFunctional::calculate_greens_ewald_self_vertex(double area)
/* ---------------------------------------------------------------------- */
double FixPolarizeFunctional::calculate_ndotgreens_ewald_self_vertex(double area,
double curvature)
double FixPolarizeFunctional::calculate_ndotgreens_ewald_self_vertex(double area, double curvature)
{
// this term is important, cannot be set to zero
// curvature = 1 / R, minus if norm is inverse of R to center.
@ -1161,16 +1155,16 @@ double FixPolarizeFunctional::calculate_ndotgreens_ewald_self_vertex(double area
where ^t is the transpose operator
-- ---------------------------------------------------------------------- */
double FixPolarizeFunctional::inner_product(double* x, double* y, int N)
double FixPolarizeFunctional::inner_product(double *x, double *y, int N)
{
double t = 0;
for (int i = 0; i < N; i++) t += x[i]*y[i];
for (int i = 0; i < N; i++) t += x[i] * y[i];
return t;
}
/* ---------------------------------------------------------------------- */
void FixPolarizeFunctional::cg_solver(double** A, double* b, double* x, int N)
void FixPolarizeFunctional::cg_solver(double **A, double *b, double *x, int N)
{
calculate_matrix_multiply_vector(A, x, cg_p, N);
for (int i = 0; i < N; i++) {
@ -1194,8 +1188,8 @@ void FixPolarizeFunctional::cg_solver(double** A, double* b, double* x, int N)
// x = x + alpha * p
// r = r - alpha * Ap
for (int i = 0; i < N; i++) {
x[i] = x[i] + alpha*cg_p[i];
cg_r[i] = cg_r[i] - alpha*cg_Ap[i];
x[i] = x[i] + alpha * cg_p[i];
cg_r[i] = cg_r[i] - alpha * cg_Ap[i];
}
double rsq_new = inner_product(cg_r, cg_r, N);
@ -1203,8 +1197,7 @@ void FixPolarizeFunctional::cg_solver(double** A, double* b, double* x, int N)
// beta = rsq_new / rsq
double beta = rsq_new / rsq;
for (int i = 0; i < N; i++)
cg_p[i] = cg_r[i] + beta*cg_p[i];
for (int i = 0; i < N; i++) cg_p[i] = cg_r[i] + beta * cg_p[i];
rsq = rsq_new;
}
}