[lammps-users] Fix Spring/Self and Energy minimization, possible Bug

Hi everybody,

I would like to use energy minimization with some atoms fixed to their
original position by harmonic restraints. For this I figured, the energy of
the springs has to be included in the potential energy and the manual suggest
that this can be done with "fix modify energy yes".
However after a few attempts I looked into the code and found that fix
spring/self had no compute_scalar() function, which I understand is used by
the potential energy compute for the energy contribution of fix.

So I added to fix_spring_self.cpp

double FixSpringSelf::compute_scalar()

    double **x = atom->x;
    int *mask = atom->mask;
    int *image = atom->image;
    int nlocal = atom->nlocal;

    double xprd = domain->xprd;
    double yprd = domain->yprd;
    double zprd = domain->zprd;
    int xbox,ybox,zbox;
   double dx,dy,dz;
   double eng = 0;
   double total_eng;
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) {
        xbox = (image[i] & 1023) - 512;
        ybox = (image[i] >> 10 & 1023) - 512;
        zbox = (image[i] >> 20) - 512;
        dx = x[i][0] + xbox*xprd - xoriginal[i][0];
        dy = x[i][1] + ybox*yprd - xoriginal[i][1];
        dz = x[i][2] + zbox*zprd - xoriginal[i][2];
        eng += k/2*(dx*dx+dy*dy+dz*dz);
    return total_eng;
FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
  scalar_flag = 1;
  scalar_vector_freq = 1;

Did I do everthing correctly? Are the forces from the springs taken into
account for the minimization, I am not sure I fully understand the
minimization code...
For my system the difference between the modified version and the original is
that minimization converges after a lower number of steps at a 10% higher
energy. The displaced atoms are indeed at a lower msqd from the original
postition, but the rate with which the msqd changes at every minimization
step is the same.

Tank you!