diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst index 44210218c0..60eded42e0 100644 --- a/doc/src/Developer_utils.rst +++ b/doc/src/Developer_utils.rst @@ -334,10 +334,11 @@ arguments of commands in LAMMPS are parsed and to make abstractions of repetitive tasks. The :cpp:class:`LAMMPS_NS::ArgInfo` class provides an abstraction -for parsing references to compute or fix styles or variables. These -would start with a "c\_", "f\_", "v\_" followed by the ID or name of -than instance and may be postfixed with one or two array indices -"[]" with numbers > 0. +for parsing references to compute or fix styles, variables or custom +integer or double properties handled by :doc:`fix property/atom `. +These would start with a "c\_", "f\_", "v\_", "d\_", "d2\_", "i\_", or "i2_" +followed by the ID or name of than instance and may be postfixed with +one or two array indices "[]" with numbers > 0. A typical code segment would look like this: diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp index 5c0602a453..23ec79fa98 100644 --- a/src/ASPHERE/pair_resquared.cpp +++ b/src/ASPHERE/pair_resquared.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -18,32 +17,31 @@ #include "pair_resquared.h" -#include -#include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp), - cr60(pow(60.0,1.0/3.0)), - b_alpha(45.0/56.0) +PairRESquared::PairRESquared(LAMMPS *lmp) : + Pair(lmp), cr60(pow(60.0, 1.0 / 3.0)), b_alpha(45.0 / 56.0) { single_enable = 0; - cr60 = pow(60.0,1.0/3.0); - b_alpha = 45.0/56.0; - solv_f_a = 3.0/(16.0*atan(1.0)*-36.0); - solv_f_r = 3.0/(16.0*atan(1.0)*2025.0); + cr60 = pow(60.0, 1.0 / 3.0); + b_alpha = 45.0 / 56.0; + solv_f_a = 3.0 / (16.0 * atan(1.0) * -36.0); + solv_f_r = 3.0 / (16.0 * atan(1.0) * 2025.0); } /* ---------------------------------------------------------------------- @@ -68,8 +66,8 @@ PairRESquared::~PairRESquared() memory->destroy(lj3); memory->destroy(lj4); memory->destroy(offset); - delete [] lshape; - delete [] setwell; + delete[] lshape; + delete[] setwell; } } @@ -77,14 +75,14 @@ PairRESquared::~PairRESquared() void PairRESquared::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double evdwl,one_eng,rsq,r2inv,r6inv,forcelj,factor_lj; - double fforce[3],ttor[3],rtor[3],r12[3]; - int *ilist,*jlist,*numneigh,**firstneigh; - RE2Vars wi,wj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double evdwl, one_eng, rsq, r2inv, r6inv, forcelj, factor_lj; + double fforce[3], ttor[3], rtor[3], r12[3]; + int *ilist, *jlist, *numneigh, **firstneigh; + RE2Vars wi, wj; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -107,7 +105,7 @@ void PairRESquared::compute(int eflag, int vflag) // not a LJ sphere - if (lshape[itype] != 0.0) precompute_i(i,wi); + if (lshape[itype] != 0.0) precompute_i(i, wi); jlist = firstneigh[i]; jnum = numneigh[i]; @@ -119,10 +117,10 @@ void PairRESquared::compute(int eflag, int vflag) // r12 = center to center vector - r12[0] = x[j][0]-x[i][0]; - r12[1] = x[j][1]-x[i][1]; - r12[2] = x[j][2]-x[i][2]; - rsq = MathExtra::dot3(r12,r12); + r12[0] = x[j][0] - x[i][0]; + r12[1] = x[j][1] - x[i][1]; + r12[2] = x[j][2] - x[i][2]; + rsq = MathExtra::dot3(r12, r12); jtype = type[j]; // compute if less than cutoff @@ -132,49 +130,50 @@ void PairRESquared::compute(int eflag, int vflag) switch (form[itype][jtype]) { - case SPHERE_SPHERE: - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); - forcelj *= -r2inv; - if (eflag) one_eng = - r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) - - offset[itype][jtype]; - fforce[0] = r12[0]*forcelj; - fforce[1] = r12[1]*forcelj; - fforce[2] = r12[2]*forcelj; - break; + case SPHERE_SPHERE: + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]); + forcelj *= -r2inv; + if (eflag) { + one_eng = r6inv * (r6inv * lj3[itype][jtype] - lj4[itype][jtype]); + one_eng -= offset[itype][jtype]; + } + fforce[0] = r12[0] * forcelj; + fforce[1] = r12[1] * forcelj; + fforce[2] = r12[2] * forcelj; + break; - case SPHERE_ELLIPSE: - precompute_i(j,wj); - if (newton_pair || j < nlocal) { - one_eng = resquared_lj(j,i,wj,r12,rsq,fforce,rtor,true); - tor[j][0] += rtor[0]*factor_lj; - tor[j][1] += rtor[1]*factor_lj; - tor[j][2] += rtor[2]*factor_lj; - } else - one_eng = resquared_lj(j,i,wj,r12,rsq,fforce,rtor,false); - break; + case SPHERE_ELLIPSE: + precompute_i(j, wj); + if (newton_pair || j < nlocal) { + one_eng = resquared_lj(j, i, wj, r12, rsq, fforce, rtor, true); + tor[j][0] += rtor[0] * factor_lj; + tor[j][1] += rtor[1] * factor_lj; + tor[j][2] += rtor[2] * factor_lj; + } else + one_eng = resquared_lj(j, i, wj, r12, rsq, fforce, rtor, false); + break; - case ELLIPSE_SPHERE: - one_eng = resquared_lj(i,j,wi,r12,rsq,fforce,ttor,true); - tor[i][0] += ttor[0]*factor_lj; - tor[i][1] += ttor[1]*factor_lj; - tor[i][2] += ttor[2]*factor_lj; - break; + case ELLIPSE_SPHERE: + one_eng = resquared_lj(i, j, wi, r12, rsq, fforce, ttor, true); + tor[i][0] += ttor[0] * factor_lj; + tor[i][1] += ttor[1] * factor_lj; + tor[i][2] += ttor[2] * factor_lj; + break; - default: - precompute_i(j,wj); - one_eng = resquared_analytic(i,j,wi,wj,r12,rsq,fforce,ttor,rtor); - tor[i][0] += ttor[0]*factor_lj; - tor[i][1] += ttor[1]*factor_lj; - tor[i][2] += ttor[2]*factor_lj; - if (newton_pair || j < nlocal) { - tor[j][0] += rtor[0]*factor_lj; - tor[j][1] += rtor[1]*factor_lj; - tor[j][2] += rtor[2]*factor_lj; - } - break; + default: + precompute_i(j, wj); + one_eng = resquared_analytic(i, j, wi, wj, r12, rsq, fforce, ttor, rtor); + tor[i][0] += ttor[0] * factor_lj; + tor[i][1] += ttor[1] * factor_lj; + tor[i][2] += ttor[2] * factor_lj; + if (newton_pair || j < nlocal) { + tor[j][0] += rtor[0] * factor_lj; + tor[j][1] += rtor[1] * factor_lj; + tor[j][2] += rtor[2] * factor_lj; + } + break; } fforce[0] *= factor_lj; @@ -190,11 +189,11 @@ void PairRESquared::compute(int eflag, int vflag) f[j][2] -= fforce[2]; } - if (eflag) evdwl = factor_lj*one_eng; + if (eflag) evdwl = factor_lj * one_eng; - if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, - evdwl,0.0,fforce[0],fforce[1],fforce[2], - -r12[0],-r12[1],-r12[2]); + if (evflag) + ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fforce[0], fforce[1], fforce[2], + -r12[0], -r12[1], -r12[2]); } } } @@ -209,31 +208,30 @@ void PairRESquared::compute(int eflag, int vflag) void PairRESquared::allocate() { allocated = 1; - int n = atom->ntypes; + const int n = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, n, n, "pair:setflag"); + for (int i = 1; i < n; i++) + for (int j = i; j < n; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, n, n, "pair:cutsq"); - memory->create(form,n+1,n+1,"pair:form"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(shape1,n+1,3,"pair:shape1"); - memory->create(shape2,n+1,3,"pair:shape2"); - memory->create(well,n+1,3,"pair:well"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(form, n, n, "pair:form"); + memory->create(epsilon, n, n, "pair:epsilon"); + memory->create(sigma, n, n, "pair:sigma"); + memory->create(shape1, n, 3, "pair:shape1"); + memory->create(shape2, n, 3, "pair:shape2"); + memory->create(well, n, 3, "pair:well"); + memory->create(cut, n, n, "pair:cut"); + memory->create(lj1, n, n, "pair:lj1"); + memory->create(lj2, n, n, "pair:lj2"); + memory->create(lj3, n, n, "pair:lj3"); + memory->create(lj4, n, n, "pair:lj4"); + memory->create(offset, n, n, "pair:offset"); - lshape = new double[n+1]; - setwell = new int[n+1]; - for (int i = 1; i <= n; i++) setwell[i] = 0; + lshape = new double[n]; + setwell = new int[n]; + for (int i = 1; i < n; i++) setwell[i] = 0; } /* ---------------------------------------------------------------------- @@ -242,14 +240,14 @@ void PairRESquared::allocate() void PairRESquared::settings(int narg, char **arg) { - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); - cut_global = utils::numeric(FLERR,arg[0],false,lmp); + cut_global = utils::numeric(FLERR, arg[0], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; @@ -262,29 +260,28 @@ void PairRESquared::settings(int narg, char **arg) void PairRESquared::coeff(int narg, char **arg) { - if (narg < 10 || narg > 11) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 10 || narg > 11) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); - double eia_one = utils::numeric(FLERR,arg[4],false,lmp); - double eib_one = utils::numeric(FLERR,arg[5],false,lmp); - double eic_one = utils::numeric(FLERR,arg[6],false,lmp); - double eja_one = utils::numeric(FLERR,arg[7],false,lmp); - double ejb_one = utils::numeric(FLERR,arg[8],false,lmp); - double ejc_one = utils::numeric(FLERR,arg[9],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); + double eia_one = utils::numeric(FLERR, arg[4], false, lmp); + double eib_one = utils::numeric(FLERR, arg[5], false, lmp); + double eic_one = utils::numeric(FLERR, arg[6], false, lmp); + double eja_one = utils::numeric(FLERR, arg[7], false, lmp); + double ejb_one = utils::numeric(FLERR, arg[8], false, lmp); + double ejc_one = utils::numeric(FLERR, arg[9], false, lmp); double cut_one = cut_global; - if (narg == 11) cut_one = utils::numeric(FLERR,arg[10],false,lmp); + if (narg == 11) cut_one = utils::numeric(FLERR, arg[10], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; @@ -292,22 +289,26 @@ void PairRESquared::coeff(int narg, char **arg) well[i][0] = eia_one; well[i][1] = eib_one; well[i][2] = eic_one; - if (eia_one == 1.0 && eib_one == 1.0 && eic_one == 1.0) setwell[i] = 2; - else setwell[i] = 1; + if (eia_one == 1.0 && eib_one == 1.0 && eic_one == 1.0) + setwell[i] = 2; + else + setwell[i] = 1; } if (eja_one != 0.0 || ejb_one != 0.0 || ejc_one != 0.0) { well[j][0] = eja_one; well[j][1] = ejb_one; well[j][2] = ejc_one; - if (eja_one == 1.0 && ejb_one == 1.0 && ejc_one == 1.0) setwell[j] = 2; - else setwell[j] = 1; + if (eja_one == 1.0 && ejb_one == 1.0 && ejc_one == 1.0) + setwell[j] = 2; + else + setwell[j] = 1; } setflag[i][j] = 1; count++; } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -317,21 +318,21 @@ void PairRESquared::coeff(int narg, char **arg) void PairRESquared::init_style() { avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); - if (!avec) error->all(FLERR,"Pair resquared requires atom style ellipsoid"); + if (!avec) error->all(FLERR, "Pair resquared requires atom style ellipsoid"); - neighbor->request(this,instance_me); + neighbor->request(this, instance_me); // per-type shape precalculations // require that atom shapes are identical within each type for (int i = 1; i <= atom->ntypes; i++) { - if (!atom->shape_consistency(i,shape1[i][0],shape1[i][1],shape1[i][2])) - error->all(FLERR,"Pair resquared requires atoms with same type have same shape"); + if (!atom->shape_consistency(i, shape1[i][0], shape1[i][1], shape1[i][2])) + error->all(FLERR, "Pair resquared requires atoms with same type have same shape"); if (setwell[i]) { - shape2[i][0] = shape1[i][0]*shape1[i][0]; - shape2[i][1] = shape1[i][1]*shape1[i][1]; - shape2[i][2] = shape1[i][2]*shape1[i][2]; - lshape[i] = shape1[i][0]*shape1[i][1]*shape1[i][2]; + shape2[i][0] = shape1[i][0] * shape1[i][0]; + shape2[i][1] = shape1[i][1] * shape1[i][1]; + shape2[i][2] = shape1[i][2] * shape1[i][2]; + lshape[i] = shape1[i][0] * shape1[i][1] * shape1[i][2]; } } } @@ -343,14 +344,12 @@ void PairRESquared::init_style() double PairRESquared::init_one(int i, int j) { if (setwell[i] == 0 || setwell[j] == 0) - error->all(FLERR,"Pair resquared epsilon a,b,c coeffs are not all set"); + error->all(FLERR, "Pair resquared epsilon a,b,c coeffs are not all set"); int ishape = 0; - if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0) - ishape = 1; + if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0) ishape = 1; int jshape = 0; - if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0) - jshape = 1; + if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0) jshape = 1; if (ishape == 0 && jshape == 0) { form[i][j] = SPHERE_SPHERE; @@ -371,28 +370,27 @@ double PairRESquared::init_one(int i, int j) if (setflag[i][j] == 0) { if (setflag[j][i] == 0) { if (ishape == 0 && jshape == 0) { - epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], - sigma[i][i],sigma[j][j]); - sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], sigma[i][i], sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i], sigma[j][j]); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); } else - error->all(FLERR, - "Pair resquared epsilon and sigma coeffs are not all set"); + error->all(FLERR, "Pair resquared epsilon and sigma coeffs are not all set"); } epsilon[i][j] = epsilon[j][i]; sigma[i][j] = sigma[j][i]; cut[i][j] = cut[j][i]; } - lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j], 12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; - offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio, 12.0) - pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; epsilon[j][i] = epsilon[i][j]; sigma[j][i] = sigma[i][j]; @@ -413,16 +411,16 @@ void PairRESquared::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) { - fwrite(&setwell[i],sizeof(int),1,fp); - if (setwell[i]) fwrite(&well[i][0],sizeof(double),3,fp); + fwrite(&setwell[i], sizeof(int), 1, fp); + if (setwell[i]) fwrite(&well[i][0], sizeof(double), 3, fp); for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); } } } @@ -437,27 +435,27 @@ void PairRESquared::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) { - if (me == 0) utils::sfread(FLERR,&setwell[i],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setwell[i],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setwell[i], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setwell[i], 1, MPI_INT, 0, world); if (setwell[i]) { - if (me == 0) utils::sfread(FLERR,&well[i][0],sizeof(double),3,fp,nullptr,error); - MPI_Bcast(&well[i][0],3,MPI_DOUBLE,0,world); + if (me == 0) utils::sfread(FLERR, &well[i][0], sizeof(double), 3, fp, nullptr, error); + MPI_Bcast(&well[i][0], 3, MPI_DOUBLE, 0, world); } for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -469,8 +467,8 @@ void PairRESquared::read_restart(FILE *fp) void PairRESquared::write_restart_settings(FILE *fp) { - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -481,34 +479,34 @@ void PairRESquared::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- Precompute per-particle temporaries for RE-squared calculation ------------------------------------------------------------------------- */ -void PairRESquared::precompute_i(const int i,RE2Vars &ws) +void PairRESquared::precompute_i(const int i, RE2Vars &ws) { - double aTs[3][3]; // A1'*S1^2 + double aTs[3][3]; // A1'*S1^2 int *ellipsoid = atom->ellipsoid; AtomVecEllipsoid::Bonus *bonus = avec->bonus; - MathExtra::quat_to_mat_trans(bonus[ellipsoid[i]].quat,ws.A); - MathExtra::transpose_diag3(ws.A,well[atom->type[i]],ws.aTe); - MathExtra::transpose_diag3(ws.A,shape2[atom->type[i]],aTs); - MathExtra::diag_times3(shape2[atom->type[i]],ws.A,ws.sa); - MathExtra::times3(aTs,ws.A,ws.gamma); - MathExtra::rotation_generator_x(ws.A,ws.lA[0]); - MathExtra::rotation_generator_y(ws.A,ws.lA[1]); - MathExtra::rotation_generator_z(ws.A,ws.lA[2]); - for (int i=0; i<3; i++) { - MathExtra::times3(aTs,ws.lA[i],ws.lAtwo[i]); - MathExtra::transpose_times3(ws.lA[i],ws.sa,ws.lAsa[i]); - MathExtra::plus3(ws.lAsa[i],ws.lAtwo[i],ws.lAsa[i]); + MathExtra::quat_to_mat_trans(bonus[ellipsoid[i]].quat, ws.A); + MathExtra::transpose_diag3(ws.A, well[atom->type[i]], ws.aTe); + MathExtra::transpose_diag3(ws.A, shape2[atom->type[i]], aTs); + MathExtra::diag_times3(shape2[atom->type[i]], ws.A, ws.sa); + MathExtra::times3(aTs, ws.A, ws.gamma); + MathExtra::rotation_generator_x(ws.A, ws.lA[0]); + MathExtra::rotation_generator_y(ws.A, ws.lA[1]); + MathExtra::rotation_generator_z(ws.A, ws.lA[2]); + for (int m = 0; m < 3; m++) { + MathExtra::times3(aTs, ws.lA[m], ws.lAtwo[m]); + MathExtra::transpose_times3(ws.lA[m], ws.sa, ws.lAsa[m]); + MathExtra::plus3(ws.lAsa[m], ws.lAtwo[m], ws.lAsa[m]); } } @@ -517,6 +515,7 @@ void PairRESquared::precompute_i(const int i,RE2Vars &ws) derivative of m (m2) ------------------------------------------------------------------------- */ +// clang-format off double PairRESquared::det_prime(const double m[3][3], const double m2[3][3]) { double ans; @@ -531,458 +530,441 @@ double PairRESquared::det_prime(const double m[3][3], const double m2[3][3]) m2[2][0]*m[0][1]*m[1][2] - m2[2][0]*m[0][2]*m[1][1]; return ans; } +// clang-format on /* ---------------------------------------------------------------------- Compute the energy, force, torque for a pair (INTEGRATED-INTEGRATED) ------------------------------------------------------------------------- */ -double PairRESquared::resquared_analytic(const int i, const int j, - const RE2Vars &wi, const RE2Vars &wj, - const double *r, const double rsq, - double *fforce, double *ttor, - double *rtor) +double PairRESquared::resquared_analytic(const int i, const int j, const RE2Vars &wi, + const RE2Vars &wj, const double *r, const double rsq, + double *fforce, double *ttor, double *rtor) { int *type = atom->type; // pair computations for energy, force, torque - double z1[3],z2[3]; // A1*rhat # don't need to store - double v1[3],v2[3]; // inv(S1^2)*z1 # don't need to store - double sigma1,sigma2; // 1/sqrt(z1'*v1) - double sigma1p2,sigma2p2; // sigma1^2 - double rnorm; // L2 norm of r - double rhat[3]; // r/rnorm - double s[3]; // inv(gamma1+gamma2)*rhat - double sigma12; // 1/sqrt(0.5*s'*rhat) - double H12[3][3]; // gamma1/sigma1+gamma2/sigma2 - double dH; // det(H12) - double lambda; // dS1/sigma1p2+dS2/sigma2p2 - double nu; // sqrt(dH/(sigma1+sigma2)) - double w[3]; // inv(A1'*E1*A1+A2'*E2*A2)*rhat - double h12; // rnorm-sigma12; - double eta; // lambda/nu - double chi; // 2*rhat'*w - double sprod; // dS1*dS2 - double sigh; // sigma/h12 - double tprod; // eta*chi*sigh - double Ua,Ur; // attractive/repulsive parts of potential + double z1[3], z2[3]; // A1*rhat # don't need to store + double v1[3], v2[3]; // inv(S1^2)*z1 # don't need to store + double sigma1, sigma2; // 1/sqrt(z1'*v1) + double sigma1p2, sigma2p2; // sigma1^2 + double rnorm; // L2 norm of r + double rhat[3]; // r/rnorm + double s[3]; // inv(gamma1+gamma2)*rhat + double sigma12; // 1/sqrt(0.5*s'*rhat) + double H12[3][3]; // gamma1/sigma1+gamma2/sigma2 + double dH; // det(H12) + double lambda; // dS1/sigma1p2+dS2/sigma2p2 + double nu; // sqrt(dH/(sigma1+sigma2)) + double w[3]; // inv(A1'*E1*A1+A2'*E2*A2)*rhat + double h12; // rnorm-sigma12; + double eta; // lambda/nu + double chi; // 2*rhat'*w + double sprod; // dS1*dS2 + double sigh; // sigma/h12 + double tprod; // eta*chi*sigh + double Ua, Ur; // attractive/repulsive parts of potential // pair computations for force, torque - double sec; // sigma*eta*chi - double sigma1p3, sigma2p3; // sigma1^3 - double vsigma1[3], vsigma2[3]; // sigma1^3*v1; - double sigma12p3; // sigma12^3 - double gsigma1[3][3], gsigma2[3][3]; // -gamma1/sigma1^2 - double tsig1sig2; // eta/(2*(sigma1+sigma2)) - double tdH; // eta/(2*dH) - double teta1,teta2; // 2*eta/lambda*dS1/sigma1p3 - double fourw[3]; // 4*w; - double spr[3]; // 0.5*sigma12^3*s - double hsec; // h12+[3,b_alpha]*sec - double dspu; // 1/h12 - 1/hsec + temp - double pbsu; // 3*sigma/hsec - double dspr; // 7/h12-1/hsec+temp - double pbsr; // b_alpha*sigma/hsec; - double u[3]; // (-rhat(i)*rhat+eye(:,i))/rnorm - double u1[3],u2[3]; // A1*u - double dsigma1,dsigma2; // u1'*vsigma1 (force) p'*vsigma1 (tor) - double dH12[3][3]; // dsigma1*gsigma1 + dsigma2*gsigma2 - double ddH; // derivative of det(H12) - double deta,dchi,dh12; // derivatives of eta,chi,h12 - double dUr,dUa; // derivatives of Ua,Ur + double sec; // sigma*eta*chi + double sigma1p3, sigma2p3; // sigma1^3 + double vsigma1[3], vsigma2[3]; // sigma1^3*v1; + double sigma12p3; // sigma12^3 + double gsigma1[3][3], gsigma2[3][3]; // -gamma1/sigma1^2 + double tsig1sig2; // eta/(2*(sigma1+sigma2)) + double tdH; // eta/(2*dH) + double teta1, teta2; // 2*eta/lambda*dS1/sigma1p3 + double fourw[3]; // 4*w; + double spr[3]; // 0.5*sigma12^3*s + double hsec; // h12+[3,b_alpha]*sec + double dspu; // 1/h12 - 1/hsec + temp + double pbsu; // 3*sigma/hsec + double dspr; // 7/h12-1/hsec+temp + double pbsr; // b_alpha*sigma/hsec; + double u[3]; // (-rhat(i)*rhat+eye(:,i))/rnorm + double u1[3], u2[3]; // A1*u + double dsigma1, dsigma2; // u1'*vsigma1 (force) p'*vsigma1 (tor) + double dH12[3][3]; // dsigma1*gsigma1 + dsigma2*gsigma2 + double ddH; // derivative of det(H12) + double deta, dchi, dh12; // derivatives of eta,chi,h12 + double dUr, dUa; // derivatives of Ua,Ur // pair computations for torque - double fwae[3]; // -fourw'*aTe - double p[3]; // lA*rhat + double fwae[3]; // -fourw'*aTe + double p[3]; // lA*rhat rnorm = sqrt(rsq); - rhat[0] = r[0]/rnorm; - rhat[1] = r[1]/rnorm; - rhat[2] = r[2]/rnorm; + rhat[0] = r[0] / rnorm; + rhat[1] = r[1] / rnorm; + rhat[2] = r[2] / rnorm; // energy double temp[3][3]; - MathExtra::plus3(wi.gamma,wj.gamma,temp); - int ierror = MathExtra::mldivide3(temp,rhat,s); - if (ierror) error->all(FLERR,"Bad matrix inversion in mldivide3"); + MathExtra::plus3(wi.gamma, wj.gamma, temp); + int ierror = MathExtra::mldivide3(temp, rhat, s); + if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3"); - sigma12 = 1.0/sqrt(0.5*MathExtra::dot3(s,rhat)); - MathExtra::matvec(wi.A,rhat,z1); - MathExtra::matvec(wj.A,rhat,z2); - v1[0] = z1[0]/shape2[type[i]][0]; - v1[1] = z1[1]/shape2[type[i]][1]; - v1[2] = z1[2]/shape2[type[i]][2]; - v2[0] = z2[0]/shape2[type[j]][0]; - v2[1] = z2[1]/shape2[type[j]][1]; - v2[2] = z2[2]/shape2[type[j]][2]; - sigma1 = 1.0/sqrt(MathExtra::dot3(z1,v1)); - sigma2 = 1.0/sqrt(MathExtra::dot3(z2,v2)); - H12[0][0] = wi.gamma[0][0]/sigma1+wj.gamma[0][0]/sigma2; - H12[0][1] = wi.gamma[0][1]/sigma1+wj.gamma[0][1]/sigma2; - H12[0][2] = wi.gamma[0][2]/sigma1+wj.gamma[0][2]/sigma2; - H12[1][0] = wi.gamma[1][0]/sigma1+wj.gamma[1][0]/sigma2; - H12[1][1] = wi.gamma[1][1]/sigma1+wj.gamma[1][1]/sigma2; - H12[1][2] = wi.gamma[1][2]/sigma1+wj.gamma[1][2]/sigma2; - H12[2][0] = wi.gamma[2][0]/sigma1+wj.gamma[2][0]/sigma2; - H12[2][1] = wi.gamma[2][1]/sigma1+wj.gamma[2][1]/sigma2; - H12[2][2] = wi.gamma[2][2]/sigma1+wj.gamma[2][2]/sigma2; - dH=MathExtra::det3(H12); - sigma1p2 = sigma1*sigma1; - sigma2p2 = sigma2*sigma2; - lambda = lshape[type[i]]/sigma1p2 + lshape[type[j]]/sigma2p2; - nu = sqrt(dH/(sigma1+sigma2)); - MathExtra::times3(wi.aTe,wi.A,temp); + sigma12 = 1.0 / sqrt(0.5 * MathExtra::dot3(s, rhat)); + MathExtra::matvec(wi.A, rhat, z1); + MathExtra::matvec(wj.A, rhat, z2); + v1[0] = z1[0] / shape2[type[i]][0]; + v1[1] = z1[1] / shape2[type[i]][1]; + v1[2] = z1[2] / shape2[type[i]][2]; + v2[0] = z2[0] / shape2[type[j]][0]; + v2[1] = z2[1] / shape2[type[j]][1]; + v2[2] = z2[2] / shape2[type[j]][2]; + sigma1 = 1.0 / sqrt(MathExtra::dot3(z1, v1)); + sigma2 = 1.0 / sqrt(MathExtra::dot3(z2, v2)); + H12[0][0] = wi.gamma[0][0] / sigma1 + wj.gamma[0][0] / sigma2; + H12[0][1] = wi.gamma[0][1] / sigma1 + wj.gamma[0][1] / sigma2; + H12[0][2] = wi.gamma[0][2] / sigma1 + wj.gamma[0][2] / sigma2; + H12[1][0] = wi.gamma[1][0] / sigma1 + wj.gamma[1][0] / sigma2; + H12[1][1] = wi.gamma[1][1] / sigma1 + wj.gamma[1][1] / sigma2; + H12[1][2] = wi.gamma[1][2] / sigma1 + wj.gamma[1][2] / sigma2; + H12[2][0] = wi.gamma[2][0] / sigma1 + wj.gamma[2][0] / sigma2; + H12[2][1] = wi.gamma[2][1] / sigma1 + wj.gamma[2][1] / sigma2; + H12[2][2] = wi.gamma[2][2] / sigma1 + wj.gamma[2][2] / sigma2; + dH = MathExtra::det3(H12); + sigma1p2 = sigma1 * sigma1; + sigma2p2 = sigma2 * sigma2; + lambda = lshape[type[i]] / sigma1p2 + lshape[type[j]] / sigma2p2; + nu = sqrt(dH / (sigma1 + sigma2)); + MathExtra::times3(wi.aTe, wi.A, temp); double temp2[3][3]; - MathExtra::times3(wj.aTe,wj.A,temp2); - MathExtra::plus3(temp,temp2,temp); - ierror = MathExtra::mldivide3(temp,rhat,w); - if (ierror) error->all(FLERR,"Bad matrix inversion in mldivide3"); + MathExtra::times3(wj.aTe, wj.A, temp2); + MathExtra::plus3(temp, temp2, temp); + ierror = MathExtra::mldivide3(temp, rhat, w); + if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3"); - h12 = rnorm-sigma12; - eta = lambda/nu; - chi = 2.0*MathExtra::dot3(rhat,w); + h12 = rnorm - sigma12; + eta = lambda / nu; + chi = 2.0 * MathExtra::dot3(rhat, w); sprod = lshape[type[i]] * lshape[type[j]]; - sigh = sigma[type[i]][type[j]]/h12; - tprod = eta*chi*sigh; + sigh = sigma[type[i]][type[j]] / h12; + tprod = eta * chi * sigh; - double stemp = h12/2.0; - Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* - (shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)* - (shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp); - Ua = (1.0+3.0*tprod)*sprod/Ua; - Ua = epsilon[type[i]][type[j]]*Ua/-36.0; + double stemp = h12 / 2.0; + Ua = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) * + (shape1[type[j]][0] + stemp) * (shape1[type[j]][1] + stemp) * (shape1[type[j]][2] + stemp); + Ua = (1.0 + 3.0 * tprod) * sprod / Ua; + Ua = epsilon[type[i]][type[j]] * Ua / -36.0; - stemp = h12/cr60; - Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* - (shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)* - (shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp); - Ur = (1.0+b_alpha*tprod)*sprod/Ur; - Ur = epsilon[type[i]][type[j]]*Ur*pow(sigh,6.0)/2025.0; + stemp = h12 / cr60; + Ur = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) * + (shape1[type[j]][0] + stemp) * (shape1[type[j]][1] + stemp) * (shape1[type[j]][2] + stemp); + Ur = (1.0 + b_alpha * tprod) * sprod / Ur; + Ur = epsilon[type[i]][type[j]] * Ur * pow(sigh, 6.0) / 2025.0; // force - sec = sigma[type[i]][type[j]]*eta*chi; - sigma12p3 = pow(sigma12,3.0); - sigma1p3 = sigma1p2*sigma1; - sigma2p3 = sigma2p2*sigma2; - vsigma1[0] = -sigma1p3*v1[0]; - vsigma1[1] = -sigma1p3*v1[1]; - vsigma1[2] = -sigma1p3*v1[2]; - vsigma2[0] = -sigma2p3*v2[0]; - vsigma2[1] = -sigma2p3*v2[1]; - vsigma2[2] = -sigma2p3*v2[2]; - gsigma1[0][0] = -wi.gamma[0][0]/sigma1p2; - gsigma1[0][1] = -wi.gamma[0][1]/sigma1p2; - gsigma1[0][2] = -wi.gamma[0][2]/sigma1p2; - gsigma1[1][0] = -wi.gamma[1][0]/sigma1p2; - gsigma1[1][1] = -wi.gamma[1][1]/sigma1p2; - gsigma1[1][2] = -wi.gamma[1][2]/sigma1p2; - gsigma1[2][0] = -wi.gamma[2][0]/sigma1p2; - gsigma1[2][1] = -wi.gamma[2][1]/sigma1p2; - gsigma1[2][2] = -wi.gamma[2][2]/sigma1p2; - gsigma2[0][0] = -wj.gamma[0][0]/sigma2p2; - gsigma2[0][1] = -wj.gamma[0][1]/sigma2p2; - gsigma2[0][2] = -wj.gamma[0][2]/sigma2p2; - gsigma2[1][0] = -wj.gamma[1][0]/sigma2p2; - gsigma2[1][1] = -wj.gamma[1][1]/sigma2p2; - gsigma2[1][2] = -wj.gamma[1][2]/sigma2p2; - gsigma2[2][0] = -wj.gamma[2][0]/sigma2p2; - gsigma2[2][1] = -wj.gamma[2][1]/sigma2p2; - gsigma2[2][2] = -wj.gamma[2][2]/sigma2p2; - tsig1sig2 = eta/(2.0*(sigma1+sigma2)); - tdH = eta/(2.0*dH); - teta1 = 2.0*eta/lambda; - teta2 = teta1*lshape[type[j]]/sigma2p3; - teta1 = teta1*lshape[type[i]]/sigma1p3; - fourw[0] = 4.0*w[0]; - fourw[1] = 4.0*w[1]; - fourw[2] = 4.0*w[2]; - spr[0] = 0.5*sigma12p3*s[0]; - spr[1] = 0.5*sigma12p3*s[1]; - spr[2] = 0.5*sigma12p3*s[2]; + sec = sigma[type[i]][type[j]] * eta * chi; + sigma12p3 = pow(sigma12, 3.0); + sigma1p3 = sigma1p2 * sigma1; + sigma2p3 = sigma2p2 * sigma2; + vsigma1[0] = -sigma1p3 * v1[0]; + vsigma1[1] = -sigma1p3 * v1[1]; + vsigma1[2] = -sigma1p3 * v1[2]; + vsigma2[0] = -sigma2p3 * v2[0]; + vsigma2[1] = -sigma2p3 * v2[1]; + vsigma2[2] = -sigma2p3 * v2[2]; + gsigma1[0][0] = -wi.gamma[0][0] / sigma1p2; + gsigma1[0][1] = -wi.gamma[0][1] / sigma1p2; + gsigma1[0][2] = -wi.gamma[0][2] / sigma1p2; + gsigma1[1][0] = -wi.gamma[1][0] / sigma1p2; + gsigma1[1][1] = -wi.gamma[1][1] / sigma1p2; + gsigma1[1][2] = -wi.gamma[1][2] / sigma1p2; + gsigma1[2][0] = -wi.gamma[2][0] / sigma1p2; + gsigma1[2][1] = -wi.gamma[2][1] / sigma1p2; + gsigma1[2][2] = -wi.gamma[2][2] / sigma1p2; + gsigma2[0][0] = -wj.gamma[0][0] / sigma2p2; + gsigma2[0][1] = -wj.gamma[0][1] / sigma2p2; + gsigma2[0][2] = -wj.gamma[0][2] / sigma2p2; + gsigma2[1][0] = -wj.gamma[1][0] / sigma2p2; + gsigma2[1][1] = -wj.gamma[1][1] / sigma2p2; + gsigma2[1][2] = -wj.gamma[1][2] / sigma2p2; + gsigma2[2][0] = -wj.gamma[2][0] / sigma2p2; + gsigma2[2][1] = -wj.gamma[2][1] / sigma2p2; + gsigma2[2][2] = -wj.gamma[2][2] / sigma2p2; + tsig1sig2 = eta / (2.0 * (sigma1 + sigma2)); + tdH = eta / (2.0 * dH); + teta1 = 2.0 * eta / lambda; + teta2 = teta1 * lshape[type[j]] / sigma2p3; + teta1 = teta1 * lshape[type[i]] / sigma1p3; + fourw[0] = 4.0 * w[0]; + fourw[1] = 4.0 * w[1]; + fourw[2] = 4.0 * w[2]; + spr[0] = 0.5 * sigma12p3 * s[0]; + spr[1] = 0.5 * sigma12p3 * s[1]; + spr[2] = 0.5 * sigma12p3 * s[2]; - stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+ - 1.0/(shape1[type[i]][1]*2.0+h12)+ - 1.0/(shape1[type[i]][2]*2.0+h12)+ - 1.0/(shape1[type[j]][0]*2.0+h12)+ - 1.0/(shape1[type[j]][1]*2.0+h12)+ - 1.0/(shape1[type[j]][2]*2.0+h12); - hsec = h12+3.0*sec; - dspu = 1.0/h12-1.0/hsec+stemp; - pbsu = 3.0*sigma[type[i]][type[j]]/hsec; + stemp = 1.0 / (shape1[type[i]][0] * 2.0 + h12) + 1.0 / (shape1[type[i]][1] * 2.0 + h12) + + 1.0 / (shape1[type[i]][2] * 2.0 + h12) + 1.0 / (shape1[type[j]][0] * 2.0 + h12) + + 1.0 / (shape1[type[j]][1] * 2.0 + h12) + 1.0 / (shape1[type[j]][2] * 2.0 + h12); + hsec = h12 + 3.0 * sec; + dspu = 1.0 / h12 - 1.0 / hsec + stemp; + pbsu = 3.0 * sigma[type[i]][type[j]] / hsec; - stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+ - 1.0/(shape1[type[i]][1]*cr60+h12)+ - 1.0/(shape1[type[i]][2]*cr60+h12)+ - 1.0/(shape1[type[j]][0]*cr60+h12)+ - 1.0/(shape1[type[j]][1]*cr60+h12)+ - 1.0/(shape1[type[j]][2]*cr60+h12); - hsec = h12+b_alpha*sec; - dspr = 7.0/h12-1.0/hsec+stemp; - pbsr = b_alpha*sigma[type[i]][type[j]]/hsec; + stemp = 1.0 / (shape1[type[i]][0] * cr60 + h12) + 1.0 / (shape1[type[i]][1] * cr60 + h12) + + 1.0 / (shape1[type[i]][2] * cr60 + h12) + 1.0 / (shape1[type[j]][0] * cr60 + h12) + + 1.0 / (shape1[type[j]][1] * cr60 + h12) + 1.0 / (shape1[type[j]][2] * cr60 + h12); + hsec = h12 + b_alpha * sec; + dspr = 7.0 / h12 - 1.0 / hsec + stemp; + pbsr = b_alpha * sigma[type[i]][type[j]] / hsec; - for (int i=0; i<3; i++) { - u[0] = -rhat[i]*rhat[0]; - u[1] = -rhat[i]*rhat[1]; - u[2] = -rhat[i]*rhat[2]; - u[i] += 1.0; + for (int m = 0; m < 3; m++) { + u[0] = -rhat[m] * rhat[0]; + u[1] = -rhat[m] * rhat[1]; + u[2] = -rhat[m] * rhat[2]; + u[m] += 1.0; u[0] /= rnorm; u[1] /= rnorm; u[2] /= rnorm; - MathExtra::matvec(wi.A,u,u1); - MathExtra::matvec(wj.A,u,u2); - dsigma1=MathExtra::dot3(u1,vsigma1); - dsigma2=MathExtra::dot3(u2,vsigma2); - dH12[0][0] = dsigma1*gsigma1[0][0]+dsigma2*gsigma2[0][0]; - dH12[0][1] = dsigma1*gsigma1[0][1]+dsigma2*gsigma2[0][1]; - dH12[0][2] = dsigma1*gsigma1[0][2]+dsigma2*gsigma2[0][2]; - dH12[1][0] = dsigma1*gsigma1[1][0]+dsigma2*gsigma2[1][0]; - dH12[1][1] = dsigma1*gsigma1[1][1]+dsigma2*gsigma2[1][1]; - dH12[1][2] = dsigma1*gsigma1[1][2]+dsigma2*gsigma2[1][2]; - dH12[2][0] = dsigma1*gsigma1[2][0]+dsigma2*gsigma2[2][0]; - dH12[2][1] = dsigma1*gsigma1[2][1]+dsigma2*gsigma2[2][1]; - dH12[2][2] = dsigma1*gsigma1[2][2]+dsigma2*gsigma2[2][2]; - ddH = det_prime(H12,dH12); - deta = (dsigma1+dsigma2)*tsig1sig2; - deta -= ddH*tdH; - deta -= dsigma1*teta1+dsigma2*teta2; - dchi = MathExtra::dot3(u,fourw); - dh12 = rhat[i]+MathExtra::dot3(u,spr); - dUa = pbsu*(eta*dchi+deta*chi)-dh12*dspu; - dUr = pbsr*(eta*dchi+deta*chi)-dh12*dspr; - fforce[i]=dUr*Ur+dUa*Ua; + MathExtra::matvec(wi.A, u, u1); + MathExtra::matvec(wj.A, u, u2); + dsigma1 = MathExtra::dot3(u1, vsigma1); + dsigma2 = MathExtra::dot3(u2, vsigma2); + dH12[0][0] = dsigma1 * gsigma1[0][0] + dsigma2 * gsigma2[0][0]; + dH12[0][1] = dsigma1 * gsigma1[0][1] + dsigma2 * gsigma2[0][1]; + dH12[0][2] = dsigma1 * gsigma1[0][2] + dsigma2 * gsigma2[0][2]; + dH12[1][0] = dsigma1 * gsigma1[1][0] + dsigma2 * gsigma2[1][0]; + dH12[1][1] = dsigma1 * gsigma1[1][1] + dsigma2 * gsigma2[1][1]; + dH12[1][2] = dsigma1 * gsigma1[1][2] + dsigma2 * gsigma2[1][2]; + dH12[2][0] = dsigma1 * gsigma1[2][0] + dsigma2 * gsigma2[2][0]; + dH12[2][1] = dsigma1 * gsigma1[2][1] + dsigma2 * gsigma2[2][1]; + dH12[2][2] = dsigma1 * gsigma1[2][2] + dsigma2 * gsigma2[2][2]; + ddH = det_prime(H12, dH12); + deta = (dsigma1 + dsigma2) * tsig1sig2; + deta -= ddH * tdH; + deta -= dsigma1 * teta1 + dsigma2 * teta2; + dchi = MathExtra::dot3(u, fourw); + dh12 = rhat[m] + MathExtra::dot3(u, spr); + dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu; + dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr; + fforce[m] = dUr * Ur + dUa * Ua; } // torque on i - MathExtra::vecmat(fourw,wi.aTe,fwae); + MathExtra::vecmat(fourw, wi.aTe, fwae); - for (int i=0; i<3; i++) { - MathExtra::matvec(wi.lA[i],rhat,p); - dsigma1 = MathExtra::dot3(p,vsigma1); - dH12[0][0] = wi.lAsa[i][0][0]/sigma1+dsigma1*gsigma1[0][0]; - dH12[0][1] = wi.lAsa[i][0][1]/sigma1+dsigma1*gsigma1[0][1]; - dH12[0][2] = wi.lAsa[i][0][2]/sigma1+dsigma1*gsigma1[0][2]; - dH12[1][0] = wi.lAsa[i][1][0]/sigma1+dsigma1*gsigma1[1][0]; - dH12[1][1] = wi.lAsa[i][1][1]/sigma1+dsigma1*gsigma1[1][1]; - dH12[1][2] = wi.lAsa[i][1][2]/sigma1+dsigma1*gsigma1[1][2]; - dH12[2][0] = wi.lAsa[i][2][0]/sigma1+dsigma1*gsigma1[2][0]; - dH12[2][1] = wi.lAsa[i][2][1]/sigma1+dsigma1*gsigma1[2][1]; - dH12[2][2] = wi.lAsa[i][2][2]/sigma1+dsigma1*gsigma1[2][2]; - ddH = det_prime(H12,dH12); - deta = tsig1sig2*dsigma1-tdH*ddH; - deta -= teta1*dsigma1; + for (int i = 0; i < 3; i++) { + MathExtra::matvec(wi.lA[i], rhat, p); + dsigma1 = MathExtra::dot3(p, vsigma1); + dH12[0][0] = wi.lAsa[i][0][0] / sigma1 + dsigma1 * gsigma1[0][0]; + dH12[0][1] = wi.lAsa[i][0][1] / sigma1 + dsigma1 * gsigma1[0][1]; + dH12[0][2] = wi.lAsa[i][0][2] / sigma1 + dsigma1 * gsigma1[0][2]; + dH12[1][0] = wi.lAsa[i][1][0] / sigma1 + dsigma1 * gsigma1[1][0]; + dH12[1][1] = wi.lAsa[i][1][1] / sigma1 + dsigma1 * gsigma1[1][1]; + dH12[1][2] = wi.lAsa[i][1][2] / sigma1 + dsigma1 * gsigma1[1][2]; + dH12[2][0] = wi.lAsa[i][2][0] / sigma1 + dsigma1 * gsigma1[2][0]; + dH12[2][1] = wi.lAsa[i][2][1] / sigma1 + dsigma1 * gsigma1[2][1]; + dH12[2][2] = wi.lAsa[i][2][2] / sigma1 + dsigma1 * gsigma1[2][2]; + ddH = det_prime(H12, dH12); + deta = tsig1sig2 * dsigma1 - tdH * ddH; + deta -= teta1 * dsigma1; double tempv[3]; - MathExtra::matvec(wi.lA[i],w,tempv); - dchi = -MathExtra::dot3(fwae,tempv); - MathExtra::matvec(wi.lAtwo[i],spr,tempv); - dh12 = -MathExtra::dot3(s,tempv); + MathExtra::matvec(wi.lA[i], w, tempv); + dchi = -MathExtra::dot3(fwae, tempv); + MathExtra::matvec(wi.lAtwo[i], spr, tempv); + dh12 = -MathExtra::dot3(s, tempv); - dUa = pbsu*(eta*dchi + deta*chi)-dh12*dspu; - dUr = pbsr*(eta*dchi + deta*chi)-dh12*dspr; - ttor[i] = -(dUa*Ua+dUr*Ur); + dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu; + dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr; + ttor[i] = -(dUa * Ua + dUr * Ur); } // torque on j - if (!(force->newton_pair || j < atom->nlocal)) - return Ua+Ur; + if (!(force->newton_pair || j < atom->nlocal)) return Ua + Ur; - MathExtra::vecmat(fourw,wj.aTe,fwae); + MathExtra::vecmat(fourw, wj.aTe, fwae); - for (int i=0; i<3; i++) { - MathExtra::matvec(wj.lA[i],rhat,p); - dsigma2 = MathExtra::dot3(p,vsigma2); - dH12[0][0] = wj.lAsa[i][0][0]/sigma2+dsigma2*gsigma2[0][0]; - dH12[0][1] = wj.lAsa[i][0][1]/sigma2+dsigma2*gsigma2[0][1]; - dH12[0][2] = wj.lAsa[i][0][2]/sigma2+dsigma2*gsigma2[0][2]; - dH12[1][0] = wj.lAsa[i][1][0]/sigma2+dsigma2*gsigma2[1][0]; - dH12[1][1] = wj.lAsa[i][1][1]/sigma2+dsigma2*gsigma2[1][1]; - dH12[1][2] = wj.lAsa[i][1][2]/sigma2+dsigma2*gsigma2[1][2]; - dH12[2][0] = wj.lAsa[i][2][0]/sigma2+dsigma2*gsigma2[2][0]; - dH12[2][1] = wj.lAsa[i][2][1]/sigma2+dsigma2*gsigma2[2][1]; - dH12[2][2] = wj.lAsa[i][2][2]/sigma2+dsigma2*gsigma2[2][2]; - ddH = det_prime(H12,dH12); - deta = tsig1sig2*dsigma2-tdH*ddH; - deta -= teta2*dsigma2; + for (int i = 0; i < 3; i++) { + MathExtra::matvec(wj.lA[i], rhat, p); + dsigma2 = MathExtra::dot3(p, vsigma2); + dH12[0][0] = wj.lAsa[i][0][0] / sigma2 + dsigma2 * gsigma2[0][0]; + dH12[0][1] = wj.lAsa[i][0][1] / sigma2 + dsigma2 * gsigma2[0][1]; + dH12[0][2] = wj.lAsa[i][0][2] / sigma2 + dsigma2 * gsigma2[0][2]; + dH12[1][0] = wj.lAsa[i][1][0] / sigma2 + dsigma2 * gsigma2[1][0]; + dH12[1][1] = wj.lAsa[i][1][1] / sigma2 + dsigma2 * gsigma2[1][1]; + dH12[1][2] = wj.lAsa[i][1][2] / sigma2 + dsigma2 * gsigma2[1][2]; + dH12[2][0] = wj.lAsa[i][2][0] / sigma2 + dsigma2 * gsigma2[2][0]; + dH12[2][1] = wj.lAsa[i][2][1] / sigma2 + dsigma2 * gsigma2[2][1]; + dH12[2][2] = wj.lAsa[i][2][2] / sigma2 + dsigma2 * gsigma2[2][2]; + ddH = det_prime(H12, dH12); + deta = tsig1sig2 * dsigma2 - tdH * ddH; + deta -= teta2 * dsigma2; double tempv[3]; - MathExtra::matvec(wj.lA[i],w,tempv); - dchi = -MathExtra::dot3(fwae,tempv); - MathExtra::matvec(wj.lAtwo[i],spr,tempv); - dh12 = -MathExtra::dot3(s,tempv); + MathExtra::matvec(wj.lA[i], w, tempv); + dchi = -MathExtra::dot3(fwae, tempv); + MathExtra::matvec(wj.lAtwo[i], spr, tempv); + dh12 = -MathExtra::dot3(s, tempv); - dUa = pbsu*(eta*dchi + deta*chi)-dh12*dspu; - dUr = pbsr*(eta*dchi + deta*chi)-dh12*dspr; - rtor[i] = -(dUa*Ua+dUr*Ur); + dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu; + dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr; + rtor[i] = -(dUa * Ua + dUr * Ur); } - return Ua+Ur; + return Ua + Ur; } /* ---------------------------------------------------------------------- Compute the energy, force, torque for a pair (INTEGRATED-LJ) ------------------------------------------------------------------------- */ -double PairRESquared::resquared_lj(const int i, const int j, - const RE2Vars &wi, const double *r, - const double rsq, double *fforce, - double *ttor, bool calc_torque) +double PairRESquared::resquared_lj(const int i, const int j, const RE2Vars &wi, const double *r, + const double rsq, double *fforce, double *ttor, bool calc_torque) { int *type = atom->type; // pair computations for energy, force, torque - double rnorm; // L2 norm of r - double rhat[3]; // r/rnorm - double s[3]; // inv(gamma1)*rhat - double sigma12; // 1/sqrt(0.5*s'*rhat) - double w[3]; // inv(A1'*E1*A1+I)*rhat - double h12; // rnorm-sigma12; - double chi; // 2*rhat'*w - double sigh; // sigma/h12 - double tprod; // chi*sigh - double Ua,Ur; // attractive/repulsive parts of potential + double rnorm; // L2 norm of r + double rhat[3]; // r/rnorm + double s[3]; // inv(gamma1)*rhat + double sigma12; // 1/sqrt(0.5*s'*rhat) + double w[3]; // inv(A1'*E1*A1+I)*rhat + double h12; // rnorm-sigma12; + double chi; // 2*rhat'*w + double sigh; // sigma/h12 + double tprod; // chi*sigh + double Ua, Ur; // attractive/repulsive parts of potential // pair computations for force, torque - double sec; // sigma*chi - double sigma12p3; // sigma12^3 - double fourw[3]; // 4*w; - double spr[3]; // 0.5*sigma12^3*s - double hsec; // h12+[3,b_alpha]*sec - double dspu; // 1/h12 - 1/hsec + temp - double pbsu; // 3*sigma/hsec - double dspr; // 7/h12-1/hsec+temp - double pbsr; // b_alpha*sigma/hsec; - double u[3]; // (-rhat(i)*rhat+eye(:,i))/rnorm - double dchi,dh12; // derivatives of chi,h12 - double dUr,dUa; // derivatives of Ua,Ur - double h12p3; // h12^3 + double sec; // sigma*chi + double sigma12p3; // sigma12^3 + double fourw[3]; // 4*w; + double spr[3]; // 0.5*sigma12^3*s + double hsec; // h12+[3,b_alpha]*sec + double dspu; // 1/h12 - 1/hsec + temp + double pbsu; // 3*sigma/hsec + double dspr; // 7/h12-1/hsec+temp + double pbsr; // b_alpha*sigma/hsec; + double u[3]; // (-rhat(i)*rhat+eye(:,i))/rnorm + double dchi, dh12; // derivatives of chi,h12 + double dUr, dUa; // derivatives of Ua,Ur + double h12p3; // h12^3 // pair computations for torque - double fwae[3]; // -fourw'*aTe - double p[3]; // lA*rhat + double fwae[3]; // -fourw'*aTe + double p[3]; // lA*rhat // distance of closest approach correction - double aTs[3][3]; // A1'*S1^2 - double gamma[3][3]; // A1'*S1^2*A - double lAtwo[3][3][3]; // A1'*S1^2*wi.lA + double aTs[3][3]; // A1'*S1^2 + double gamma[3][3]; // A1'*S1^2*A + double lAtwo[3][3][3]; // A1'*S1^2*wi.lA double scorrect[3]; - double half_sigma=sigma[type[i]][type[j]] / 2.0; - scorrect[0] = shape1[type[i]][0]+half_sigma; - scorrect[1] = shape1[type[i]][1]+half_sigma; - scorrect[2] = shape1[type[i]][2]+half_sigma; + double half_sigma = sigma[type[i]][type[j]] / 2.0; + scorrect[0] = shape1[type[i]][0] + half_sigma; + scorrect[1] = shape1[type[i]][1] + half_sigma; + scorrect[2] = shape1[type[i]][2] + half_sigma; scorrect[0] = scorrect[0] * scorrect[0] / 2.0; scorrect[1] = scorrect[1] * scorrect[1] / 2.0; scorrect[2] = scorrect[2] * scorrect[2] / 2.0; - MathExtra::transpose_diag3(wi.A,scorrect,aTs); - MathExtra::times3(aTs,wi.A,gamma); - for (int ii=0; ii<3; ii++) - MathExtra::times3(aTs,wi.lA[ii],lAtwo[ii]); + MathExtra::transpose_diag3(wi.A, scorrect, aTs); + MathExtra::times3(aTs, wi.A, gamma); + for (int ii = 0; ii < 3; ii++) MathExtra::times3(aTs, wi.lA[ii], lAtwo[ii]); - rnorm=sqrt(rsq); - rhat[0] = r[0]/rnorm; - rhat[1] = r[1]/rnorm; - rhat[2] = r[2]/rnorm; + rnorm = sqrt(rsq); + rhat[0] = r[0] / rnorm; + rhat[1] = r[1] / rnorm; + rhat[2] = r[2] / rnorm; // energy - int ierror = MathExtra::mldivide3(gamma,rhat,s); - if (ierror) error->all(FLERR,"Bad matrix inversion in mldivide3"); + int ierror = MathExtra::mldivide3(gamma, rhat, s); + if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3"); - sigma12 = 1.0/sqrt(0.5*MathExtra::dot3(s,rhat)); + sigma12 = 1.0 / sqrt(0.5 * MathExtra::dot3(s, rhat)); double temp[3][3]; - MathExtra::times3(wi.aTe,wi.A,temp); + MathExtra::times3(wi.aTe, wi.A, temp); temp[0][0] += 1.0; temp[1][1] += 1.0; temp[2][2] += 1.0; - ierror = MathExtra::mldivide3(temp,rhat,w); - if (ierror) error->all(FLERR,"Bad matrix inversion in mldivide3"); + ierror = MathExtra::mldivide3(temp, rhat, w); + if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3"); - h12 = rnorm-sigma12; - chi = 2.0*MathExtra::dot3(rhat,w); - sigh = sigma[type[i]][type[j]]/h12; - tprod = chi*sigh; + h12 = rnorm - sigma12; + chi = 2.0 * MathExtra::dot3(rhat, w); + sigh = sigma[type[i]][type[j]] / h12; + tprod = chi * sigh; - h12p3 = pow(h12,3.0); - double sigmap3 = pow(sigma[type[i]][type[j]],3.0); - double stemp = h12/2.0; - Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* - (shape1[type[i]][2]+stemp)*h12p3/8.0; - Ua = (1.0+3.0*tprod)*lshape[type[i]]/Ua; - Ua = epsilon[type[i]][type[j]]*Ua*sigmap3*solv_f_a; + h12p3 = pow(h12, 3.0); + double sigmap3 = pow(sigma[type[i]][type[j]], 3.0); + double stemp = h12 / 2.0; + Ua = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) * + h12p3 / 8.0; + Ua = (1.0 + 3.0 * tprod) * lshape[type[i]] / Ua; + Ua = epsilon[type[i]][type[j]] * Ua * sigmap3 * solv_f_a; - stemp = h12/cr60; - Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)* - (shape1[type[i]][2]+stemp)*h12p3/60.0; - Ur = (1.0+b_alpha*tprod)*lshape[type[i]]/Ur; - Ur = epsilon[type[i]][type[j]]*Ur*sigmap3*pow(sigh,6.0)*solv_f_r; + stemp = h12 / cr60; + Ur = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) * + h12p3 / 60.0; + Ur = (1.0 + b_alpha * tprod) * lshape[type[i]] / Ur; + Ur = epsilon[type[i]][type[j]] * Ur * sigmap3 * pow(sigh, 6.0) * solv_f_r; // force - sec = sigma[type[i]][type[j]]*chi; - sigma12p3 = pow(sigma12,3.0); - fourw[0] = 4.0*w[0]; - fourw[1] = 4.0*w[1]; - fourw[2] = 4.0*w[2]; - spr[0] = 0.5*sigma12p3*s[0]; - spr[1] = 0.5*sigma12p3*s[1]; - spr[2] = 0.5*sigma12p3*s[2]; + sec = sigma[type[i]][type[j]] * chi; + sigma12p3 = pow(sigma12, 3.0); + fourw[0] = 4.0 * w[0]; + fourw[1] = 4.0 * w[1]; + fourw[2] = 4.0 * w[2]; + spr[0] = 0.5 * sigma12p3 * s[0]; + spr[1] = 0.5 * sigma12p3 * s[1]; + spr[2] = 0.5 * sigma12p3 * s[2]; - stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+ - 1.0/(shape1[type[i]][1]*2.0+h12)+ - 1.0/(shape1[type[i]][2]*2.0+h12)+ - 3.0/h12; - hsec = h12+3.0*sec; - dspu = 1.0/h12-1.0/hsec+stemp; - pbsu = 3.0*sigma[type[i]][type[j]]/hsec; + stemp = 1.0 / (shape1[type[i]][0] * 2.0 + h12) + 1.0 / (shape1[type[i]][1] * 2.0 + h12) + + 1.0 / (shape1[type[i]][2] * 2.0 + h12) + 3.0 / h12; + hsec = h12 + 3.0 * sec; + dspu = 1.0 / h12 - 1.0 / hsec + stemp; + pbsu = 3.0 * sigma[type[i]][type[j]] / hsec; - stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+ - 1.0/(shape1[type[i]][1]*cr60+h12)+ - 1.0/(shape1[type[i]][2]*cr60+h12)+ - 3.0/h12; - hsec = h12+b_alpha*sec; - dspr = 7.0/h12-1.0/hsec+stemp; - pbsr = b_alpha*sigma[type[i]][type[j]]/hsec; + stemp = 1.0 / (shape1[type[i]][0] * cr60 + h12) + 1.0 / (shape1[type[i]][1] * cr60 + h12) + + 1.0 / (shape1[type[i]][2] * cr60 + h12) + 3.0 / h12; + hsec = h12 + b_alpha * sec; + dspr = 7.0 / h12 - 1.0 / hsec + stemp; + pbsr = b_alpha * sigma[type[i]][type[j]] / hsec; - for (int i=0; i<3; i++) { - u[0] = -rhat[i]*rhat[0]; - u[1] = -rhat[i]*rhat[1]; - u[2] = -rhat[i]*rhat[2]; - u[i] += 1.0; + for (int m = 0; m < 3; m++) { + u[0] = -rhat[m] * rhat[0]; + u[1] = -rhat[m] * rhat[1]; + u[2] = -rhat[m] * rhat[2]; + u[m] += 1.0; u[0] /= rnorm; u[1] /= rnorm; u[2] /= rnorm; - dchi = MathExtra::dot3(u,fourw); - dh12 = rhat[i]+MathExtra::dot3(u,spr); - dUa = pbsu*dchi-dh12*dspu; - dUr = pbsr*dchi-dh12*dspr; - fforce[i]=dUr*Ur+dUa*Ua; + dchi = MathExtra::dot3(u, fourw); + dh12 = rhat[m] + MathExtra::dot3(u, spr); + dUa = pbsu * dchi - dh12 * dspu; + dUr = pbsr * dchi - dh12 * dspr; + fforce[m] = dUr * Ur + dUa * Ua; } // torque on i if (calc_torque) { - MathExtra::vecmat(fourw,wi.aTe,fwae); + MathExtra::vecmat(fourw, wi.aTe, fwae); - for (int i=0; i<3; i++) { - MathExtra::matvec(wi.lA[i],rhat,p); + for (int m = 0; m < 3; m++) { + MathExtra::matvec(wi.lA[m], rhat, p); double tempv[3]; - MathExtra::matvec(wi.lA[i],w,tempv); - dchi = -MathExtra::dot3(fwae,tempv); - MathExtra::matvec(lAtwo[i],spr,tempv); - dh12 = -MathExtra::dot3(s,tempv); + MathExtra::matvec(wi.lA[m], w, tempv); + dchi = -MathExtra::dot3(fwae, tempv); + MathExtra::matvec(lAtwo[m], spr, tempv); + dh12 = -MathExtra::dot3(s, tempv); - dUa = pbsu*dchi-dh12*dspu; - dUr = pbsr*dchi-dh12*dspr; - ttor[i] = -(dUa*Ua+dUr*Ur); + dUa = pbsu * dchi - dh12 * dspu; + dUr = pbsr * dchi - dh12 * dspr; + ttor[m] = -(dUa * Ua + dUr * Ur); } } - return Ua+Ur; + return Ua + Ur; } diff --git a/src/BODY/pair_body_rounded_polygon.cpp b/src/BODY/pair_body_rounded_polygon.cpp index 7dc7956586..23a4da7f07 100644 --- a/src/BODY/pair_body_rounded_polygon.cpp +++ b/src/BODY/pair_body_rounded_polygon.cpp @@ -21,21 +21,21 @@ #include "pair_body_rounded_polygon.h" -#include -#include -#include "math_extra.h" #include "atom.h" #include "atom_vec_body.h" #include "body_rounded_polygon.h" #include "comm.h" -#include "force.h" -#include "fix.h" -#include "modify.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" #include "error.h" +#include "fix.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include using namespace LAMMPS_NS; diff --git a/src/BODY/pair_body_rounded_polyhedron.cpp b/src/BODY/pair_body_rounded_polyhedron.cpp index 65179ecb80..6deb06aafb 100644 --- a/src/BODY/pair_body_rounded_polyhedron.cpp +++ b/src/BODY/pair_body_rounded_polyhedron.cpp @@ -23,22 +23,22 @@ #include "pair_body_rounded_polyhedron.h" -#include -#include #include "atom.h" #include "atom_vec_body.h" #include "body_rounded_polyhedron.h" #include "comm.h" -#include "force.h" -#include "fix.h" -#include "modify.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" #include "error.h" -#include "math_extra.h" +#include "fix.h" +#include "force.h" #include "math_const.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -1212,10 +1212,7 @@ int PairBodyRoundedPolyhedron::interaction_edge_to_edge(int ibody, contact_list[num_contacts].unique = 1; num_contacts++; } - } else { - } - return interact; } diff --git a/src/DPD-REACT/pair_exp6_rx.cpp b/src/DPD-REACT/pair_exp6_rx.cpp index 985fd551fd..1d242de90b 100644 --- a/src/DPD-REACT/pair_exp6_rx.cpp +++ b/src/DPD-REACT/pair_exp6_rx.cpp @@ -617,8 +617,7 @@ void PairExp6rx::coeff(int narg, char **arg) if (strcmp(site1,"1fluid") == 0) isite1 = oneFluidApproxParameter; - else - { + else { int isp; for (isp = 0; isp < nspecies; isp++) if (strcmp(site1, &atom->dvname[isp][0]) == 0) break; @@ -631,8 +630,7 @@ void PairExp6rx::coeff(int narg, char **arg) if (strcmp(site2,"1fluid") == 0) isite2 = oneFluidApproxParameter; - else - { + else { int isp; for (isp = 0; isp < nspecies; isp++) if (strcmp(site2, &atom->dvname[isp][0]) == 0) break; @@ -644,8 +642,7 @@ void PairExp6rx::coeff(int narg, char **arg) } // Set the interaction potential type to the enumerated type. - for (int iparam = 0; iparam < nparams; ++iparam) - { + for (int iparam = 0; iparam < nparams; ++iparam) { if (strcmp( params[iparam].potential, "exp6") == 0) params[iparam].potentialType = exp6PotentialType; else @@ -663,7 +660,6 @@ void PairExp6rx::coeff(int narg, char **arg) scalingFlag = EXPONENT; exponentR = utils::numeric(FLERR,arg[6],false,lmp); exponentEpsilon = utils::numeric(FLERR,arg[7],false,lmp); - if (narg > 9) error->all(FLERR,"Incorrect args for pair coefficients"); if (narg == 9) cut_one = utils::numeric(FLERR,arg[8],false,lmp); } else if (strcmp(arg[5],"polynomial") == 0) { scalingFlag = POLYNOMIAL; diff --git a/src/EXTRA-FIX/fix_gld.cpp b/src/EXTRA-FIX/fix_gld.cpp index 34ed43301e..2e30882206 100644 --- a/src/EXTRA-FIX/fix_gld.cpp +++ b/src/EXTRA-FIX/fix_gld.cpp @@ -132,6 +132,7 @@ FixGLD::FixGLD(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Illegal fix gld command"); } if (strcmp(arg[iarg+1],"no") == 0) { + zeroflag = 0; } else if (strcmp(arg[iarg+1],"yes") == 0) { zeroflag = 1; } else { @@ -144,6 +145,7 @@ FixGLD::FixGLD(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Illegal fix gld command"); } if (strcmp(arg[iarg+1],"no") == 0) { + freezeflag = 0; } else if (strcmp(arg[iarg+1],"yes") == 0) { freezeflag = 1; for (int i = 0; i < atom->nlocal; i++) { diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp index b98b105045..8e3ec2ace8 100644 --- a/src/GPU/pppm_gpu.cpp +++ b/src/GPU/pppm_gpu.cpp @@ -106,6 +106,8 @@ PPPMGPU::PPPMGPU(LAMMPS *lmp) : PPPM(lmp) PPPMGPU::~PPPMGPU() { PPPM_GPU_API(clear)(poisson_time); + destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out); + destroy_3d_offset(vd_brick,nzlo_out,nylo_out); } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index a604782993..5f16176e49 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -437,7 +437,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // register with Atom class history_one = nullptr; - grow_arrays(atom->nmax); + FixWallGran::grow_arrays(atom->nmax); atom->add_callback(Atom::GROW); atom->add_callback(Atom::RESTART); @@ -1555,8 +1555,7 @@ double FixWallGran::memory_usage() void FixWallGran::grow_arrays(int nmax) { - if (use_history) memory->grow(history_one,nmax,size_history, - "fix_wall_gran:history_one"); + if (use_history) memory->grow(history_one,nmax,size_history,"fix_wall_gran:history_one"); if (peratom_flag) { memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); } diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 0a92bac543..37afe9632c 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -62,7 +62,7 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : ncontact = nullptr; walls = nullptr; history_many = nullptr; - grow_arrays(atom->nmax); + FixWallGranRegion::grow_arrays(atom->nmax); // initialize shear history as if particle is not touching region @@ -355,8 +355,7 @@ void FixWallGranRegion::grow_arrays(int nmax) if (use_history) { memory->grow(ncontact,nmax,"fix_wall_gran:ncontact"); memory->grow(walls,nmax,tmax,"fix_wall_gran:walls"); - memory->grow(history_many,nmax,tmax,size_history, - "fix_wall_gran:history_many"); + memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many"); } if (peratom_flag) memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index d0e75d4812..2846403e4c 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -109,7 +109,8 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) PairGranular::~PairGranular() { - delete [] svector; + delete[] svector; + delete[] history_transfer_factors; if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY"); else modify->delete_fix("NEIGH_HISTORY_GRANULAR"); diff --git a/src/MANYBODY/pair_tersoff_table.cpp b/src/MANYBODY/pair_tersoff_table.cpp index e64e5bce02..94c4d88b09 100644 --- a/src/MANYBODY/pair_tersoff_table.cpp +++ b/src/MANYBODY/pair_tersoff_table.cpp @@ -91,7 +91,7 @@ PairTersoffTable::~PairTersoffTable() memory->destroy(cutsq); } deallocateGrids(); - deallocatePreLoops(); + PairTersoffTable::deallocatePreLoops(); } /* ---------------------------------------------------------------------- */ diff --git a/src/OPENMP/pair_tersoff_table_omp.cpp b/src/OPENMP/pair_tersoff_table_omp.cpp index 47ee2feca8..5367a3d043 100644 --- a/src/OPENMP/pair_tersoff_table_omp.cpp +++ b/src/OPENMP/pair_tersoff_table_omp.cpp @@ -53,7 +53,7 @@ PairTersoffTableOMP::PairTersoffTableOMP(LAMMPS *lmp) : PairTersoffTableOMP::~PairTersoffTableOMP() { if (allocated) { - deallocatePreLoops(); + PairTersoffTableOMP::deallocatePreLoops(); } } diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index fd197a156a..8c3b7df190 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -165,7 +165,7 @@ FixQEqReaxFF::~FixQEqReaxFF() memory->destroy(s_hist); memory->destroy(t_hist); - deallocate_storage(); + FixQEqReaxFF::deallocate_storage(); deallocate_matrix(); memory->destroy(shld); diff --git a/src/REPLICA/fix_pimd.cpp b/src/REPLICA/fix_pimd.cpp index ce1d41302f..bb4af16f1b 100644 --- a/src/REPLICA/fix_pimd.cpp +++ b/src/REPLICA/fix_pimd.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -24,69 +23,29 @@ #include "fix_pimd.h" -#include -#include - -#include "universe.h" -#include "comm.h" -#include "force.h" #include "atom.h" +#include "comm.h" #include "domain.h" -#include "update.h" +#include "error.h" +#include "force.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "universe.h" +#include "update.h" + +#include +#include using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{PIMD,NMPIMD,CMD}; +enum { PIMD, NMPIMD, CMD }; /* ---------------------------------------------------------------------- */ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - method = PIMD; - fmass = 1.0; - nhc_temp = 298.15; - nhc_nchain = 2; - sp = 1.0; - - for (int i=3; iuniverse_all(FLERR,"Unknown method parameter for fix pimd"); - } - else if (strcmp(arg[i],"fmass")==0) - { - fmass = atof(arg[i+1]); - if (fmass<0.0 || fmass>1.0) error->universe_all(FLERR,"Invalid fmass value for fix pimd"); - } - else if (strcmp(arg[i],"sp")==0) - { - sp = atof(arg[i+1]); - if (fmass<0.0) error->universe_all(FLERR,"Invalid sp value for fix pimd"); - } - else if (strcmp(arg[i],"temp")==0) - { - nhc_temp = atof(arg[i+1]); - if (nhc_temp<0.0) error->universe_all(FLERR,"Invalid temp value for fix pimd"); - } - else if (strcmp(arg[i],"nhc")==0) - { - nhc_nchain = atoi(arg[i+1]); - if (nhc_nchain<2) error->universe_all(FLERR,"Invalid nhc value for fix pimd"); - } - else error->universe_all(arg[i],i+1,"Unknown keyword for fix pimd"); - } - - /* Initiation */ - max_nsend = 0; tag_send = nullptr; buf_send = nullptr; @@ -110,25 +69,60 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) nhc_eta_dotdot = nullptr; nhc_eta_mass = nullptr; + method = PIMD; + fmass = 1.0; + nhc_temp = 298.15; + nhc_nchain = 2; + sp = 1.0; + + for (int i = 3; i < narg - 1; i += 2) { + if (strcmp(arg[i], "method") == 0) { + if (strcmp(arg[i + 1], "pimd") == 0) + method = PIMD; + else if (strcmp(arg[i + 1], "nmpimd") == 0) + method = NMPIMD; + else if (strcmp(arg[i + 1], "cmd") == 0) + method = CMD; + else + error->universe_all(FLERR, "Unknown method parameter for fix pimd"); + } else if (strcmp(arg[i], "fmass") == 0) { + fmass = utils::numeric(FLERR, arg[i + 1], false, lmp); + if (fmass < 0.0 || fmass > 1.0) + error->universe_all(FLERR, "Invalid fmass value for fix pimd"); + } else if (strcmp(arg[i], "sp") == 0) { + sp = utils::numeric(FLERR, arg[i + 1], false, lmp); + if (fmass < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd"); + } else if (strcmp(arg[i], "temp") == 0) { + nhc_temp = utils::numeric(FLERR, arg[i + 1], false, lmp); + if (nhc_temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd"); + } else if (strcmp(arg[i], "nhc") == 0) { + nhc_nchain = utils::inumeric(FLERR, arg[i + 1], false, lmp); + if (nhc_nchain < 2) error->universe_all(FLERR, "Invalid nhc value for fix pimd"); + } else + error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd", arg[i])); + } + + /* Initiation */ + size_peratom_cols = 12 * nhc_nchain + 3; nhc_offset_one_1 = 3 * nhc_nchain; - nhc_offset_one_2 = 3 * nhc_nchain +3; + nhc_offset_one_2 = 3 * nhc_nchain + 3; nhc_size_one_1 = sizeof(double) * nhc_offset_one_1; nhc_size_one_2 = sizeof(double) * nhc_offset_one_2; restart_peratom = 1; - peratom_flag = 1; - peratom_freq = 1; + peratom_flag = 1; + peratom_freq = 1; global_freq = 1; vector_flag = 1; size_vector = 2; - extvector = 1; + extvector = 1; comm_forward = 3; - atom->add_callback(Atom::GROW); // Call LAMMPS to allocate memory for per-atom array - atom->add_callback(Atom::RESTART); // Call LAMMPS to re-assign restart-data for per-atom array + atom->add_callback(Atom::GROW); // Call LAMMPS to allocate memory for per-atom array + atom->add_callback(Atom::RESTART); // Call LAMMPS to re-assign restart-data for per-atom array grow_arrays(atom->nmax); @@ -138,7 +132,37 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } /* ---------------------------------------------------------------------- */ +FixPIMD::~FixPIMD() +{ + delete[] mass; + atom->delete_callback(id, Atom::GROW); + atom->delete_callback(id, Atom::RESTART); + memory->destroy(M_x2xp); + memory->destroy(M_xp2x); + memory->destroy(M_f2fp); + memory->destroy(M_fp2f); + memory->sfree(lam); + + if (buf_beads) + for (int i = 0; i < np; i++) memory->sfree(buf_beads[i]); + delete[] buf_beads; + delete[] plan_send; + delete[] plan_recv; + delete[] mode_index; + + memory->sfree(tag_send); + memory->sfree(buf_send); + memory->sfree(buf_recv); + + memory->destroy(array_atom); + memory->destroy(nhc_eta); + memory->destroy(nhc_eta_dot); + memory->destroy(nhc_eta_dotdot); + memory->destroy(nhc_eta_mass); +} + +/* ---------------------------------------------------------------------- */ int FixPIMD::setmask() { int mask = 0; @@ -153,9 +177,10 @@ int FixPIMD::setmask() void FixPIMD::init() { if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Fix pimd requires an atom map, see atom_modify"); + error->all(FLERR, "Fix pimd requires an atom map, see atom_modify"); - if (universe->me==0 && screen) fprintf(screen,"Fix pimd initializing Path-Integral ...\n"); + if (universe->me == 0 && universe->uscreen) + fprintf(universe->uscreen, "Fix pimd initializing Path-Integral ...\n"); // prepare the constants @@ -181,16 +206,16 @@ void FixPIMD::init() /* The current solution, using LAMMPS internal real units */ const double Boltzmann = force->boltz; - const double Plank = force->hplanck; + const double Plank = force->hplanck; - double hbar = Plank / ( 2.0 * MY_PI ); - double beta = 1.0 / (Boltzmann * nhc_temp); - double _fbond = 1.0 * np / (beta*beta*hbar*hbar) ; + double hbar = Plank / (2.0 * MY_PI); + double beta = 1.0 / (Boltzmann * nhc_temp); + double _fbond = 1.0 * np / (beta * beta * hbar * hbar); omega_np = sqrt(np) / (hbar * beta) * sqrt(force->mvv2e); - fbond = - _fbond * force->mvv2e; + fbond = -_fbond * force->mvv2e; - if (universe->me==0) + if (universe->me == 0) printf("Fix pimd -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)\n\n", fbond); dtv = update->dt; @@ -198,10 +223,12 @@ void FixPIMD::init() comm_init(); - mass = new double [atom->ntypes+1]; + mass = new double[atom->ntypes + 1]; - if (method==CMD || method==NMPIMD) nmpimd_init(); - else for (int i=1; i<=atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass; + if (method == CMD || method == NMPIMD) + nmpimd_init(); + else + for (int i = 1; i <= atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass; if (!nhc_ready) nhc_init(); } @@ -210,7 +237,8 @@ void FixPIMD::init() void FixPIMD::setup(int vflag) { - if (universe->me==0 && screen) fprintf(screen,"Setting up Path-Integral ...\n"); + if (universe->me == 0 && universe->uscreen) + fprintf(universe->uscreen, "Setting up Path-Integral ...\n"); post_force(vflag); } @@ -234,13 +262,13 @@ void FixPIMD::final_integrate() void FixPIMD::post_force(int /*flag*/) { - for (int i=0; inlocal; i++) for(int j=0; j<3; j++) atom->f[i][j] /= np; + for (int i = 0; i < atom->nlocal; i++) + for (int j = 0; j < 3; j++) atom->f[i][j] /= np; comm_exec(atom->x); spring_force(); - if (method==CMD || method==NMPIMD) - { + if (method == CMD || method == NMPIMD) { /* forward comm for the force on ghost atoms */ nmpimd_fill(atom->f); @@ -262,34 +290,38 @@ void FixPIMD::post_force(int /*flag*/) void FixPIMD::nhc_init() { double tau = 1.0 / omega_np; - double KT = force->boltz * nhc_temp; + double KT = force->boltz * nhc_temp; double mass0 = KT * tau * tau; int max = 3 * atom->nlocal; - for (int i=0; iiworld==0) ; else nhc_eta_mass[i][ichain] *= fmass; + nhc_eta_mass[i][ichain] = mass0; + if ((method == CMD || method == NMPIMD) && universe->iworld == 0) + ; + else + nhc_eta_mass[i][ichain] *= fmass; } - nhc_eta_dot[i][nhc_nchain] = 0.0; + nhc_eta_dot[i][nhc_nchain] = 0.0; - for (int ichain=1; ichainmvv2e - KT) / nhc_eta_mass[i][ichain]; + for (int ichain = 1; ichain < nhc_nchain; ichain++) + nhc_eta_dotdot[i][ichain] = (nhc_eta_mass[i][ichain - 1] * nhc_eta_dot[i][ichain - 1] * + nhc_eta_dot[i][ichain - 1] * force->mvv2e - + KT) / + nhc_eta_mass[i][ichain]; } // Zero NH acceleration for CMD - if (method==CMD && universe->iworld==0) for (int i=0; iiworld == 0) + for (int i = 0; i < max; i++) + for (int ichain = 0; ichain < nhc_nchain; ichain++) nhc_eta_dotdot[i][ichain] = 0.0; nhc_ready = true; } @@ -302,8 +334,7 @@ void FixPIMD::nhc_update_x() double **x = atom->x; double **v = atom->v; - if (method==CMD || method==NMPIMD) - { + if (method == CMD || method == NMPIMD) { nmpimd_fill(atom->v); comm_exec(atom->v); @@ -313,8 +344,7 @@ void FixPIMD::nhc_update_x() nmpimd_transform(buf_beads, v, M_xp2x[universe->iworld]); } - for (int i=0; iv; double **f = atom->f; - for (int i=0; iiworld==0) return; + if (method == CMD && universe->iworld == 0) return; double expfac; int nmax = 3 * atom->nlocal; double KT = force->boltz * nhc_temp; double kecurrent, t_current; - double dthalf = 0.5 * update->dt; - double dt4 = 0.25 * update->dt; - double dt8 = 0.125 * update->dt; + double dthalf = 0.5 * update->dt; + double dt4 = 0.25 * update->dt; + double dt8 = 0.125 * update->dt; - for (int i=0; imvv2e; + kecurrent = mass[type[iatm]] * vv[idim] * vv[idim] * force->mvv2e; t_current = kecurrent / force->boltz; double *eta = nhc_eta[i]; @@ -366,9 +394,8 @@ void FixPIMD::nhc_update_v() eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0]; - for (int ichain=nhc_nchain-1; ichain>0; ichain--) - { - expfac = exp(-dt8 * eta_dot[ichain+1]); + for (int ichain = nhc_nchain - 1; ichain > 0; ichain--) { + expfac = exp(-dt8 * eta_dot[ichain + 1]); eta_dot[ichain] *= expfac; eta_dot[ichain] += eta_dotdot[ichain] * dt4; eta_dot[ichain] *= expfac; @@ -388,19 +415,18 @@ void FixPIMD::nhc_update_v() kecurrent = force->boltz * t_current; eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0]; - for (int ichain=0; ichaincreate(M_f2fp, np, np, "fix_feynman:M_f2fp"); memory->create(M_fp2f, np, np, "fix_feynman:M_fp2f"); - lam = (double*) memory->smalloc(sizeof(double)*np, "FixPIMD::lam"); + lam = (double *) memory->smalloc(sizeof(double) * np, "FixPIMD::lam"); // Set up eigenvalues lam[0] = 0.0; - if (np%2==0) lam[np-1] = 4.0 * np; + if (np % 2 == 0) lam[np - 1] = 4.0 * np; - for (int i=2; i<=np/2; i++) - { - lam[2*i-3] = lam[2*i-2] = 2.0 * np * (1.0 - 1.0 *cos(2.0*MY_PI*(i-1)/np)); + for (int i = 2; i <= np / 2; i++) { + lam[2 * i - 3] = lam[2 * i - 2] = 2.0 * np * (1.0 - 1.0 * cos(2.0 * MY_PI * (i - 1) / np)); } // Set up eigenvectors for non-degenerated modes - for (int i=0; iiworld; - for (int i=1; i<=atom->ntypes; i++) - { + for (int i = 1; i <= atom->ntypes; i++) { mass[i] = atom->mass[i]; - if (iworld) - { + if (iworld) { mass[i] *= lam[iworld]; mass[i] *= fmass; } @@ -486,17 +507,17 @@ void FixPIMD::nmpimd_fill(double **ptr) /* ---------------------------------------------------------------------- */ -void FixPIMD::nmpimd_transform(double** src, double** des, double *vector) +void FixPIMD::nmpimd_transform(double **src, double **des, double *vector) { int n = atom->nlocal; int m = 0; - for (int i=0; ix; double **f = atom->f; - double* _mass = atom->mass; - int* type = atom->type; + double *_mass = atom->mass; + int *type = atom->type; int nlocal = atom->nlocal; - double* xlast = buf_beads[x_last]; - double* xnext = buf_beads[x_next]; + double *xlast = buf_beads[x_last]; + double *xnext = buf_beads[x_next]; - for (int i=0; ime - comm->nprocs; int rank_next = universe->me + comm->nprocs; - if (rank_last<0) rank_last += universe->nprocs; - if (rank_next>=universe->nprocs) rank_next -= universe->nprocs; + if (rank_last < 0) rank_last += universe->nprocs; + if (rank_next >= universe->nprocs) rank_next -= universe->nprocs; - plan_send[0] = rank_next; plan_send[1] = rank_last; - plan_recv[0] = rank_last; plan_recv[1] = rank_next; + plan_send[0] = rank_next; + plan_send[1] = rank_last; + plan_recv[0] = rank_last; + plan_recv[1] = rank_next; - mode_index[0] = 0; mode_index[1] = 1; - x_last = 1; x_next = 0; - } - else - { + mode_index[0] = 0; + mode_index[1] = 1; + x_last = 1; + x_next = 0; + } else { size_plan = np - 1; - plan_send = new int [size_plan]; - plan_recv = new int [size_plan]; - mode_index = new int [size_plan]; + plan_send = new int[size_plan]; + plan_recv = new int[size_plan]; + mode_index = new int[size_plan]; - for (int i=0; ime + comm->nprocs * (i+1); - if (plan_send[i]>=universe->nprocs) plan_send[i] -= universe->nprocs; + for (int i = 0; i < size_plan; i++) { + plan_send[i] = universe->me + comm->nprocs * (i + 1); + if (plan_send[i] >= universe->nprocs) plan_send[i] -= universe->nprocs; - plan_recv[i] = universe->me - comm->nprocs * (i+1); - if (plan_recv[i]<0) plan_recv[i] += universe->nprocs; + plan_recv[i] = universe->me - comm->nprocs * (i + 1); + if (plan_recv[i] < 0) plan_recv[i] += universe->nprocs; - mode_index[i]=(universe->iworld+i+1)%(universe->nworlds); + mode_index[i] = (universe->iworld + i + 1) % (universe->nworlds); } - x_next = (universe->iworld+1+universe->nworlds)%(universe->nworlds); - x_last = (universe->iworld-1+universe->nworlds)%(universe->nworlds); + x_next = (universe->iworld + 1 + universe->nworlds) % (universe->nworlds); + x_last = (universe->iworld - 1 + universe->nworlds) % (universe->nworlds); } - if (buf_beads) - { - for (int i=0; inlocal; - if (nlocal > max_nlocal) - { - max_nlocal = nlocal+200; + if (nlocal > max_nlocal) { + max_nlocal = nlocal + 200; int size = sizeof(double) * max_nlocal * 3; - buf_recv = (double*) memory->srealloc(buf_recv, size, "FixPIMD:x_recv"); + buf_recv = (double *) memory->srealloc(buf_recv, size, "FixPIMD:x_recv"); - for (int i=0; isrealloc(buf_beads[i], size, "FixPIMD:x_beads[i]"); + for (int i = 0; i < np; i++) + buf_beads[i] = (double *) memory->srealloc(buf_beads[i], size, "FixPIMD:x_beads[i]"); } // copy local positions - memcpy(buf_beads[universe->iworld], &(ptr[0][0]), sizeof(double)*nlocal*3); + memcpy(buf_beads[universe->iworld], &(ptr[0][0]), sizeof(double) * nlocal * 3); // go over comm plans - for (int iplan = 0; iplanuworld, MPI_STATUS_IGNORE); + MPI_Sendrecv(&(nlocal), 1, MPI_INT, plan_send[iplan], 0, &(nsend), 1, MPI_INT, plan_recv[iplan], + 0, universe->uworld, MPI_STATUS_IGNORE); // allocate arrays - if (nsend > max_nsend) - { - max_nsend = nsend+200; - tag_send = (tagint*) memory->srealloc(tag_send, sizeof(tagint)*max_nsend, "FixPIMD:tag_send"); - buf_send = (double*) memory->srealloc(buf_send, sizeof(double)*max_nsend*3, "FixPIMD:x_send"); + if (nsend > max_nsend) { + max_nsend = nsend + 200; + tag_send = + (tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMD:tag_send"); + buf_send = + (double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3, "FixPIMD:x_send"); } // send tags - MPI_Sendrecv( atom->tag, nlocal, MPI_LMP_TAGINT, plan_send[iplan], 0, - tag_send, nsend, MPI_LMP_TAGINT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE); + MPI_Sendrecv(atom->tag, nlocal, MPI_LMP_TAGINT, plan_send[iplan], 0, tag_send, nsend, + MPI_LMP_TAGINT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE); // wrap positions double *wrap_ptr = buf_send; - int ncpy = sizeof(double)*3; + int ncpy = sizeof(double) * 3; - for (int i=0; imap(tag_send[i]); - if (index<0) - { + if (index < 0) { char error_line[256]; - sprintf(error_line, "Atom " TAGINT_FORMAT " is missing at world [%d] " - "rank [%d] required by rank [%d] (" TAGINT_FORMAT ", " - TAGINT_FORMAT ", " TAGINT_FORMAT ").\n", tag_send[i], - universe->iworld, comm->me, plan_recv[iplan], - atom->tag[0], atom->tag[1], atom->tag[2]); + sprintf(error_line, + "Atom " TAGINT_FORMAT " is missing at world [%d] " + "rank [%d] required by rank [%d] (" TAGINT_FORMAT ", " TAGINT_FORMAT + ", " TAGINT_FORMAT ").\n", + tag_send[i], universe->iworld, comm->me, plan_recv[iplan], atom->tag[0], + atom->tag[1], atom->tag[2]); - error->universe_one(FLERR,error_line); + error->universe_one(FLERR, error_line); } memcpy(wrap_ptr, ptr[index], ncpy); @@ -677,21 +694,20 @@ void FixPIMD::comm_exec(double **ptr) // sendrecv x - MPI_Sendrecv( buf_send, nsend*3, MPI_DOUBLE, plan_recv[iplan], 0, - buf_recv, nlocal*3, MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, MPI_STATUS_IGNORE); + MPI_Sendrecv(buf_send, nsend * 3, MPI_DOUBLE, plan_recv[iplan], 0, buf_recv, nlocal * 3, + MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, MPI_STATUS_IGNORE); // copy x - memcpy(buf_beads[mode_index[iplan]], buf_recv, sizeof(double)*nlocal*3); + memcpy(buf_beads[mode_index[iplan]], buf_recv, sizeof(double) * nlocal * 3); } } /* ---------------------------------------------------------------------- */ -int FixPIMD::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int FixPIMD::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - int i,j,m; + int i, j, m; m = 0; @@ -709,7 +725,7 @@ int FixPIMD::pack_forward_comm(int n, int *list, double *buf, void FixPIMD::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -726,47 +742,51 @@ void FixPIMD::unpack_forward_comm(int n, int first, double *buf) double FixPIMD::memory_usage() { - return (double)atom->nmax * size_peratom_cols * sizeof(double); + return (double) atom->nmax * size_peratom_cols * sizeof(double); } /* ---------------------------------------------------------------------- */ void FixPIMD::grow_arrays(int nmax) { - if (nmax==0) return; - int count = nmax*3; + if (nmax == 0) return; + int count = nmax * 3; memory->grow(array_atom, nmax, size_peratom_cols, "FixPIMD::array_atom"); - memory->grow(nhc_eta, count, nhc_nchain, "FixPIMD::nh_eta"); - memory->grow(nhc_eta_dot, count, nhc_nchain+1, "FixPIMD::nh_eta_dot"); - memory->grow(nhc_eta_dotdot, count, nhc_nchain, "FixPIMD::nh_eta_dotdot"); - memory->grow(nhc_eta_mass, count, nhc_nchain, "FixPIMD::nh_eta_mass"); + memory->grow(nhc_eta, count, nhc_nchain, "FixPIMD::nh_eta"); + memory->grow(nhc_eta_dot, count, nhc_nchain + 1, "FixPIMD::nh_eta_dot"); + memory->grow(nhc_eta_dotdot, count, nhc_nchain, "FixPIMD::nh_eta_dotdot"); + memory->grow(nhc_eta_mass, count, nhc_nchain, "FixPIMD::nh_eta_mass"); } /* ---------------------------------------------------------------------- */ void FixPIMD::copy_arrays(int i, int j, int /*delflag*/) { - int i_pos = i*3; - int j_pos = j*3; + int i_pos = i * 3; + int j_pos = j * 3; - memcpy(nhc_eta [j_pos], nhc_eta [i_pos], nhc_size_one_1); - memcpy(nhc_eta_dot [j_pos], nhc_eta_dot [i_pos], nhc_size_one_2); + memcpy(nhc_eta[j_pos], nhc_eta[i_pos], nhc_size_one_1); + memcpy(nhc_eta_dot[j_pos], nhc_eta_dot[i_pos], nhc_size_one_2); memcpy(nhc_eta_dotdot[j_pos], nhc_eta_dotdot[i_pos], nhc_size_one_1); - memcpy(nhc_eta_mass [j_pos], nhc_eta_mass [i_pos], nhc_size_one_1); + memcpy(nhc_eta_mass[j_pos], nhc_eta_mass[i_pos], nhc_size_one_1); } /* ---------------------------------------------------------------------- */ int FixPIMD::pack_exchange(int i, double *buf) { - int offset=0; + int offset = 0; int pos = i * 3; - memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1; - memcpy(buf+offset, nhc_eta_dot[pos], nhc_size_one_2); offset += nhc_offset_one_2; - memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1; - memcpy(buf+offset, nhc_eta_mass[pos], nhc_size_one_1); offset += nhc_offset_one_1; + memcpy(buf + offset, nhc_eta[pos], nhc_size_one_1); + offset += nhc_offset_one_1; + memcpy(buf + offset, nhc_eta_dot[pos], nhc_size_one_2); + offset += nhc_offset_one_2; + memcpy(buf + offset, nhc_eta_dotdot[pos], nhc_size_one_1); + offset += nhc_offset_one_1; + memcpy(buf + offset, nhc_eta_mass[pos], nhc_size_one_1); + offset += nhc_offset_one_1; return size_peratom_cols; } @@ -775,13 +795,17 @@ int FixPIMD::pack_exchange(int i, double *buf) int FixPIMD::unpack_exchange(int nlocal, double *buf) { - int offset=0; - int pos = nlocal*3; + int offset = 0; + int pos = nlocal * 3; - memcpy(nhc_eta[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1; - memcpy(nhc_eta_dot[pos], buf+offset, nhc_size_one_2); offset += nhc_offset_one_2; - memcpy(nhc_eta_dotdot[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1; - memcpy(nhc_eta_mass[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1; + memcpy(nhc_eta[pos], buf + offset, nhc_size_one_1); + offset += nhc_offset_one_1; + memcpy(nhc_eta_dot[pos], buf + offset, nhc_size_one_2); + offset += nhc_offset_one_2; + memcpy(nhc_eta_dotdot[pos], buf + offset, nhc_size_one_1); + offset += nhc_offset_one_1; + memcpy(nhc_eta_mass[pos], buf + offset, nhc_size_one_1); + offset += nhc_offset_one_1; return size_peratom_cols; } @@ -790,17 +814,21 @@ int FixPIMD::unpack_exchange(int nlocal, double *buf) int FixPIMD::pack_restart(int i, double *buf) { - int offset=0; + int offset = 0; int pos = i * 3; // pack buf[0] this way because other fixes unpack it - buf[offset++] = size_peratom_cols+1; + buf[offset++] = size_peratom_cols + 1; - memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1; - memcpy(buf+offset, nhc_eta_dot[pos], nhc_size_one_2); offset += nhc_offset_one_2; - memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1; - memcpy(buf+offset, nhc_eta_mass[pos], nhc_size_one_1); offset += nhc_offset_one_1; + memcpy(buf + offset, nhc_eta[pos], nhc_size_one_1); + offset += nhc_offset_one_1; + memcpy(buf + offset, nhc_eta_dot[pos], nhc_size_one_2); + offset += nhc_offset_one_2; + memcpy(buf + offset, nhc_eta_dotdot[pos], nhc_size_one_1); + offset += nhc_offset_one_1; + memcpy(buf + offset, nhc_eta_mass[pos], nhc_size_one_1); + offset += nhc_offset_one_1; - return size_peratom_cols+1; + return size_peratom_cols + 1; } /* ---------------------------------------------------------------------- */ @@ -813,15 +841,19 @@ void FixPIMD::unpack_restart(int nlocal, int nth) // unpack the Nth first values this way because other fixes pack them int m = 0; - for (int i=0; i (extra[nlocal][m]); + for (int i = 0; i < nth; i++) m += static_cast(extra[nlocal][m]); m++; int pos = nlocal * 3; - memcpy(nhc_eta[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1; - memcpy(nhc_eta_dot[pos], extra[nlocal]+m, nhc_size_one_2); m += nhc_offset_one_2; - memcpy(nhc_eta_dotdot[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1; - memcpy(nhc_eta_mass[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1; + memcpy(nhc_eta[pos], extra[nlocal] + m, nhc_size_one_1); + m += nhc_offset_one_1; + memcpy(nhc_eta_dot[pos], extra[nlocal] + m, nhc_size_one_2); + m += nhc_offset_one_2; + memcpy(nhc_eta_dotdot[pos], extra[nlocal] + m, nhc_size_one_1); + m += nhc_offset_one_1; + memcpy(nhc_eta_mass[pos], extra[nlocal] + m, nhc_size_one_1); + m += nhc_offset_one_1; nhc_ready = true; } @@ -830,21 +862,21 @@ void FixPIMD::unpack_restart(int nlocal, int nth) int FixPIMD::maxsize_restart() { - return size_peratom_cols+1; + return size_peratom_cols + 1; } /* ---------------------------------------------------------------------- */ int FixPIMD::size_restart(int /*nlocal*/) { - return size_peratom_cols+1; + return size_peratom_cols + 1; } /* ---------------------------------------------------------------------- */ double FixPIMD::compute_vector(int n) { - if (n==0) { return spring_energy; } - if (n==1) { return t_sys; } + if (n == 0) { return spring_energy; } + if (n == 1) { return t_sys; } return 0.0; } diff --git a/src/REPLICA/fix_pimd.h b/src/REPLICA/fix_pimd.h index 2c13663f23..05a24dedce 100644 --- a/src/REPLICA/fix_pimd.h +++ b/src/REPLICA/fix_pimd.h @@ -27,6 +27,7 @@ namespace LAMMPS_NS { class FixPIMD : public Fix { public: FixPIMD(class LAMMPS *, int, char **); + virtual ~FixPIMD(); int setmask(); diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 66a1dd83d5..a445ace756 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -94,7 +94,8 @@ NEB::~NEB() { MPI_Comm_free(&roots); memory->destroy(all); - delete [] rdist; + delete[] rdist; + if (fp) fclose(fp); } /* ---------------------------------------------------------------------- @@ -379,8 +380,8 @@ void NEB::readfile(char *file, int flag) char line[MAXLINE]; double xx,yy,zz,delx,dely,delz; - if (me_universe == 0 && screen) - fprintf(screen,"Reading NEB coordinate file(s) ...\n"); + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Reading NEB coordinate file(s) ...\n"); // flag = 0, universe root reads header of file, bcast to universe // flag = 1, each replica's root reads header of file, bcast to world @@ -534,6 +535,7 @@ void NEB::readfile(char *file, int flag) else fclose(fp); } } + fp = nullptr; } /* ---------------------------------------------------------------------- diff --git a/src/SPIN/neb_spin.cpp b/src/SPIN/neb_spin.cpp index a07a16350a..1b798166a1 100644 --- a/src/SPIN/neb_spin.cpp +++ b/src/SPIN/neb_spin.cpp @@ -78,7 +78,8 @@ NEBSpin::~NEBSpin() { MPI_Comm_free(&roots); memory->destroy(all); - delete [] rdist; + delete[] rdist; + if (fp) fclose(fp); } /* ---------------------------------------------------------------------- @@ -373,8 +374,8 @@ void NEBSpin::readfile(char *file, int flag) double xx,yy,zz; double musp,spx,spy,spz; - if (me_universe == 0 && screen) - fprintf(screen,"Reading NEBSpin coordinate file(s) ...\n"); + if (me_universe == 0 && universe->uscreen) + fprintf(universe->uscreen,"Reading NEBSpin coordinate file(s) ...\n"); // flag = 0, universe root reads header of file, bcast to universe // flag = 1, each replica's root reads header of file, bcast to world @@ -546,8 +547,8 @@ void NEBSpin::readfile(char *file, int flag) // clean up - delete [] buffer; - delete [] values; + delete[] buffer; + delete[] values; if (flag == 0) { if (me_universe == 0) { @@ -560,6 +561,7 @@ void NEBSpin::readfile(char *file, int flag) else fclose(fp); } } + fp = nullptr; } /* ---------------------------------------------------------------------- diff --git a/src/arg_info.cpp b/src/arg_info.cpp index a274b47b92..5063b6b471 100644 --- a/src/arg_info.cpp +++ b/src/arg_info.cpp @@ -30,7 +30,8 @@ using namespace LAMMPS_NS; ArgInfo::ArgInfo(const std::string &arg, int allowed) : type(NONE), dim(0), index1(-1), index2(-1) { - if ((arg.size() > 2) && (arg[1] == '_')) { + if (((arg.size() > 3) && (arg[1] == '2') && (arg[2] == '_')) + || ((arg.size() > 2) && (arg[1] == '_'))) { if ((arg[0] == 'c') && (allowed & COMPUTE)) type = COMPUTE; else if ((arg[0] == 'f') && (allowed & FIX)) @@ -46,10 +47,11 @@ ArgInfo::ArgInfo(const std::string &arg, int allowed) : type(NONE), dim(0), inde name = arg; return; } + const int offset = (arg[1] == '_') ? 2 : 3; - std::size_t has_idx1 = arg.find('[', 2); + std::size_t has_idx1 = arg.find('[', offset); if (has_idx1 != std::string::npos) { - name = arg.substr(2, has_idx1 - 2); + name = arg.substr(offset, has_idx1 - offset); dim = 1; std::size_t has_idx2 = arg.find('[', has_idx1 + 1); @@ -79,7 +81,7 @@ ArgInfo::ArgInfo(const std::string &arg, int allowed) : type(NONE), dim(0), inde } } else { index1 = 0; - name = arg.substr(2); + name = arg.substr(offset); } } else { index1 = 0; diff --git a/src/compute_dipole.cpp b/src/compute_dipole.cpp index 053f11e6f9..b67c3df3ce 100644 --- a/src/compute_dipole.cpp +++ b/src/compute_dipole.cpp @@ -29,7 +29,7 @@ enum { MASSCENTER, GEOMCENTER }; ComputeDipole::ComputeDipole(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg != 3) error->all(FLERR, "Illegal compute com command"); + if ((narg < 3) || (narg > 4)) error->all(FLERR, "Illegal compute dipole command"); scalar_flag = 1; vector_flag = 1; @@ -39,11 +39,8 @@ ComputeDipole::ComputeDipole(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, n vector = new double[size_vector]; vector[0] = vector[1] = vector[2] = 0.0; - usecenter = MASSCENTER; - if ((narg != 3) && (narg != 4)) error->all(FLERR, "Illegal compute dipole command"); - if (narg == 4) { if (utils::strmatch(arg[3], "^geom")) usecenter = GEOMCENTER; diff --git a/src/read_restart.cpp b/src/read_restart.cpp index faf9fd1620..e42bb59272 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -44,7 +44,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -ReadRestart::ReadRestart(LAMMPS *lmp) : Command(lmp) {} +ReadRestart::ReadRestart(LAMMPS *lmp) : Command(lmp), mpiio(nullptr) {} /* ---------------------------------------------------------------------- */ @@ -110,8 +110,7 @@ void ReadRestart::command(int narg, char **arg) } fp = fopen(hfile.c_str(),"rb"); if (fp == nullptr) - error->one(FLERR,"Cannot open restart file {}: {}", - hfile, utils::getsyserror()); + error->one(FLERR,"Cannot open restart file {}: {}", hfile, utils::getsyserror()); } // read magic string, endian flag, format revision @@ -336,8 +335,7 @@ void ReadRestart::command(int narg, char **arg) procfile.replace(procfile.find("%"),1,fmt::format("{}",icluster)); fp = fopen(procfile.c_str(),"rb"); if (fp == nullptr) - error->one(FLERR,"Cannot open restart file {}: {}", - procfile, utils::getsyserror()); + error->one(FLERR,"Cannot open restart file {}: {}", procfile, utils::getsyserror()); } int procsperfile; @@ -401,7 +399,7 @@ void ReadRestart::command(int narg, char **arg) // clean-up memory - delete [] file; + delete[] file; memory->destroy(buf); // for multiproc or MPI-IO files: @@ -525,8 +523,9 @@ void ReadRestart::command(int narg, char **arg) MPI_Barrier(world); if (comm->me == 0) - utils::logmesg(lmp," read_restart CPU = {:.3f} seconds\n", - MPI_Wtime()-time1); + utils::logmesg(lmp," read_restart CPU = {:.3f} seconds\n",MPI_Wtime()-time1); + + delete mpiio; } /* ---------------------------------------------------------------------- @@ -599,7 +598,7 @@ void ReadRestart::header() if (me == 0) utils::logmesg(lmp," restart file = {}, LAMMPS = {}\n", version,lmp->version); - delete [] version; + delete[] version; // we have no forward compatibility, thus exit with error @@ -638,7 +637,7 @@ void ReadRestart::header() } else if (flag == UNITS) { char *style = read_string(); if (strcmp(style,update->unit_style) != 0) update->set_units(style); - delete [] style; + delete[] style; } else if (flag == NTIMESTEP) { update->ntimestep = read_bigint(); @@ -772,9 +771,9 @@ void ReadRestart::header() atom->create_avec(style,nargcopy,argcopy,1); if (comm->me ==0) utils::logmesg(lmp," restoring atom style {} from restart\n",style); - for (int i = 0; i < nargcopy; i++) delete [] argcopy[i]; - delete [] argcopy; - delete [] style; + for (int i = 0; i < nargcopy; i++) delete[] argcopy[i]; + delete[] argcopy; + delete[] style; } else if (flag == NATOMS) { atom->natoms = read_bigint(); @@ -889,7 +888,7 @@ void ReadRestart::type_arrays() double *mass = new double[atom->ntypes+1]; read_double_vec(atom->ntypes,&mass[1]); atom->set_mass(mass); - delete [] mass; + delete[] mass; } else error->all(FLERR, "Invalid flag in type arrays section of restart file"); @@ -910,7 +909,7 @@ void ReadRestart::force_fields() if (flag == PAIR) { style = read_string(); force->create_pair(style,1); - delete [] style; + delete[] style; if (comm->me ==0) utils::logmesg(lmp," restoring pair style {} from restart\n", force->pair_style); @@ -926,7 +925,7 @@ void ReadRestart::force_fields() } else if (flag == BOND) { style = read_string(); force->create_bond(style,1); - delete [] style; + delete[] style; if (comm->me ==0) utils::logmesg(lmp," restoring bond style {} from restart\n", force->bond_style); @@ -935,7 +934,7 @@ void ReadRestart::force_fields() } else if (flag == ANGLE) { style = read_string(); force->create_angle(style,1); - delete [] style; + delete[] style; if (comm->me ==0) utils::logmesg(lmp," restoring angle style {} from restart\n", force->angle_style); @@ -944,7 +943,7 @@ void ReadRestart::force_fields() } else if (flag == DIHEDRAL) { style = read_string(); force->create_dihedral(style,1); - delete [] style; + delete[] style; if (comm->me ==0) utils::logmesg(lmp," restoring dihedral style {} from restart\n", force->dihedral_style); @@ -953,7 +952,7 @@ void ReadRestart::force_fields() } else if (flag == IMPROPER) { style = read_string(); force->create_improper(style,1); - delete [] style; + delete[] style; if (comm->me ==0) utils::logmesg(lmp," restoring improper style {} from restart\n", force->improper_style); @@ -1117,7 +1116,7 @@ void ReadRestart::magic_string() MPI_Bcast(str,n,MPI_CHAR,0,world); if (strcmp(str,MAGIC_STRING) != 0) error->all(FLERR,"Invalid LAMMPS restart file"); - delete [] str; + delete[] str; } /* ---------------------------------------------------------------------- @@ -1164,7 +1163,7 @@ void ReadRestart::check_eof_magic() if (strcmp(str,MAGIC_STRING) != 0) error->all(FLERR,"Incomplete or corrupted LAMMPS restart file"); - delete [] str; + delete[] str; } /* ---------------------------------------------------------------------- diff --git a/src/utils.cpp b/src/utils.cpp index 3e27d97e5c..8da2efd03d 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -177,20 +177,23 @@ char *utils::fgets_trunc(char *buf, int size, FILE *fp) char dummy[MAXDUMMY]; char *ptr = fgets(buf, size, fp); - // EOF + // EOF? if (!ptr) return nullptr; int n = strlen(buf); - // line is shorter than buffer, append newline if needed, - if (n < size - 2) { + // check the string being read in: + // - if string is shorter than the buffer make sure it has a final newline and return + // - if string is exactly the size of the buffer and has a final newline return + // - otherwise truncate with final newline and read into dummy buffer until EOF or newline is found + if (n < size - 1) { if (buf[n - 1] != '\n') { buf[n] = '\n'; buf[n + 1] = '\0'; } return buf; - - // line fits exactly. overwrite last but one character. + } else if (buf[n - 1] == '\n') { + return buf; } else buf[size - 2] = '\n'; @@ -203,15 +206,15 @@ char *utils::fgets_trunc(char *buf, int size, FILE *fp) n = 0; } while (n == MAXDUMMY - 1 && ptr[MAXDUMMY - 1] != '\n'); - // return first chunk + // return truncated chunk return buf; } -#define MAXPATHLENBUF 1024 /* like fgets() but aborts with an error or EOF is encountered */ void utils::sfgets(const char *srcname, int srcline, char *s, int size, FILE *fp, const char *filename, Error *error) { + constexpr int MAXPATHLENBUF=1024; char *rv = fgets(s, size, fp); if (rv == nullptr) { // something went wrong char buf[MAXPATHLENBUF]; @@ -240,6 +243,7 @@ void utils::sfgets(const char *srcname, int srcline, char *s, int size, FILE *fp void utils::sfread(const char *srcname, int srcline, void *s, size_t size, size_t num, FILE *fp, const char *filename, Error *error) { + constexpr int MAXPATHLENBUF=1024; size_t rv = fread(s, size, num, fp); if (rv != num) { // something went wrong char buf[MAXPATHLENBUF]; @@ -1093,11 +1097,6 @@ bool utils::file_is_readable(const std::string &path) search current directory and the LAMMPS_POTENTIALS directory if specified ------------------------------------------------------------------------- */ -#if defined(_WIN32) -#define OS_PATH_VAR_SEP ";" -#else -#define OS_PATH_VAR_SEP ":" -#endif std::string utils::get_potential_file_path(const std::string &path) { @@ -1111,8 +1110,11 @@ std::string utils::get_potential_file_path(const std::string &path) const char *var = getenv("LAMMPS_POTENTIALS"); if (var != nullptr) { - Tokenizer dirs(var, OS_PATH_VAR_SEP); - +#if defined(_WIN32) + Tokenizer dirs(var, ";"); +#else + Tokenizer dirs(var, ":"); +#endif while (dirs.has_next()) { auto pot = utils::path_basename(filepath); auto dir = dirs.next(); @@ -1124,7 +1126,6 @@ std::string utils::get_potential_file_path(const std::string &path) } return ""; } -#undef OS_PATH_VAR_SEP /* ---------------------------------------------------------------------- read first line of potential file diff --git a/unittest/formats/test_file_operations.cpp b/unittest/formats/test_file_operations.cpp index af2023cebb..0f7362674c 100644 --- a/unittest/formats/test_file_operations.cpp +++ b/unittest/formats/test_file_operations.cpp @@ -74,7 +74,8 @@ protected: } }; -#define MAX_BUF_SIZE 128 +static constexpr int MAX_BUF_SIZE=128; + TEST_F(FileOperationsTest, safe_fgets) { char buf[MAX_BUF_SIZE]; @@ -110,7 +111,6 @@ TEST_F(FileOperationsTest, safe_fgets) fclose(fp); } -#define MAX_BUF_SIZE 128 TEST_F(FileOperationsTest, fgets_trunc) { char buf[MAX_BUF_SIZE]; @@ -119,13 +119,15 @@ TEST_F(FileOperationsTest, fgets_trunc) FILE *fp = fopen("safe_file_read_test.txt", "rb"); ASSERT_NE(fp, nullptr); + // read line shorter than buffer memset(buf, 0, MAX_BUF_SIZE); ptr = utils::fgets_trunc(buf, MAX_BUF_SIZE, fp); ASSERT_THAT(buf, StrEq("one line\n")); ASSERT_NE(ptr,nullptr); + // read line of exactly the buffer length memset(buf, 0, MAX_BUF_SIZE); - ptr = utils::fgets_trunc(buf, MAX_BUF_SIZE, fp); + ptr = utils::fgets_trunc(buf, sizeof("two_lines\n"), fp); ASSERT_THAT(buf, StrEq("two_lines\n")); ASSERT_NE(ptr,nullptr); @@ -172,7 +174,6 @@ TEST_F(FileOperationsTest, fgets_trunc) fclose(fp); } -#define MAX_BUF_SIZE 128 TEST_F(FileOperationsTest, safe_fread) { char buf[MAX_BUF_SIZE]; diff --git a/unittest/utils/test_argutils.cpp b/unittest/utils/test_argutils.cpp index a527df6fe7..c5ef19bc10 100644 --- a/unittest/utils/test_argutils.cpp +++ b/unittest/utils/test_argutils.cpp @@ -169,6 +169,16 @@ TEST(ArgInfo, variable2) ASSERT_THAT(arg.get_name(), StrEq("x")); } +TEST(ArgInfo, variable3) +{ + ArgInfo arg("v_x[11][5]"); + ASSERT_EQ(arg.get_dim(), 2); + ASSERT_EQ(arg.get_type(), ArgInfo::VARIABLE); + ASSERT_EQ(arg.get_index1(), 11); + ASSERT_EQ(arg.get_index2(), 5); + ASSERT_THAT(arg.get_name(), StrEq("x")); +} + TEST(ArgInfo, dname0) { ArgInfo arg("d_text", ArgInfo::DNAME); @@ -179,6 +189,36 @@ TEST(ArgInfo, dname0) ASSERT_THAT(arg.get_name(), StrEq("text")); } +TEST(ArgInfo, dname1) +{ + ArgInfo arg("d2_text", ArgInfo::DNAME | ArgInfo::INAME); + ASSERT_EQ(arg.get_dim(), 0); + ASSERT_EQ(arg.get_type(), ArgInfo::DNAME); + ASSERT_EQ(arg.get_index1(), 0); + ASSERT_EQ(arg.get_index2(), -1); + ASSERT_THAT(arg.get_name(), StrEq("text")); +} + +TEST(ArgInfo, dname2) +{ + ArgInfo arg("d2_text[11]", ArgInfo::DNAME); + ASSERT_EQ(arg.get_dim(), 1); + ASSERT_EQ(arg.get_type(), ArgInfo::DNAME); + ASSERT_EQ(arg.get_index1(), 11); + ASSERT_EQ(arg.get_index2(), -1); + ASSERT_THAT(arg.get_name(), StrEq("text")); +} + +TEST(ArgInfo, dname3) +{ + ArgInfo arg("d2_text[24][11]", ArgInfo::DNAME); + ASSERT_EQ(arg.get_dim(), 2); + ASSERT_EQ(arg.get_type(), ArgInfo::DNAME); + ASSERT_EQ(arg.get_index1(), 24); + ASSERT_EQ(arg.get_index2(), 11); + ASSERT_THAT(arg.get_name(), StrEq("text")); +} + TEST(ArgInfo, iname0) { ArgInfo arg("i_text", ArgInfo::INAME); @@ -189,6 +229,36 @@ TEST(ArgInfo, iname0) ASSERT_THAT(arg.get_name(), StrEq("text")); } +TEST(ArgInfo, iname1) +{ + ArgInfo arg("i2_text", ArgInfo::INAME); + ASSERT_EQ(arg.get_dim(), 0); + ASSERT_EQ(arg.get_type(), ArgInfo::INAME); + ASSERT_EQ(arg.get_index1(), 0); + ASSERT_EQ(arg.get_index2(), -1); + ASSERT_THAT(arg.get_name(), StrEq("text")); +} + +TEST(ArgInfo, iname2) +{ + ArgInfo arg("i2_text[2]", ArgInfo::INAME | ArgInfo::DNAME); + ASSERT_EQ(arg.get_dim(), 1); + ASSERT_EQ(arg.get_type(), ArgInfo::INAME); + ASSERT_EQ(arg.get_index1(), 2); + ASSERT_EQ(arg.get_index2(), -1); + ASSERT_THAT(arg.get_name(), StrEq("text")); +} + +TEST(ArgInfo, iname3) +{ + ArgInfo arg("i2_text[2][100]", ArgInfo::INAME | ArgInfo::DNAME); + ASSERT_EQ(arg.get_dim(), 2); + ASSERT_EQ(arg.get_type(), ArgInfo::INAME); + ASSERT_EQ(arg.get_index1(), 2); + ASSERT_EQ(arg.get_index2(), 100); + ASSERT_THAT(arg.get_name(), StrEq("text")); +} + TEST(ArgInfo, unsupported1) { ArgInfo arg("v_text[02][05]", ArgInfo::COMPUTE | ArgInfo::FIX);