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);

}

MPI_Allreduce(&eng,&total_eng,1,MPI_DOUBLE,MPI_SUM,world);

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!

-Sebastian