diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index 7008d81241..f7a87f6913 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -21,14 +21,17 @@ #include "domain.h" #include "force.h" #include "bond.h" +#include "math_extra.h" +#include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define DELTA 10000 +#define EPSILON 1.0e-12 -enum{DIST,ENG,FORCE}; +enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE}; /* ---------------------------------------------------------------------- */ @@ -42,6 +45,8 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Compute bond/local used when bonds are not allowed"); local_flag = 1; + comm_forward = 3; + nvalues = narg - 3; if (nvalues == 1) size_local_cols = 0; else size_local_cols = nvalues; @@ -51,16 +56,26 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = 0; for (int iarg = 3; iarg < narg; iarg++) { if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST; - else if (strcmp(arg[iarg],"eng") == 0) bstyle[nvalues++] = ENG; + else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT; else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE; + else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB; + else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT; + else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS; + else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA; + else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB; else error->all(FLERR,"Invalid keyword in compute bond/local command"); } // set singleflag if need to call bond->single() + // set velflag if compute any quantities based on velocities singleflag = 0; - for (int i = 0; i < nvalues; i++) - if (bstyle[i] != DIST) singleflag = 1; + ghostvelflag = 0; + for (int i = 0; i < nvalues; i++) { + if (bstyle[i] == ENGPOT || bstyle[i] == FORCE) singleflag = 1; + if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || + bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1; + } nmax = 0; } @@ -81,9 +96,17 @@ void ComputeBondLocal::init() if (force->bond == NULL) error->all(FLERR,"No bond style is defined for compute bond/local"); + // set ghostvelflag if need to acquire ghost atom velocities + + if (velflag && !comm->ghost_velocity) ghostvelflag = 1; + else ghostvelflag = 0; + // do initial memory allocation so that memory_usage() is correct + initflag = 1; ncount = compute_bonds(0); + initflag = 0; + if (ncount > nmax) reallocate(ncount); size_local_rows = ncount; } @@ -117,10 +140,22 @@ int ComputeBondLocal::compute_bonds(int flag) { int i,m,n,nb,atom1,atom2,imol,iatom,btype; tagint tagprev; - double delx,dely,delz,rsq; + double dx,dy,dz,rsq; + double mass1,mass2,masstotal,invmasstotal; + double xcm[3],vcm[3]; + double delr1[3],delr2[3],delv1[3],delv2[3]; + double r12[3],vpar1,vpar2; + double vvib,vrotsq; + double inertia,omegasq; + double mvv2e; + double engpot,engtrans,engvib,engrot,engtot,fbond; double *ptr; double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; tagint *tag = atom->tag; int *num_bond = atom->num_bond; tagint **bond_atom = atom->bond_atom; @@ -136,7 +171,12 @@ int ComputeBondLocal::compute_bonds(int flag) int molecular = atom->molecular; Bond *bond = force->bond; - double eng,fbond; + + // communicate ghost velocities if needed + + if (ghostvelflag && !initflag) comm->forward_comm_compute(this); + + // loop over all atoms and their bonds m = n = 0; for (atom1 = 0; atom1 < nlocal; atom1++) { @@ -164,16 +204,93 @@ int ComputeBondLocal::compute_bonds(int flag) if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; if (btype == 0) continue; - if (flag) { - delx = x[atom1][0] - x[atom2][0]; - dely = x[atom1][1] - x[atom2][1]; - delz = x[atom1][2] - x[atom2][2]; - domain->minimum_image(delx,dely,delz); - rsq = delx*delx + dely*dely + delz*delz; + if (!flag) { + m++; + continue; + } - if (singleflag && (btype > 0)) - eng = bond->single(btype,rsq,atom1,atom2,fbond); - else eng = fbond = 0.0; + dx = x[atom1][0] - x[atom2][0]; + dy = x[atom1][1] - x[atom2][1]; + dz = x[atom1][2] - x[atom2][2]; + domain->minimum_image(dx,dy,dz); + rsq = dx*dx + dy*dy + dz*dz; + + if (btype == 0) { + engpot = fbond = 0.0; + engvib = engrot = engtrans = omegasq = vvib = 0.0; + } else { + + if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond); + + if (velflag) { + if (rmass) { + mass1 = rmass[atom1]; + mass2 = rmass[atom2]; + } + else { + mass1 = mass[type[atom1]]; + mass2 = mass[type[atom2]]; + } + masstotal = mass1+mass2; + invmasstotal = 1.0 / (masstotal); + xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal; + xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal; + xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal; + vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal; + vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal; + vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal; + + engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm); + + // r12 = unit bond vector from atom1 to atom2 + + MathExtra::sub3(x[atom2],x[atom1],r12); + MathExtra::norm3(r12); + + // delr = vector from COM to each atom + // delv = velocity of each atom relative to COM + + MathExtra::sub3(x[atom1],xcm,delr1); + MathExtra::sub3(x[atom2],xcm,delr2); + MathExtra::sub3(v[atom1],vcm,delv1); + MathExtra::sub3(v[atom2],vcm,delv2); + + // vpar = component of delv parallel to bond vector + + vpar1 = MathExtra::dot3(delv1,r12); + vpar2 = MathExtra::dot3(delv2,r12); + engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2); + + // vvib = relative velocity of 2 atoms along bond direction + // vvib < 0 for 2 atoms moving towards each other + // vvib > 0 for 2 atoms moving apart + + vvib = vpar2 - vpar1; + + // vrotsq = tangential speed squared of atom1 only + // omegasq = omega squared, and is the same for atom1 and atom2 + + inertia = mass1*MathExtra::lensq3(delr1) + + mass2*MathExtra::lensq3(delr2); + vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1; + omegasq = vrotsq / MathExtra::lensq3(delr1); + + engrot = 0.5 * inertia * omegasq; + + // sanity check: engtotal = engtrans + engvib + engrot + + //engtot = 0.5 * (mass1*MathExtra::lensq3(v[atom1]) + + // mass2*MathExtra::lensq3(v[atom2])); + //if (fabs(engtot-engtrans-engvib-engrot) > EPSILON) + // error->one(FLERR,"Sanity check on 3 energy components failed"); + + // scale energies by units + + mvv2e = force->mvv2e; + engtrans *= mvv2e; + engvib *= mvv2e; + engrot *= mvv2e; + } if (nvalues == 1) ptr = &vector[m]; else ptr = array[m]; @@ -183,12 +300,27 @@ int ComputeBondLocal::compute_bonds(int flag) case DIST: ptr[n] = sqrt(rsq); break; - case ENG: - ptr[n] = eng; + case ENGPOT: + ptr[n] = engpot; break; case FORCE: ptr[n] = sqrt(rsq)*fbond; break; + case ENGVIB: + ptr[n] = engvib; + break; + case ENGROT: + ptr[n] = engrot; + break; + case ENGTRANS: + ptr[n] = engtrans; + break; + case OMEGA: + ptr[n] = sqrt(omegasq); + break; + case VELVIB: + ptr[n] = vvib; + break; } } } @@ -202,6 +334,43 @@ int ComputeBondLocal::compute_bonds(int flag) /* ---------------------------------------------------------------------- */ +int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + double **v = atom->v; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf) +{ + int i,m,last; + + double **v = atom->v; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputeBondLocal::reallocate(int n) { // grow vector or array and indices array diff --git a/src/compute_bond_local.h b/src/compute_bond_local.h index 7d04515007..64cafc75e0 100644 --- a/src/compute_bond_local.h +++ b/src/compute_bond_local.h @@ -30,13 +30,15 @@ class ComputeBondLocal : public Compute { ~ComputeBondLocal(); void init(); void compute_local(); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); double memory_usage(); private: int nvalues; int ncount; int *bstyle; - int singleflag; + int singleflag,velflag,ghostvelflag,initflag; int nmax; diff --git a/src/dump_local.cpp b/src/dump_local.cpp index 3586a16518..7e1381b05f 100644 --- a/src/dump_local.cpp +++ b/src/dump_local.cpp @@ -21,6 +21,7 @@ #include "compute.h" #include "domain.h" #include "update.h" +#include "input.h" #include "memory.h" #include "error.h" #include "force.h" @@ -45,7 +46,18 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) : nevery = force->inumeric(FLERR,arg[3]); if (nevery <= 0) error->all(FLERR,"Illegal dump local command"); - size_one = nfield = narg-5; + nfield = narg - 5; + + // expand args if any have wildcard character "*" + + int expand = 0; + char **earg; + nfield = input->expand_args(nfield,&arg[5],1,earg); + + if (earg != &arg[5]) expand = 1; + + // allocate field vectors + pack_choice = new FnPtrPack[nfield]; vtype = new int[nfield]; @@ -67,7 +79,8 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) : // process attributes - parse_fields(narg,arg); + parse_fields(nfield,earg); + size_one = nfield; // setup format strings @@ -88,11 +101,11 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) : // setup column string int n = 0; - for (int iarg = 5; iarg < narg; iarg++) n += strlen(arg[iarg]) + 2; + for (int iarg = 0; iarg < nfield; iarg++) n += strlen(earg[iarg]) + 2; columns = new char[n]; columns[0] = '\0'; - for (int iarg = 5; iarg < narg; iarg++) { - strcat(columns,arg[iarg]); + for (int iarg = 0; iarg < nfield; iarg++) { + strcat(columns,earg[iarg]); strcat(columns," "); } @@ -102,6 +115,13 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) : n = strlen(str) + 1; label = new char[n]; strcpy(label,str); + + // if wildcard expansion occurred, free earg memory from exapnd_args() + + if (expand) { + for (int i = 0; i < nfield; i++) delete [] earg[i]; + memory->sfree(earg); + } } /* ---------------------------------------------------------------------- */ @@ -376,8 +396,8 @@ void DumpLocal::parse_fields(int narg, char **arg) // customize by adding to if statement int i; - for (int iarg = 5; iarg < narg; iarg++) { - i = iarg-5; + for (int iarg = 0; iarg < narg; iarg++) { + i = iarg; if (strcmp(arg[iarg],"index") == 0) { pack_choice[i] = &DumpLocal::pack_index; diff --git a/src/fix_ave_atom.cpp b/src/fix_ave_atom.cpp index ce904df85d..cf2842f441 100644 --- a/src/fix_ave_atom.cpp +++ b/src/fix_ave_atom.cpp @@ -36,7 +36,8 @@ enum{X,V,F,COMPUTE,FIX,VARIABLE}; FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - nvalues(0), which(NULL), argindex(NULL), value2index(NULL), ids(NULL), array(NULL) + nvalues(0), which(NULL), argindex(NULL), value2index(NULL), + ids(NULL), array(NULL) { if (narg < 7) error->all(FLERR,"Illegal fix ave/atom command");