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