Re-enable colvars from binary file using read_restart

LAMMPS LAMMPS Mailing List Mirror
1

LAMMPS version 11Aug17

Hi all,

We are currently using LAMMPS together with the USER-COLVARS package. We are not interested in computing the energy of the colvars, but to enforce the proximity of many particle pairs using harmonics. So, stretching the usual applications of the COLVARS package, we wish to load many colvars (~10,000 if possible). However to load all the variables it takes a lot of time, sometimes even more than the simulation itself.

Is there the possibility to speed the loading process?

We were thinking to try to load them from the binary restart file, because eventually the loading process may be faster. However, checking the mailing list or the manual, we couldn’t find a way to re-enable the colvars once the restart file has been loaded. Can you help us on this?

Thanks for all in advance,

Julen

2

LAMMPS version 11Aug17

please update to a newer version. there have been numerous updates and bugfixes since.

colvars can be restarted (at least in the current version), but that won’t help you much because what the restart procedure in fix colvars does is the do the same serialization that happens when writing to a colvars restart file, but to a string buffer and then store this string buffer. so you won’t gain anything. you still need to read/parse a text file and that cannot be changed with colvars. in general, since colvars is not parallelized, you also have a significant load imbalance caused by colvars.

what keeps you from using explicit harmonic bonds?

axel.

3

I second what Axel suggests. This type of restraint goes a bit outside the scope of Colvars, which assumes a high degree of dimensionality reduction.

If there is any reason why harmonic bonds wouldn’t work well (mind that there are lots of things you can do in LAMMPS!) you could look into the distancePairs variable:
https://colvars.github.io/colvars-refman-lammps/colvars-refman-lammps.html#sec:cvc_distancePairs
and apply a 10,000-dimensional harmonic restraint to it. This would also not help with the parallel performance but most definitely with the initialization time (a single variable instead of many).

Again, this would also get really tedious, because you would need to write “centers (x0, y0, z0, …)” as a single line in the restraint… Use this route only if regular bonds really cannot be used.

Giacomo

4

Hi Axel and Giacomo,

Originally we used colvars for another project in which there was time dependence, but luckily that isn’t the case in mine. We’ll see how to include the harmonic bounds instead of colvars and come back to you if we encounter any issue that we cannot solve with the manual.

Thanks a lot,
Julen

5

Hi again,

Sorry for coming back so soon, but we use a lot of lowerbound harmonics, do they exist in LAMMPS, and if not what is the best way to implement them?

Thanks in advance,
Julen

6

sorry, but i don’t understand what you are asking. it is extremely difficult to give advice while “in the blind”.
why don’t you set up a (very) small demo example representing what it is that you want to do with the pieces you are currently doing with colvars marked properly, so that we can propose proper alternatives. keep in mind, that other people cannot read your mind and do not instantly know what you are trying to achieve without a sufficient amount of context and background information. advice in general is easiest given, demonstrated and tested with small, tangible, and working examples.

axel.

7

Hi Julen, if you mean that you need half-harmonic or flat-bottom potentials, I don’t know if there is a specific LAMMPS code for that you could always add specially-designed bonds based on bond_style table.

Giacomo

8

Dear Axel and Giacomo,

Thanks for your replies and sorry for the late reply

We have been coding the types of restraints that we needed. In particular:

- we extended the “fix restrain bond” adding a new argument:

fix ID group-ID restrain keyword args

bond args = atom1 atom2 Kstart Kstop r0start r0stop
atom1,atom2 = IDs of 2 atoms in bond
Kstart,Kstop = restraint coefficients at start/end of run (energy units)

r0start,r0stop = equilibrium bond distance (distance units)

- we defined a new fix for lowerbound harmonic restraints. It is a harmonic restraint to keep atoms far apart from each other, and vanishes above the equilibrium bond distance.

lbound args = atom1 atom2 Kstart Kstop r0start r0stop

atom1,atom2 = IDs of 2 atoms in bond
Kstart,Kstop = restraint coefficients at start/end of run (energy units)

r0start,r0stop = equilibrium bond distance (distance units)

These harmonics are useful to implement restraint-based modelling approaches for the structural determination of biological molecules using molecular dynamics. They were first implemented using Colvars (see https://www.biorxiv.org/content/10.1101/642009v1), but we recently discovered that using the “fix restraint” option they can be used much more efficiently.

We have implemented them with a modified “fix restrain” class of the Lammps code:

https://github.com/lammps/lammps/blob/fa7085be077d9271f8400c27c63613cf1a652144/src/fix_restrain.cpp#L222

So it looks like:

/* ----------------------------------------------------------------------
apply harmonic lower-bound bond restraints
---------------------------------------------------------------------- */

void FixRestrain::restrain_lbound(int m)
{
int i1,i2;
double delx,dely,delz,fbond;
double rsq,r,dr,rk;

double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;

double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
double k = kstart[m] + delta * (kstop[m] - kstart[m]);
// New
double deq = deqstart[m] + delta * (deqstop[m] - deqstart[m]);

i1 = atom->map(ids[m][0]);
i2 = atom->map(ids[m][1]);

// newton_bond on: only processor owning i2 computes restraint
// newton_bond off: only processors owning either of i1,i2 computes restraint

if (newton_bond) {
if (i2 == -1 || i2 >= nlocal) return;
if (i1 == -1) {
char str[128];
sprintf(str,
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT,
ids[m][0],ids[m][1],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
} else {
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
if (i1 == -1 || i2 == -1) {
char str[128];
sprintf(str,
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT,
ids[m][0],ids[m][1],
comm->me,update->ntimestep);
error->one(FLERR,str);
}
}

delx = x[i1][0] - x[i2][0];
dely = x[i1][1] - x[i2][1];
delz = x[i1][2] - x[i2][2];
domain->minimum_image(delx,dely,delz);

rsq = delxdelx + delydely + delz*delz;
r = sqrt(rsq);
// Old
//dr = r - target[m];
// New
dr = r - deq;
rk = k * dr;

// force & energy
// Old
//if (r > 0.0) fbond = -2.0*rk/r;
//else fbond = 0.0;

//ebond += rk*dr;
//energy += rk*dr;
// New
if (dr < 0) {
if (r > 0.0) fbond = -2.0*rk/r;
else fbond = 0.0;

elbound += rk*dr;
energy += rk*dr;
} else {
fbond = 0.0;

elbound += 0.0;
energy += 0.0;
}

// apply force to each of 2 atoms

if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fbond;
f[i1][1] += dely*fbond;
f[i1][2] += delz*fbond;
}

if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx*fbond;
f[i2][1] -= dely*fbond;
f[i2][2] -= delz*fbond;
}
}

Among other smalls parts of fix_restrain.cpp and fix_restrain.h (we can send you the added code if needed). Then we can use the restraints in lammps together with harmonic restraints:

fix bond1 all restrain bond 1 50 0.0 1.0 50.0 1.0
fix bond2 all restrain bond 1 25 0.0 1.0 50.0 1.0
fix lbound1 all restrain lbound 1 5 0.0 1000.0 5.0 5.0

This code has been tested and works very well. We are wondering if it is possible to share this code with other Lammps users.

Could we include these changes in Lammps? Or it is better to include them as a USER plugin?

Best regards,
Julen

9

>We are wondering if it is possible to share this code with other Lammps users.

>Could we include these changes in Lammps?

Hi - You can submit your new fix to the LAMMPS GitHub site as a pull request,
to be considered for inclusion in the LAMMPS distro.

Thanks,
Steve