Error: fix widom region extends outside simulation box

LAMMPS LAMMPS General Discussion
1

Hi folks-

I am currently running the 2 August 2023 Development version of LAMMPS built from source - this question pertains specifically to the MC package.

I have a trajectory output from an atomistic MD simulation of a polymer system, and am trying to calculate Henry’s Law coefficients using the fix widom command. I’ve written a code using the LAMMPS Python API that creates a “free volume” map of my simulation cell using probe atoms, essentially testing whether there are overlaps nearby, to try and reduce the high deltaE insertions. Essentially, I am creating this “map” for a frame, and want to create a small cubic region around each probe atom that is in the free space, then take the union of all those regions (on the order of 10^3 - 10^4 of them), and use that for the region flag input to fix widom.

My issue is coming from a much less sophisticated test case. I have a molecule template built correctly (though this is reproducible with particle insertions), and am passing a region in to fix widom which is the dimensions of my simulation box, but slightly smaller (i.e. I am increasing the low boundaries and decreasing the high boundaries by a length of 1.2 Angstroms). When I try to run this, I get the following error (this is after reading in my data file to get atom types, topology, potentials, etc.):

read_dump configs.lammpstrj 0 x y z wrapped no format native
Scanning dump file …
Reading snapshot from dump file …
orthogonal box = (0.71337975 0.71337975 0.71337975) to (34.464498 34.464498 34.464498)
2825 atoms before read
2825 atoms in snapshot
0 atoms purged
2825 atoms replaced
0 atoms trimmed
0 atoms added
2825 atoms after read
region thebox block 1.9133797515978204 33.26449800840162 1.9133797515978204 33.26449800840162 1.9133797515978204 33.26449800840162 side in units box
fix mywidom all widom 1 1000 0 2917 308.15 mol ch4 full_energy region thebox
ERROR: Fix widom region thebox extends outside simulation box (src/MC/fix_widom.cpp:118)

Now, I can see from the echoed orthogonal box output after reading the snapshot that the region bounds I provided are, in fact, NOT outside the simulation box. My buffer is also sufficient such that, if the center of mass of my molecule were placed exactly on the edge of the region, no part of it would extend outside the simulation box.

Two things that are worth noting: this is an inconsistent error, in that I will occasionally run this without changing anything and it works without error. To me that doesn’t make sense.

Additionally, I noticed that the line of the error in fix_widom.cpp which is echoed, line 118, seems to refer specifically to the case of a triclinic simulation box, which I don’t have here.

Any input or similar experiences to this case would be helpful. I saw this post from earlier this year which had a similar issue, but didn’t really come to a resolution.

Thanks,

Sam

2

The problem is likely due to the information that you are not showing us. For example, the dimensions of an initial box may not be accurate by the time you start the run with fix widom.

I disagree with your assessment. That actually had multiple viable suggestions and resolutions.

3

For example, the dimensions of an initial box may not be accurate by the time you start the run with fix widom.

The issue here that I’m stuck on is that what’s failing is not a run command, but the actual initialization of fix widom. Meaning even though I can see in the line-by-line output above that the region’s bounds are well within the box dimensions but there is an error being thrown from this part of FixWidom() in fix_widom.cpp:

// error checks on region and its extent being inside simulation box
region_xlo = region_xhi = region_ylo = region_yhi = region_zlo = region_zhi = 0.0;
if (region) {
if (region->bboxflag == 0)
error->all(FLERR,“Fix widom region {} does not support a bounding box”, region->id);
if (region->dynamic_check())
error->all(FLERR,“Fix widom region {} cannot be dynamic”, region->id);
region_xlo = region->extent_xlo;
region_xhi = region->extent_xhi;
region_ylo = region->extent_ylo;
region_yhi = region->extent_yhi;
region_zlo = region->extent_zlo;
region_zhi = region->extent_zhi;
if (triclinic) {
if ((region_xlo < domain->boxlo_bound[0]) || (region_xhi > domain->boxhi_bound[0]) ||
(region_ylo < domain->boxlo_bound[1]) || (region_yhi > domain->boxhi_bound[1]) ||
(region_zlo < domain->boxlo_bound[2]) || (region_zhi > domain->boxhi_bound[2]))
error->all(FLERR,“Fix widom region {} extends outside simulation box”, region->id);
} else {
if ((region_xlo < domain->boxlo[0]) || (region_xhi > domain->boxhi[0]) ||
(region_ylo < domain->boxlo[1]) || (region_yhi > domain->boxhi[1]) ||
(region_zlo < domain->boxlo[2]) || (region_zhi > domain->boxhi[2]))
error->all(FLERR,“Fix widom region {} extends outside simulation box”, region->id);
}

The line numbering given in the error tells me that this comes from the if (triclinic) part of the loop, which is strange because I’ve defined an orthogonal simulation box as was shown in the output from the command on the previous line.

I’ll try to fiddle with this a bit more to see if I can figure out where the issue is. Thank you for the quick reply, Axel!