[lammps-users] fix shake for 4 atom molecule

LAMMPS LAMMPS Mailing List Mirror
1

Hi Lammps users!

I’d like to run a simulation using a force field that includes a rigid carbonate molecule. Lammps currently doesn’t include the ability to fix both bonds and angles for a 4-atom cluster, only three (e.g. water). So the first option is to use fix rigid, but I’d like to run at an npt or nvt ensemble and that option currently only works with nve. When I’ve tried running a mixed ensemble (CO3 under fix rigid with nve and the rest of the atoms with npt) I’m getting some very large fluctuations in temperature (spikes of several hundred degrees for a 300 degree simulation). I don’t know if this is due to the force field or from running a mixed ensemble though. Has anyone else tried running under a mixed ensemble before?

Another option I thought of is to modify fix_shake.cpp to allow for 4-atom clusters with fixed angles. That way I could run all in npt or nvt. Currently under find_clusters, for a fixed angle on a water molecule it does this:

if (nshake[i] == 2){
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
if (n < 0) continue;
if (angle_type[i][n] < 0) continue;
if (angle_flag[angle_type[i][n]]) {
shake_flag[i] = 1;
shake_type[i][2] = angle_type[i][n];
}
}

I could write something like this for a three atom cluster, but I’m a little concerned about what this will do to the degrees of freedom for the temperature calculation. In dof it contains this code:

for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
if (shake_flag[i] == 0) continue;
if (shake_atom[i][0] != tag[i]) continue;
if (shake_flag[i] == 1) n += 3;
else if (shake_flag[i] == 2) n += 1;
else if (shake_flag[i] == 3) n += 2;
else if (shake_flag[i] == 4) n += 3;
}

which looks correct even for a fixed four atom cluster. Is that right?

Thanks a lot for your help!

Cheers,
Andrew

2

You'd have to add a lot more to fix shake to do 4-body interactions,
like formulate the matrix equation and solve it.

In a couple weeks we'll release a fix rigid/nvt and fix rigid/npt
that will hopefully do what you want. We're just typing up
some loose ends with the new code.

Steve

3

Thanks for the response. I got into the code yesterday and realized the reason why lammps has this restriction. Perhaps I could create a larger square matrix in order to solve for the determinant but it’ll be easier to just wait for the new release. :slight_smile:

Cheers,
Andrew