diff --git a/src/USER-DIELECTRIC/fix_polarize_bem_gmres.cpp b/src/USER-DIELECTRIC/fix_polarize_bem_gmres.cpp index 296c9c3835..2fd3d14b3f 100644 --- a/src/USER-DIELECTRIC/fix_polarize_bem_gmres.cpp +++ b/src/USER-DIELECTRIC/fix_polarize_bem_gmres.cpp @@ -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; diff --git a/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp b/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp index 04f66a40e1..325e3044ff 100644 --- a/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp +++ b/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp @@ -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 @@ -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; } diff --git a/src/USER-DIELECTRIC/fix_polarize_functional.cpp b/src/USER-DIELECTRIC/fix_polarize_functional.cpp index 51d8dc9d94..55fcd98bbe 100644 --- a/src/USER-DIELECTRIC/fix_polarize_functional.cpp +++ b/src/USER-DIELECTRIC/fix_polarize_functional.cpp @@ -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; } }