diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp index 048dddf2a1..8cc2815453 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -37,8 +36,8 @@ #include "neigh_request.h" #include "neighbor.h" #include "pair.h" -#include "update.h" #include "universe.h" +#include "update.h" #include #include @@ -115,11 +114,10 @@ static int constexpr albemunu[21][4] = { /* ---------------------------------------------------------------------- */ -ComputeBornMatrix::ComputeBornMatrix(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), id_virial(nullptr), temp_x(nullptr), - temp_f(nullptr) +ComputeBornMatrix::ComputeBornMatrix(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), id_virial(nullptr), temp_x(nullptr), temp_f(nullptr) { - if (narg < 3) error->all(FLERR,"Illegal compute born/matrix command"); + if (narg < 3) error->all(FLERR, "Illegal compute born/matrix command"); MPI_Comm_rank(world, &me); @@ -142,125 +140,139 @@ ComputeBornMatrix::ComputeBornMatrix(LAMMPS *lmp, int narg, char **arg) : } else { int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"numdiff") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal compute born/matrix command"); + if (strcmp(arg[iarg], "numdiff") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal compute born/matrix command"); numflag = 1; - numdelta = utils::numeric(FLERR,arg[iarg+1],false,lmp); + numdelta = utils::numeric(FLERR, arg[iarg + 1], false, lmp); if (numdelta <= 0.0) error->all(FLERR, "Illegal compute born/matrix command"); - id_virial = utils::strdup(arg[iarg+2]); + id_virial = utils::strdup(arg[iarg + 2]); int icompute = modify->find_compute(id_virial); - if (icompute < 0) error->all(FLERR,"Could not find compute born/matrix pressure ID"); + if (icompute < 0) error->all(FLERR, "Could not find compute born/matrix pressure ID"); compute_virial = modify->compute[icompute]; if (compute_virial->pressflag == 0) - error->all(FLERR,"Compute born/matrix pressure ID does not compute pressure"); + error->all(FLERR, "Compute born/matrix pressure ID does not compute pressure"); iarg += 3; - } else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; - else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1; - else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1; - else if (strcmp(arg[iarg],"dihedral") == 0) dihedflag = 1; - else if (strcmp(arg[iarg],"improper") == 0) impflag = 1; - else error->all(FLERR,"Illegal compute born/matrix command"); + } else if (strcmp(arg[iarg], "pair") == 0) { + pairflag = 1; + } else if (strcmp(arg[iarg], "bond") == 0) { + bondflag = 1; + } else if (strcmp(arg[iarg], "angle") == 0) { + angleflag = 1; + } else if (strcmp(arg[iarg], "dihedral") == 0) { + dihedflag = 1; + } else if (strcmp(arg[iarg], "improper") == 0) { + impflag = 1; + } else { + error->all(FLERR, "Illegal compute born/matrix command"); + } ++iarg; } } if (pairflag) { - if (numflag) error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); if (force->pair) { - if (force->pair->born_matrix_enable == 0) { + if (force->pair->born_matrix_enable == 0) if (comm->me == 0) error->warning(FLERR, "Pair style does not support compute born/matrix"); - } - } else { - pairflag = 0; } - } - if (bondflag) { - if (numflag) error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); - if (force->bond) { - if (force->bond->born_matrix_enable == 0) { - if (comm->me == 0) error->warning(FLERR, "Bond style does not support compute born/matrix"); - } - } else { - bondflag = 0; - } - } - if (angleflag) { - if (numflag) error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); - if (force->angle) { - if (force->angle->born_matrix_enable == 0) { - if (comm->me == 0) error->warning(FLERR, "Angle style does not support compute born/matrix"); - } - } else { - angleflag = 0; - } - } - if (dihedflag) { - if (numflag) error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); - if (force->dihedral) { - if (force->dihedral->born_matrix_enable == 0) { - if (comm->me == 0) error->warning(FLERR, "Dihedral style does not support compute born/matrix"); - } - } else { - dihedflag = 0; - } - } - if (impflag) { - if (numflag) error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); - if (force->improper) { - if (force->improper->born_matrix_enable == 0) { - if (comm->me == 0) error->warning(FLERR, "Improper style does not support compute born/matrix"); - } - } else { - impflag = 0; - } - } - if (force->kspace) { - if (!numflag) error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix"); - } - - // Initialize some variables - - values_local = values_global = vector = nullptr; - - // this fix produces a global vector - - memory->create(vector, nvalues, "born_matrix:vector"); - memory->create(values_global, nvalues, "born_matrix:values_global"); - size_vector = nvalues; - - vector_flag = 1; - extvector = 0; - maxatom = 0; - - if (!numflag) { - memory->create(values_local, nvalues, "born_matrix:values_local"); } else { - - reallocate(); - - // set fixed-point to default = center of cell - - fixedpoint[0] = 0.5 * (domain->boxlo[0] + domain->boxhi[0]); - fixedpoint[1] = 0.5 * (domain->boxlo[1] + domain->boxhi[1]); - fixedpoint[2] = 0.5 * (domain->boxlo[2] + domain->boxhi[2]); - - // define the cartesian indices for each strain (Voigt order) - - dirlist[0][0] = 0; - dirlist[0][1] = 0; - dirlist[1][0] = 1; - dirlist[1][1] = 1; - dirlist[2][0] = 2; - dirlist[2][1] = 2; - - dirlist[3][0] = 1; - dirlist[3][1] = 2; - dirlist[4][0] = 0; - dirlist[4][1] = 2; - dirlist[5][0] = 0; - dirlist[5][1] = 1; + pairflag = 0; } } +if (bondflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->bond) { + if (force->bond->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Bond style does not support compute born/matrix"); + } + } else { + bondflag = 0; + } +} +if (angleflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->angle) { + if (force->angle->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Angle style does not support compute born/matrix"); + } + } else { + angleflag = 0; + } +} +if (dihedflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->dihedral) { + if (force->dihedral->born_matrix_enable == 0) { + if (comm->me == 0) + error->warning(FLERR, "Dihedral style does not support compute born/matrix"); + } + } else { + dihedflag = 0; + } +} +if (impflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->improper) { + if (force->improper->born_matrix_enable == 0) { + if (comm->me == 0) + error->warning(FLERR, "Improper style does not support compute born/matrix"); + } + } else { + impflag = 0; + } +} +if (force->kspace) { + if (!numflag) error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix"); +} + +// Initialize some variables + +values_local = values_global = vector = nullptr; + +// this fix produces a global vector + +memory->create(vector, nvalues, "born_matrix:vector"); +memory->create(values_global, nvalues, "born_matrix:values_global"); +size_vector = nvalues; + +vector_flag = 1; +extvector = 0; +maxatom = 0; + +if (!numflag) { + memory->create(values_local, nvalues, "born_matrix:values_local"); +} else { + + reallocate(); + + // set fixed-point to default = center of cell + + fixedpoint[0] = 0.5 * (domain->boxlo[0] + domain->boxhi[0]); + fixedpoint[1] = 0.5 * (domain->boxlo[1] + domain->boxhi[1]); + fixedpoint[2] = 0.5 * (domain->boxlo[2] + domain->boxhi[2]); + + // define the cartesian indices for each strain (Voigt order) + + dirlist[0][0] = 0; + dirlist[0][1] = 0; + dirlist[1][0] = 1; + dirlist[1][1] = 1; + dirlist[2][0] = 2; + dirlist[2][1] = 2; + + dirlist[3][0] = 1; + dirlist[3][1] = 2; + dirlist[4][0] = 0; + dirlist[4][1] = 2; + dirlist[5][0] = 0; + dirlist[5][1] = 1; +} +} /* ---------------------------------------------------------------------- */ @@ -297,8 +309,7 @@ void ComputeBornMatrix::init() // check for virial compute int icompute = modify->find_compute(id_virial); - if (icompute < 0) error->all(FLERR, - "Virial compute ID for compute born/matrix does not exist"); + if (icompute < 0) error->all(FLERR, "Virial compute ID for compute born/matrix does not exist"); compute_virial = modify->compute[icompute]; // set up reverse index lookup @@ -318,7 +329,6 @@ void ComputeBornMatrix::init() virialVtoV[3] = 5; virialVtoV[4] = 4; virialVtoV[5] = 3; - } } @@ -362,15 +372,12 @@ void ComputeBornMatrix::compute_vector() // convert from pressure to energy units - double inv_nktv2p = 1.0/force->nktv2p; + double inv_nktv2p = 1.0 / force->nktv2p; double volume = domain->xprd * domain->yprd * domain->zprd; - for (int m = 0; m < nvalues; m++) { - values_global[m] *= inv_nktv2p * volume; - } + for (int m = 0; m < nvalues; m++) { values_global[m] *= inv_nktv2p * volume; } } for (int m = 0; m < nvalues; m++) vector[m] = values_global[m]; - } /* ---------------------------------------------------------------------- @@ -451,23 +458,18 @@ void ComputeBornMatrix::compute_pairs() // Add contribution to Born tensor - pair->born_matrix(i, j, itype, jtype, rsq, factor_coul, - factor_lj, dupair, du2pair); + pair->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, dupair, du2pair); pair_pref = du2pair - dupair * rinv; // See albemunu in compute_born_matrix.h for indices order. - a = 0; - b = 0; - c = 0; - d = 0; + a = b = c = d = 0; for (int m = 0; m < nvalues; m++) { a = albemunu[m][0]; b = albemunu[m][1]; c = albemunu[m][2]; d = albemunu[m][3]; - values_local[m] += pair_pref * rij[a] * rij[b] * - rij[c] * rij[d] * r2inv; + values_local[m] += pair_pref * rij[a] * rij[b] * rij[c] * rij[d] * r2inv; } } } @@ -507,7 +509,7 @@ void ComputeBornMatrix::compute_numdiff() for (int idir = 0; idir < NDIR_VIRIAL; idir++) { // forward - + displace_atoms(nall, idir, 1.0); force_clear(nall); update_virial(); @@ -518,7 +520,7 @@ void ComputeBornMatrix::compute_numdiff() restore_atoms(nall, idir); // backward - + displace_atoms(nall, idir, -1.0); force_clear(nall); update_virial(); @@ -540,15 +542,13 @@ void ComputeBornMatrix::compute_numdiff() update_virial(); // add on virial terms - + virial_addon(); // restore original forces for owned and ghost atoms for (int i = 0; i < nall; i++) - for (int k = 0; k < 3; k++) - f[i][k] = temp_f[i][k]; - + for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k]; } /* ---------------------------------------------------------------------- @@ -561,11 +561,11 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) // NOTE: transposing k and l would seem to be equivalent but it is not // only this version matches analytic results for lj/cut - + int k = dirlist[idir][0]; int l = dirlist[idir][1]; for (int i = 0; i < nall; i++) - x[i][k] = temp_x[i][k] + numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); + x[i][k] = temp_x[i][k] + numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); } /* ---------------------------------------------------------------------- @@ -577,10 +577,9 @@ void ComputeBornMatrix::restore_atoms(int nall, int idir) // reset only idir coord - int k = dirlist[idir][0]; + int k = dirlist[idir][0]; double **x = atom->x; - for (int i = 0; i < nall; i++) - x[i][k] = temp_x[i][k]; + for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k]; } /* ---------------------------------------------------------------------- @@ -595,7 +594,7 @@ void ComputeBornMatrix::update_virial() // this may not be completely general // but it works for lj/cut and sw pair styles // and compute stress/atom output is unaffected - + int vflag = VIRIAL_PAIR; if (force->pair) force->pair->compute(eflag, vflag); @@ -625,32 +624,32 @@ void ComputeBornMatrix::virial_addon() int kd, nd, id, jd; int m; - double* sigv = compute_virial->vector; + double *sigv = compute_virial->vector; double modefactor[6] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5}; // you can compute these factor by looking at // every Dijkl terms and adding the proper virials // Take into account the symmetries. For example: // B2323 = s33+D2323; B3232= s22+D3232; - // but D3232=D2323 + // but D3232=D2323 // and Cijkl = (Bijkl+Bjikl+Bijlk+Bjilk)/4. = (Bijkl+Bjilk)/2. - - // these values have been verified correct to about 1e-9 - // against the analytic expressions for lj/cut - values_global[0] += 2.0*sigv[0]; - values_global[1] += 2.0*sigv[1]; - values_global[2] += 2.0*sigv[2]; - values_global[3] += sigv[2]; - values_global[4] += sigv[2]; - values_global[5] += sigv[1]; - values_global[6] += 0.0; - values_global[7] += 0.0; - values_global[8] += 0.0; - values_global[9] += 2.0*sigv[4]; - values_global[10] += 2.0*sigv[3]; + // these values have been verified correct to about 1e-9 + // against the analytic expressions for lj/cut + + values_global[0] += 2.0 * sigv[0]; + values_global[1] += 2.0 * sigv[1]; + values_global[2] += 2.0 * sigv[2]; + values_global[3] += sigv[2]; + values_global[4] += sigv[2]; + values_global[5] += sigv[1]; + values_global[6] += 0.0; + values_global[7] += 0.0; + values_global[8] += 0.0; + values_global[9] += 2.0 * sigv[4]; + values_global[10] += 2.0 * sigv[3]; values_global[11] += 0.0; - values_global[12] += 2.0*sigv[5]; + values_global[12] += 2.0 * sigv[5]; values_global[13] += 0.0; values_global[14] += 0.0; values_global[15] += 0.0; @@ -659,7 +658,6 @@ void ComputeBornMatrix::virial_addon() values_global[18] += 0.0; values_global[19] += 0.0; values_global[20] += sigv[5]; - } /* ---------------------------------------------------------------------- @@ -708,9 +706,9 @@ double ComputeBornMatrix::memory_usage() void ComputeBornMatrix::compute_bonds() { - int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar; + int i, m, n, nb, atom1, atom2, imol, iatom, btype, ivar; tagint tagprev; - double dx,dy,dz,rsq; + double dx, dy, dz, rsq; double **x = atom->x; double **v = atom->v; @@ -731,7 +729,7 @@ void ComputeBornMatrix::compute_bonds() Bond *bond = force->bond; - int a,b,c,d; + int a, b, c, d; double rij[3]; double rinv, r2inv; double pair_pref, dupair, du2pair; @@ -739,12 +737,13 @@ void ComputeBornMatrix::compute_bonds() // loop over all atoms and their bonds m = 0; - while (mbond_type[iatom][i]; - atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev); + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev); } if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; @@ -769,30 +768,25 @@ void ComputeBornMatrix::compute_bonds() dx = x[atom2][0] - x[atom1][0]; dy = x[atom2][1] - x[atom1][1]; dz = x[atom2][2] - x[atom1][2]; - domain->minimum_image(dx,dy,dz); - rsq = dx*dx + dy*dy + dz*dz; + domain->minimum_image(dx, dy, dz); + rsq = dx * dx + dy * dy + dz * dz; rij[0] = dx; rij[1] = dy; rij[2] = dz; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - pair_pref = 0.0; - dupair = 0.0; - du2pair = 0.0; - bond->born_matrix(btype,rsq,atom1,atom2,dupair,du2pair); - pair_pref = du2pair - dupair*rinv; + pair_pref = dupair = du2pair = 0.0; + bond->born_matrix(btype, rsq, atom1, atom2, dupair, du2pair); + pair_pref = du2pair - dupair * rinv; - a = 0; - b = 0; - c = 0; - d = 0; - for (i = 0; i<21; i++) { + a = b = c = d = 0; + for (i = 0; i < 21; i++) { a = albemunu[i][0]; b = albemunu[i][1]; c = albemunu[i][2]; d = albemunu[i][3]; - values_local[m+i] += pair_pref*rij[a]*rij[b]*rij[c]*rij[d]*r2inv; + values_local[m + i] += pair_pref * rij[a] * rij[b] * rij[c] * rij[d] * r2inv; } } } @@ -811,12 +805,12 @@ void ComputeBornMatrix::compute_bonds() void ComputeBornMatrix::compute_angles() { - int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar; + int i, m, n, na, atom1, atom2, atom3, imol, iatom, atype, ivar; tagint tagprev; - double delx1,dely1,delz1,delx2,dely2,delz2; - double rsq1,rsq2,r1,r2,cost; + double delx1, dely1, delz1, delx2, dely2, delz2; + double rsq1, rsq2, r1, r2, cost; double r1r2, r1r2inv; - double rsq1inv,rsq2inv,r1inv,r2inv,cinv; + double rsq1inv, rsq2inv, r1inv, r2inv, cinv; double *ptr; double **x = atom->x; @@ -839,7 +833,7 @@ void ComputeBornMatrix::compute_angles() Angle *angle = force->angle; - int a,b,c,d,e,f; + int a, b, c, d, e, f; double duang, du2ang; double del1[3], del2[3]; double dcos[6]; @@ -848,20 +842,16 @@ void ComputeBornMatrix::compute_angles() // Initializing array for intermediate cos derivatives // w regard to strain - for (i = 0; i < 6; i++) { - dcos[i] = 0; - } - for (i = 0; i < 21; i++) { - d2cos[i] = 0; - d2lncos[i] = 0; - } + for (i = 0; i < 6; i++) dcos[i] = 0; + for (i = 0; i < 21; i++) d2cos[i] = d2lncos[i] = 0; m = 0; while (m < nvalues) { for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == 1) na = num_angle[atom2]; + if (molecular == 1) + na = num_angle[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -879,8 +869,8 @@ void ComputeBornMatrix::compute_angles() if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; atype = onemols[imol]->angle_type[atom2][i]; tagprev = tag[atom2] - iatom - 1; - atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev); - atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev); + atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev); } if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; @@ -890,53 +880,51 @@ void ComputeBornMatrix::compute_angles() delx1 = x[atom1][0] - x[atom2][0]; dely1 = x[atom1][1] - x[atom2][1]; delz1 = x[atom1][2] - x[atom2][2]; - domain->minimum_image(delx1,dely1,delz1); + domain->minimum_image(delx1, dely1, delz1); del1[0] = delx1; del1[1] = dely1; del1[2] = delz1; - rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; - rsq1inv = 1.0/rsq1; + rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1; + rsq1inv = 1.0 / rsq1; r1 = sqrt(rsq1); - r1inv = 1.0/r1; + r1inv = 1.0 / r1; delx2 = x[atom3][0] - x[atom2][0]; dely2 = x[atom3][1] - x[atom2][1]; delz2 = x[atom3][2] - x[atom2][2]; - domain->minimum_image(delx2,dely2,delz2); + domain->minimum_image(delx2, dely2, delz2); del2[0] = delx2; del2[1] = dely2; del2[2] = delz2; - rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; - rsq2inv = 1.0/rsq2; + rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2; + rsq2inv = 1.0 / rsq2; r2 = sqrt(rsq2); - r2inv = 1.0/r2; + r2inv = 1.0 / r2; - r1r2 = delx1*delx2 + dely1*dely2 + delz1*delz2; - r1r2inv = 1/r1r2; + r1r2 = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + r1r2inv = 1 / r1r2; // cost = cosine of angle - cost = delx1*delx2 + dely1*dely2 + delz1*delz2; - cost /= r1*r2; + cost = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + cost /= r1 * r2; if (cost > 1.0) cost = 1.0; if (cost < -1.0) cost = -1.0; - cinv = 1.0/cost; + cinv = 1.0 / cost; // The method must return derivative with regards // to cos(theta)! // Use the chain rule if needed: // dU(t)/de = dt/dcos(t)*dU(t)/dt*dcos(t)/de // with dt/dcos(t) = -1/sin(t) - angle->born_matrix(atype,atom1,atom2,atom3,duang,du2ang); + angle->born_matrix(atype, atom1, atom2, atom3, duang, du2ang); // Voigt notation // 1 = 11, 2 = 22, 3 = 33 // 4 = 23, 5 = 13, 6 = 12 - a = 0; - b = 0; - c = 0; - d = 0; + a = b = c = d = 0; + // clang-format off for (i = 0; i<6; i++) { a = sigma_albe[i][0]; b = sigma_albe[i][1]; @@ -958,9 +946,10 @@ void ComputeBornMatrix::compute_angles() d2cos[i] = cost*d2lncos[i] + dcos[e]*dcos[f]*cinv; values_local[m+i] += duang*d2cos[i] + du2ang*dcos[e]*dcos[f]; } + // clang-format on } } - m+=21; + m += 21; } } @@ -974,11 +963,11 @@ void ComputeBornMatrix::compute_angles() void ComputeBornMatrix::compute_dihedrals() { - int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,dtype,ivar; + int i, m, n, nd, atom1, atom2, atom3, atom4, imol, iatom, dtype, ivar; tagint tagprev; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv; - double si,co,phi; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, ra2inv, rb2inv, rabinv; + double si, co, phi; double *ptr; double **x = atom->x; @@ -1001,8 +990,8 @@ void ComputeBornMatrix::compute_dihedrals() Dihedral *dihedral = force->dihedral; - double dudih,du2dih; - int a,b,c,d,e,f; + double dudih, du2dih; + int a, b, c, d, e, f; double b1sq; double b2sq; double b3sq; @@ -1025,25 +1014,17 @@ void ComputeBornMatrix::compute_dihedrals() double dcos[6]; double d2cos[21]; - for (i = 0; i < 6; i++) { - dmn[i] =0; - dmm[i] = 0; - dnn[i] = 0; - dcos[i] = 0; - } - for (i = 0; i < 21; i++) { - d2mn[i] = 0; - d2mm[i] = 0; - d2nn[i] = 0; - d2cos[i] = 0; - } + for (i = 0; i < 6; i++) dmn[i] = dmm[i] = dnn[i] = dcos[i] = 0; + + for (i = 0; i < 21; i++) d2mn[i] = d2mm[i] = d2nn[i] = d2cos[i] = 0; m = 0; while (m < nvalues) { for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == 1) nd = num_dihedral[atom2]; + if (molecular == 1) + nd = num_dihedral[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -1060,9 +1041,9 @@ void ComputeBornMatrix::compute_dihedrals() } else { if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; tagprev = tag[atom2] - iatom - 1; - atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev); - atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev); - atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev); + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev); } if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; @@ -1077,76 +1058,72 @@ void ComputeBornMatrix::compute_dihedrals() // dU(t)/de = dt/dcos(t)*dU(t)/dt*dcos(t)/de // with dt/dcos(t) = -1/sin(t) - dihedral->born_matrix(nd,atom1,atom2,atom3,atom4,dudih,du2dih); + dihedral->born_matrix(nd, atom1, atom2, atom3, atom4, dudih, du2dih); vb1x = x[atom1][0] - x[atom2][0]; vb1y = x[atom1][1] - x[atom2][1]; vb1z = x[atom1][2] - x[atom2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); + domain->minimum_image(vb1x, vb1y, vb1z); b1[0] = vb1x; b1[1] = vb1y; b1[2] = vb1z; - b1sq = b1[0]*b1[0]+b1[1]*b1[1]+b1[2]*b1[2]; + b1sq = b1[0] * b1[0] + b1[1] * b1[1] + b1[2] * b1[2]; vb2x = x[atom3][0] - x[atom2][0]; vb2y = x[atom3][1] - x[atom2][1]; vb2z = x[atom3][2] - x[atom2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); + domain->minimum_image(vb2x, vb2y, vb2z); b2[0] = vb2x; b2[1] = vb2y; b2[2] = vb2z; - b2sq = b2[0]*b2[0]+b2[1]*b2[1]+b2[2]*b2[2]; + b2sq = b2[0] * b2[0] + b2[1] * b2[1] + b2[2] * b2[2]; vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); + domain->minimum_image(vb2xm, vb2ym, vb2zm); vb3x = x[atom4][0] - x[atom3][0]; vb3y = x[atom4][1] - x[atom3][1]; vb3z = x[atom4][2] - x[atom3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); + domain->minimum_image(vb3x, vb3y, vb3z); b3[0] = vb3x; b3[1] = vb3y; b3[2] = vb3z; - b3sq = b3[0]*b3[0]+b3[1]*b3[1]+b3[2]*b3[2]; + b3sq = b3[0] * b3[0] + b3[1] * b3[1] + b3[2] * b3[2]; - b1b2 = b1[0]*b2[0]+b1[1]*b2[1]+b1[2]*b2[2]; - b1b3 = b1[0]*b3[0]+b1[1]*b3[1]+b1[2]*b3[2]; - b2b3 = b2[0]*b3[0]+b2[1]*b3[1]+b2[2]*b3[2]; + b1b2 = b1[0] * b2[0] + b1[1] * b2[1] + b1[2] * b2[2]; + b1b3 = b1[0] * b3[0] + b1[1] * b3[1] + b1[2] * b3[2]; + b2b3 = b2[0] * b3[0] + b2[1] * b3[1] + b2[2] * b3[2]; - ax = vb1y*vb2zm - vb1z*vb2ym; - ay = vb1z*vb2xm - vb1x*vb2zm; - az = vb1x*vb2ym - vb1y*vb2xm; - bx = vb3y*vb2zm - vb3z*vb2ym; - by = vb3z*vb2xm - vb3x*vb2zm; - bz = vb3x*vb2ym - vb3y*vb2xm; + ax = vb1y * vb2zm - vb1z * vb2ym; + ay = vb1z * vb2xm - vb1x * vb2zm; + az = vb1x * vb2ym - vb1y * vb2xm; + bx = vb3y * vb2zm - vb3z * vb2ym; + by = vb3z * vb2xm - vb3x * vb2zm; + bz = vb3x * vb2ym - vb3y * vb2xm; - rasq = ax*ax + ay*ay + az*az; - rbsq = bx*bx + by*by + bz*bz; - rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm; rg = sqrt(rgsq); ra2inv = rb2inv = 0.0; - if (rasq > 0) ra2inv = 1.0/rasq; - if (rbsq > 0) rb2inv = 1.0/rbsq; - rabinv = sqrt(ra2inv*rb2inv); + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + rabinv = sqrt(ra2inv * rb2inv); - co = (ax*bx + ay*by + az*bz)*rabinv; - si = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + co = (ax * bx + ay * by + az * bz) * rabinv; + si = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); if (co > 1.0) co = 1.0; if (co < -1.0) co = -1.0; - phi = atan2(si,co); + phi = atan2(si, co); // above a and b are m and n vectors // here they are integers indices - a = 0; - b = 0; - c = 0; - d = 0; - e = 0; - f = 0; + a = b = c = d = e = f = 0; + // clang-format off for (i = 0; i<6; i++) { a = sigma_albe[i][0]; b = sigma_albe[i][1]; @@ -1180,8 +1157,9 @@ void ComputeBornMatrix::compute_dihedrals() - rb2inv*d2nn[i]); values_local[m+i] += dudih*d2cos[i] + du2dih*dcos[e]*dcos[f]; } + // clang-format on } } - m+=21; + m += 21; } } diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.h b/src/EXTRA-COMPUTE/compute_born_matrix.h index 8e58b24ce2..71fe41bedf 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.h +++ b/src/EXTRA-COMPUTE/compute_born_matrix.h @@ -28,57 +28,56 @@ ComputeStyle(born/matrix,ComputeBornMatrix); namespace LAMMPS_NS { - class ComputeBornMatrix : public Compute { - public: - ComputeBornMatrix(class LAMMPS *, int, char **); - virtual ~ComputeBornMatrix() override; - void init() override; - void init_list(int, class NeighList *) override; - void compute_vector() override; - double memory_usage() override; +class ComputeBornMatrix : public Compute { + public: + ComputeBornMatrix(class LAMMPS *, int, char **); + virtual ~ComputeBornMatrix() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_vector() override; + double memory_usage() override; - private: + private: + // Born matrix contributions - // Born matrix contributions + void compute_pairs(); // pair and manybody + void compute_bonds(); // bonds + void compute_angles(); // angles + void compute_dihedrals(); // dihedrals + void compute_numdiff(); // stress virial finite differences + void displace_atoms(int, int, double); // displace atoms + void force_clear(int); // zero out force array + void update_virial(); // recalculate the virial + void restore_atoms(int, int); // restore atom positions + void virial_addon(); // restore atom positions + void reallocate(); // grow the atom arrays - void compute_pairs(); // pair and manybody - void compute_bonds(); // bonds - void compute_angles(); // angles - void compute_dihedrals(); // dihedrals - void compute_numdiff(); // stress virial finite differences - void displace_atoms(int, int, double); // displace atoms - void force_clear(int); // zero out force array - void update_virial(); // recalculate the virial - void restore_atoms(int, int); // restore atom positions - void virial_addon(); // restore atom positions - void reallocate(); // grow the atom arrays + int me; // process rank + int nvalues; // length of elastic tensor + int numflag; // 1 if using finite differences + double numdelta; // size of finite strain + int maxatom; // allocated size of atom arrays - int me; // process rank - int nvalues; // length of elastic tensor - int numflag; // 1 if using finite differences - double numdelta; // size of finite strain - int maxatom; // allocated size of atom arrays + int pairflag, bondflag, angleflag; + int dihedflag, impflag, kspaceflag; - int pairflag, bondflag, angleflag; - int dihedflag, impflag, kspaceflag; + double *values_local, *values_global; + double pos, pos1, dt, nktv2p, ftm2v; + class NeighList *list; - double *values_local,*values_global; - double pos,pos1,dt,nktv2p,ftm2v; - class NeighList *list; + char *id_virial; // name of virial compute + class Compute *compute_virial; // pointer to virial compute - char *id_virial; // name of virial compute - class Compute *compute_virial; // pointer to virial compute - - static constexpr int NDIR_VIRIAL = 6; // dimension of virial and strain vectors - static constexpr int NXYZ_VIRIAL = 3; // number of Cartesian coordinates - int revalbe[NDIR_VIRIAL][NDIR_VIRIAL]; - int virialVtoV[NDIR_VIRIAL]; - double **temp_x; // original coords - double **temp_f; // original forces - double fixedpoint[NXYZ_VIRIAL]; // displacement field origin - int dirlist[NDIR_VIRIAL][2]; // strain cartesian indices - }; -} + static constexpr int NDIR_VIRIAL = 6; // dimension of virial and strain vectors + static constexpr int NXYZ_VIRIAL = 3; // number of Cartesian coordinates + int revalbe[NDIR_VIRIAL][NDIR_VIRIAL]; + int virialVtoV[NDIR_VIRIAL]; + double **temp_x; // original coords + double **temp_f; // original forces + double fixedpoint[NXYZ_VIRIAL]; // displacement field origin + int dirlist[NDIR_VIRIAL][2]; // strain cartesian indices +}; +} // namespace LAMMPS_NS #endif #endif