PBC for negative coordinates

[SSchott]

Hi Leandro,

I was testing the pbc implementation for packmol-memgen, and noted that the pbc function does not seem to play well with negative coordinates, or at least that is the impression I am getting from my tests, where all molecules end up in a small "positive corner". Did you play around with this?

[lmiq]

What exactly to you have tested? The PBCs in packmol are somewhat "special" in the sense that, appart the periodicity, they also define constraints for the coordinates. That is, if we set

tolerance 2.0
pbc 40.0 40.0 40.0

This will automatically add a constraint of the type inside box -1.0 -1.0 -1.0 41.0 41.0 41.0 to all atoms of the system. This is necessary for the other possible constraints of the system to "make sense", since they will not be periodic.

So, if for example we set explicitly a constraint like inside box -10. -10. -10 10. 10. 10, and we want PBCs of side 40, we should use pbc -20. -20. -20. 20. 20. 20. to avoid the intersection of the pbc with the box defined.

Maybe that is the issue, and it can be clarified in the documentation.

[SSchott]

I am using something like the following:

pbc -41.84 -41.84 -44.61 41.84 41.84 48.75

structure PROT0.pdb
  number 1
  fixed 0. 0. 0 0. 0. 0.
  radius 1.5
end structure

structure DOPE.pdb
  nloop 20
  number 75
  atoms 1 11
    below plane 0. 0. 1. -18.0
  end atoms
  atoms 79 126
    over plane 0. 0. 1. -4.0
  end atoms
end structure

structure WAT.pdb
  nloop 20
  number 10093
  outside box -41.84 -41.84 -23.0 41.84 41.84 23.0
end structure

I get all molecules initially in the positive corner, and they do not move towards the negative side. Additionally, the GENCAN function increases in value rather than decreasing.

I noticed the use of floor here, which might be the culprit (?):

packmol/src/setibox.f90

Line 18 in d236cb4

x = x - floor(x / pbc_length(1)) * pbc_length(1)

That should be fine for positive values, but for negative values, it would cause a wrong shift? I compiled it changing that for nint, and though it is slow, it goes "in the right direction" apparently (though still running).

[lmiq]

You're right, that is not correct, but the support for negative coordinates must be fixed in other places as well. I'll try to do it asap.

[SSchott]

Ok, thanks! Changing those 3 coordinates to nint is packing something sensical at least 👍

[SSchott]

Reaching to the WAT section also caused the function value to increase, so I let you to it.

BTW, noted the missing pbc_box(1) here:

packmol/src/getinp.f90

Line 187 in d236cb4

write(*,"(a, 6f8.2)") ' PBC Reference box: ', pbc_box(2), pbc_box(2), pbc_box(3), pbc_box(4), pbc_box(5), pbc_box(6)

[lmiq]

Indeed, I spotted that one now as well. Thanks.

I don´t think you'll be able to pack with negative coordinates now. In the mean time, I would suggest just moving everything to the positive side.

[SSchott]

Even after shifting everything to positive coordinates, the packing behaves weirdly.
Current input:

pbc 0. 0. 0. 83.68 83.68 93.36

structure PROT0.pdb
  number 1
  fixed 41.84 41.84 44.61 0. 0. 0.
  radius 1.5
end structure

structure DOPE.pdb
  nloop 20
  number 75
  atoms 1 11
    below plane 0. 0. 1. 26.61
  end atoms
  atoms 79 126
    over plane 0. 0. 1. 40.61
  end atoms
end structure

structure DOPG.pdb
  nloop 20
  number 25
  atoms 1 13
    below plane 0. 0. 1. 26.61
  end atoms
  atoms 81 128
    over plane 0. 0. 1. 40.61
  end atoms
end structure

structure DOPE.pdb
  nloop 20
  number 69
  atoms 1 11
    over plane 0. 0. 1. 62.61
  end atoms
  atoms 79 126
    below plane 0. 0. 1. 48.61
  end atoms
end structure

structure DOPG.pdb
  nloop 20
  number 23
  atoms 1 13
    over plane 0. 0. 1. 62.61
  end atoms
  atoms 81 128
    below plane 0. 0. 1. 48.61
  end atoms
end structure

structure WAT.pdb
  nloop 20
  number 10093
  outside box 0. 0. 21.61 83.68 83.68 67.61
end structure

structure K+.pdb
  nloop 20
  number 44
  outside box 0. 0. 21.61 83.68 83.68 67.61
end structure

Here GENCAN value vs loop nr:

And here how the "best packed" version looks for now:

[lmiq]

Can you test the version in this PR branch? (or alternatively, send me your input files)?

git clone https://github.com/m3g/packmol
cd packmol
git fetch origin pbc_negative_coordinates
git switch --track origin/pbc_negative_coordinates
make
./packmol < inputfile.inp

In parallel: In the meantime I verified that the outside box restriction might behave badly in these cases, if cutting the region where the molecules are put in half. Maybe you want to use inside box or over plane / below plane instead.

[SSchott]

I gave it a try, and it runs as expected now, with or without negative coordinates now. The outside box def indeed causes trouble though, so had to use above/below plane:
Black: negative coordinates Red: outside box Green: positive coordinates

[lmiq]

Great. I will release this fix and probably work on improving the outside box definition. Thanks for the feedback.