diff --git a/src/fix_add_force.cpp b/src/fix_add_force.cpp index 4ac65b8bfb..252b450b63 100644 --- a/src/fix_add_force.cpp +++ b/src/fix_add_force.cpp @@ -28,9 +28,11 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) : { if (narg != 6) error->all("Illegal fix addforce command"); + scalar_flag = 1; vector_flag = 1; size_vector = 3; scalar_vector_freq = 1; + extscalar = 1; extvector = 1; xvalue = atof(arg[3]); @@ -38,7 +40,7 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) : zvalue = atof(arg[5]); force_flag = 0; - foriginal[0] = foriginal[1] = foriginal[2] = 0.0; + foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0; } /* ---------------------------------------------------------------------- */ @@ -47,6 +49,7 @@ int FixAddForce::setmask() { int mask = 0; mask |= POST_FORCE; + mask |= THERMO_ENERGY; mask |= POST_FORCE_RESPA; mask |= MIN_POST_FORCE; return mask; @@ -84,18 +87,22 @@ void FixAddForce::min_setup(int vflag) void FixAddForce::post_force(int vflag) { + double **x = atom->x; double **f = atom->f; int *mask = atom->mask; int nlocal = atom->nlocal; - foriginal[0] = foriginal[1] = foriginal[2] = 0.0; + foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0; force_flag = 0; + // foriginal[0] = - x dot f = "potential" for added force + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - foriginal[0] += f[i][0]; - foriginal[1] += f[i][1]; - foriginal[2] += f[i][2]; + foriginal[0] -= xvalue*x[i][0] + yvalue*x[i][1] + zvalue[i][2]; + foriginal[1] += f[i][0]; + foriginal[2] += f[i][1]; + foriginal[3] += f[i][2]; f[i][0] += xvalue; f[i][1] += yvalue; f[i][2] += zvalue; @@ -116,6 +123,21 @@ void FixAddForce::min_post_force(int vflag) post_force(vflag); } +/* ---------------------------------------------------------------------- + potential energy of added force +------------------------------------------------------------------------- */ + +double FixAddForce::compute_scalar() +{ + // only sum across procs one time + + if (force_flag == 0) { + MPI_Allreduce(foriginal,foriginal_all,4,MPI_DOUBLE,MPI_SUM,world); + force_flag = 1; + } + return foriginal_all[0]; +} + /* ---------------------------------------------------------------------- return components of total force on fix group before force was changed ------------------------------------------------------------------------- */ @@ -125,8 +147,8 @@ double FixAddForce::compute_vector(int n) // only sum across procs one time if (force_flag == 0) { - MPI_Allreduce(foriginal,foriginal_all,3,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(foriginal,foriginal_all,4,MPI_DOUBLE,MPI_SUM,world); force_flag = 1; } - return foriginal_all[n]; + return foriginal_all[n+1]; } diff --git a/src/fix_add_force.h b/src/fix_add_force.h index 7b5bd2ec31..52be0c55f8 100644 --- a/src/fix_add_force.h +++ b/src/fix_add_force.h @@ -28,11 +28,12 @@ class FixAddForce : public Fix { void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); + double compute_scalar(); double compute_vector(int); private: double xvalue,yvalue,zvalue; - double foriginal[3],foriginal_all[3]; + double foriginal[4],foriginal_all[4]; int force_flag; int nlevels_respa; }; diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp index 457f9dbcc3..65c2d41b11 100644 --- a/src/fix_spring_self.cpp +++ b/src/fix_spring_self.cpp @@ -94,6 +94,7 @@ int FixSpringSelf::setmask() { int mask = 0; mask |= POST_FORCE; + mask |= THERMO_ENERGY; mask |= POST_FORCE_RESPA; mask |= MIN_POST_FORCE; return mask; @@ -122,6 +123,13 @@ void FixSpringSelf::setup(int vflag) /* ---------------------------------------------------------------------- */ +void FixSpringSelf::min_setup(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + void FixSpringSelf::post_force(int vflag) { double **x = atom->x; @@ -135,7 +143,7 @@ void FixSpringSelf::post_force(int vflag) double zprd = domain->zprd; int xbox,ybox,zbox; double dx,dy,dz; - double espring = 0.0; + espring = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -150,6 +158,8 @@ void FixSpringSelf::post_force(int vflag) f[i][2] -= k*dz; espring += k * (dx*dx + dy*dy + dz*dz); } + + espring *= 0.5; } /* ---------------------------------------------------------------------- */ @@ -165,7 +175,6 @@ void FixSpringSelf::post_force_respa(int vflag, int ilevel, int iloop) double FixSpringSelf::compute_scalar() { - espring *= 0.5; double all; MPI_Allreduce(&espring,&all,1,MPI_DOUBLE,MPI_SUM,world); return all; diff --git a/src/fix_spring_self.h b/src/fix_spring_self.h index d217f0a89b..64dd01cf56 100644 --- a/src/fix_spring_self.h +++ b/src/fix_spring_self.h @@ -25,6 +25,7 @@ class FixSpringSelf : public Fix { int setmask(); void init(); void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int);