QUIP, OpenKIM, and boundary conditions

OpenKIM
1

I started working on making QUIP potentials into OpenKIM models. Most of the code is straightforward, except for the boundary conditions. Perhaps this is better handled by skype (maybe next week), but here's a description of what I've come up with so far.

Right now QUIP has a connection object and an atoms object which are used by the potential calculation routines. The calculation routines loop over all atoms (most of them, including most relevantly GAP, understand ghost atoms), loops over the number of neighbors of each atom, and asks for the index and relative position to each neighbor. The connection object doesn't actually store the relative position. It has an integer lattice vector shift, and the neighbor-getting routines calculates the relative position from the absolute positions, the lattice vectors, and the integer shift.

I initially thought it'd be easiest to create a QUIP connection object from the OpenKIM neighbor list on the fly, and then call the usual QUIP calculation routines. I'm thinking, however, that this isn't possible in general, because there's no way to get the lattice vectors in KIM.

QUIP also assumes a full neighbor list.

Does it looks like the following statements are correct?

1. To support PURE I can create a connection object where the integer lattice shift is always 0. For half lists I can just complete to a full list, and for a full list I can just ignore i > j, then pretend it's a half list and do the same thing.
2. To support RVEC I don't think I can do this. I'd have to put a pointer to the KIM object in the QUIP connection, and directly return relative positions. This would preclude ITERATOR mode, because QUIP assumes LOCATOR.
3. To support MI_PBC I can construct a lattice and create a connection object like in 1
4. To support CLUSTER I can just ignore the lattice and create a connection object

How much of a restriction (in practice, given how tests are actually coded) would it be if RVEC wasn't supported (which is ironic, since RVEC is most like what QUIP's actual calculation routines want), or if RVEC was supported but only with LOCATOR?

                                                    Noam

2

And, as a followup, suppose I think I’ve done the PBC stuff correctly. Is there any pre-existing code that will allow me to test it, e.g. a series of tests that pass the same configuration to the model with differently encoded neighbor lists?

3

Indeed, thanks to Ryan’s help in the other thread my QUIP MODEL_DRIVER is compiling successfully. I can start with the various example tests that run the test with various periodic boundary conditions. Has there been any thought of developing a general PBC test of the models, so that if a model claim to support multiple NBC methods it really gives correct, or at least consistent, answers?

Noam

4

Indeed, thanks to Ryan's help in the other thread my QUIP MODEL_DRIVER is
compiling successfully. I can start with the various example tests that run
the test with various periodic boundary conditions. Has there been any
thought of developing a general PBC test of the models, so that if a model
claim to support multiple NBC methods it really gives correct, or at least
consistent, answers?

The numerical derivatice check vc_forces_numer_deriv performs a numerical derviative check on a model for all the following NBC modes:

         NEIGH_RVEC_H
         NEIGH_PURE_H
         NEIGH_RVEC_F
         NEIGH_PURE_F
         MI_OPBC_H
         MI_OPBC_F
         CLUSTER

This would be a good way to check whether all are supported correctly. There is more information in the README file in TESTS/vc_forces_numer_deriv

-Ellad

5

Hi Noam,

Good. Glad it is working.

One thing you can do is use the vc_forces_numer_deriv Test. This runs consistency checks on the computed forces by doing a numerical derivative of the energy. The test checks EVERY NBC method that the Model claims to support. It doesn't do an automatic comparison between NBCs but the data is there for you to compare manually.

Ryan

6

Ellad beat me to it!

Of course, we are also always open to contributions from the community of new ways to better verify the Models...

Cheers,

Ryan

7

Great, thanks - I’ll try that.

Noam

8

I started working on making QUIP potentials into OpenKIM models. Most of the code is straightforward, except for the boundary conditions. Perhaps this is better handled by skype (maybe next week), but here's a description of what I've come up with so far.

Right now QUIP has a connection object and an atoms object which are used by the potential calculation routines. The calculation routines loop over all atoms (most of them, including most relevantly GAP, understand ghost atoms), loops over the number of neighbors of each atom, and asks for the index and relative position to each neighbor. The connection object doesn't actually store the relative position. It has an integer lattice vector shift, and the neighbor-getting routines calculates the relative position from the absolute positions, the lattice vectors, and the integer shift.

I initially thought it'd be easiest to create a QUIP connection object from the OpenKIM neighbor list on the fly, and then call the usual QUIP calculation routines. I'm thinking, however, that this isn't possible in general, because there's no way to get the lattice vectors in KIM.

QUIP also assumes a full neighbor list.

Does it looks like the following statements are correct?

1. To support PURE I can create a connection object where the integer lattice shift is always 0. For half lists I can just complete to a full list

Yes.

, and for a full list I can just ignore i > j, then pretend it's a half list and do the same thing.

No. Full lists are not required to be "symmetric". So, if j is a neighbor of i, it is not necessarily true that i is a neighbor of j. In particular, ghost atoms have no neighbors. In the full list case ghost atoms can have any index (they are not required to be listed after the "real atoms").

2. To support RVEC I don't think I can do this. I'd have to put a pointer to the KIM object in the QUIP connection, and directly return relative positions. This would preclude ITERATOR mode, because QUIP assumes LOCATOR.

I don't follow. The major differences between PURE and RVEC are that (i) the rij vectors are returned and (ii) "image" atoms are supported by RVEC.

Image lead to the situation where particle i can be a neighbor of itself. For example, the neighbor list for particle 15 may include particle 15. However, the rij vector for this neighbor must be non-zero, indicating that there is an "image" of particle 15 located relative to particle 15 by the vector rij.

(i) seems to pose no real problem; if you can support PURE, then you can return the rij too.

(ii) might prove to be a difficulty, but I don't know enough about your connection object to be able to tell.

3. To support MI_PBC I can construct a lattice and create a connection object like in 1

Yes.

4. To support CLUSTER I can just ignore the lattice and create a connection object

Ok, seems right.

How much of a restriction (in practice, given how tests are actually coded) would it be if RVEC wasn't supported (which is ironic, since RVEC is most like what QUIP's actual calculation routines want), or if RVEC was supported but only with LOCATOR?

RVEC is pretty important. Currently about half of all Tests on openkim.org support ONLY RVEC.

Most of the current Tests on openkim.org support both LOCATOR and ITERATOR mode. So, RVEC with only LOCATOR would be significantly better than no RVEC at all.

Cheers,

Ryan

9

, and for a full list I can just ignore i > j, then pretend it’s a
half list and do the same thing.

No. Full lists are not required to be “symmetric”. So, if j is a
neighbor of i, it is not necessarily true that i is a neighbor of j. In
particular, ghost atoms have no neighbors. In the full list case ghost
atoms can have any index (they are not required to be listed after the
“real atoms”).

Hmm. That’s going to require some thought. Other than doing unnecessary computation, is pretending that full lists are symmetric and that ghost atoms are the same as real atoms a problem?

  1. To support RVEC I don’t think I can do this. I’d have to put a
    pointer to the KIM object in the QUIP connection, and directly return
    relative positions. This would preclude ITERATOR mode, because QUIP
    assumes LOCATOR.

I don’t follow. The major differences between PURE and RVEC are that (i)
the rij vectors are returned and (ii) “image” atoms are supported by RVEC.

Image lead to the situation where particle i can be a neighbor of itself.
For example, the neighbor list for particle 15 may include particle 15.
However, the rij vector for this neighbor must be non-zero, indicating
that there is an “image” of particle 15 located relative to particle 15 by
the vector rij.

(i) seems to pose no real problem; if you can support PURE, then you can
return the rij too.

(ii) might prove to be a difficulty, but I don’t know enough about your
connection object to be able to tell.

Image atoms are fine, QUIP does those all the time. But it doesn’t actually store Rij, just |Rij| and integer 3-vectors (shift) to be multiplied by the lattice vectors such that
Rij = r_j + lattice.shift - r_i
And since I don’t have access to lattice in the OpenKIM model-driver, I can’t create the QUIP connection object from a RVEC OpenKIM neighbors - I don’t know what to set shift to. I’ll have to modify the QUIP object. It’s not an insurmountable problem, but it’s just a bit more coding.

RVEC is pretty important. Currently about half of all Tests on
openkim.org support ONLY RVEC.

Most of the current Tests on openkim.org support both LOCATOR and ITERATOR
mode. So, RVEC with only LOCATOR would be significantly better than no
RVEC at all.

OK, I guess that’s motivation to actually do it, RVEC that is. I’m not sure about ITERATOR yet. Maybe I can do it at the same time (by converting to a locator, basically, as I create the QUIP connection object).

Thanks for the clarifications.

Noam

10

, and for a full list I can just ignore i > j, then pretend it's a
half list and do the same thing.

No. Full lists are not required to be "symmetric". So, if j is a
neighbor of i, it is not necessarily true that i is a neighbor of j. In
particular, ghost atoms have no neighbors. In the full list case ghost
atoms can have any index (they are not required to be listed after the
"real atoms").

Hmm. That's going to require some thought. Other than doing unnecessary
computation, is pretending that full lists are symmetric and that ghost
atoms are the same as real atoms a problem?

I would think, yes. The ghost atoms do not contribute their energy to the total energy of the configuration. Also, if a ghost atom is given neighbors it will generate forces contributions on those neighbors and this would be incorrect according to the way the KIM API defines the force contribution.

2. To support RVEC I don't think I can do this. I'd have to put a
pointer to the KIM object in the QUIP connection, and directly return
relative positions. This would preclude ITERATOR mode, because QUIP
assumes LOCATOR.

I don't follow. The major differences between PURE and RVEC are that (i)
the rij vectors are returned and (ii) "image" atoms are supported by RVEC.

Image lead to the situation where particle i can be a neighbor of itself.
For example, the neighbor list for particle 15 may include particle 15.
However, the rij vector for this neighbor must be non-zero, indicating
that there is an "image" of particle 15 located relative to particle 15 by
the vector rij.

(i) seems to pose no real problem; if you can support PURE, then you can
return the rij too.

(ii) might prove to be a difficulty, but I don't know enough about your
connection object to be able to tell.

Image atoms are fine, QUIP does those all the time. But it doesn't
actually store Rij, just |Rij| and integer 3-vectors (shift) to be
multiplied by the lattice vectors such that
  Rij = r_j + lattice.shift - r_i
And since I don't have access to lattice in the OpenKIM model-driver, I
can't create the QUIP connection object from a RVEC OpenKIM neighbors - I
don't know what to set shift to. I'll have to modify the QUIP object.
It's not an insurmountable problem, but it's just a bit more coding.

Can't you follow the same procedure that you would for the PURE case but take the rij vector from the KIM neighbor list and pass it along to the model? Ah, I see, it just doesn't work with the currenct framework of your connection object.... OK, yeah. More coding...

Cheers,

Ryan

11

So I now have a model_driver_QUIP that gives the same forces, basically consistent numerical and analytical, for every NBC method in vc_forces_numer_deriv, but I have 2 related problems/questions:

  1. Most of the forces fail the quantitative Ridder’s method tests by a factor of 2 or so. Is that probably just an issue with the QUIP implementation itself?

  2. vc_forces_numer_deriv doesn’t actually generate any periodic images, so it’s not fully testing the different NBCs. Is there any test that does?

Noam

12

This would be a good way to check whether all are supported correctly.
There is more information in the README file in
TESTS/vc_forces_numer_deriv

So I now have a model_driver_QUIP that gives the same forces, basically
consistent numerical and analytical, for every NBC method in
vc_forces_numer_deriv, but I have 2 related problems/questions:

Great.

1. Most of the forces fail the quantitative Ridder's method tests by a
factor of 2 or so. Is that probably just an issue with the QUIP
implementation itself?

I'm not exactly sure what you mean by this. If you are referring to the "pass/fail" that the verification check prints at the end of each line, then it is a known issue that in most cases even correct Models have nearly every line listed as "fail". Ellad Tadmor implemented this routine based on the discussion in Numerical Recipes. However, it seems that the criterion used to determine the pass/fail state is too strict in practice.

2. vc_forces_numer_deriv doesn't actually generate any periodic images, so
it's not fully testing the different NBCs. Is there any test that does?

That is a good point. We don't have a better "all around stress-test" type of check. However, in our examples are some individual Tests that use these features:

ex_test_Ar_FCCcohesive_NEIGH_PURE: computes the cohesive energy of a
   periodic crystal using ghost-atoms

ex_test_Ar_FCCcohesive_NEIGH_RVEC: computes the cohesive energy of a
   periodic crystal using a one-atom configuration and many
   image-atoms

These are not particularly general tests of these features, but they are something.

Cheers,

Ryan

13

I guess I’d relax the criterion for failure, so as to distinguish typical rounding errors (whether or not Numerical Recipes thinks they are small enough) from clear coding errors. I’d suggest taking the criterion and multiplying by 100 or so?

Jim

14

So I now have a model_driver_QUIP that gives the same forces, basically
consistent numerical and analytical, for every NBC method in
vc_forces_numer_deriv, but I have 2 related problems/questions:

Great.

I’ve also confirmed that it’s consistent with the same structure evaluated directly by QUIP, so not only is it internally consistent, it’s also validated against the original implementation.

I’m not exactly sure what you mean by this. If you are referring to the
“pass/fail” that the verification check prints at the end of each line,
then it is a known issue that in most cases even correct Models have
nearly every line listed as “fail”. Ellad Tadmor implemented this routine
based on the discussion in Numerical Recipes. However, it seems that the
criterion used to determine the pass/fail state is too strict in practice.

Yes, exactly. So I won’t worry about the “fail”.

That is a good point. We don’t have a better “all around stress-test”
type of check. However, in our examples are some individual Tests that
use these features:

ex_test_Ar_FCCcohesive_NEIGH_PURE: computes the cohesive energy of a
periodic crystal using ghost-atoms

ex_test_Ar_FCCcohesive_NEIGH_RVEC: computes the cohesive energy of a
periodic crystal using a one-atom configuration and many
image-atoms

I’ll try those. One thing that might be useful for others, which I’ve implemented in a particular and not necessarily general way is to print out (I used QUIP’s format, but perhaps you could do it in some standard format if you can find one) the input geometries and output quantities that are being evaluated by these tests. That’s very useful for validation, i.e. comparing to previous or independent implementations.

Noam

15

I agree that having this capability without so many false positives is useful, even if it requires some empiricism. I think that a much smaller factor would be sufficient in my case. I’ll post an update once I figure out what that factor is.

Noam

16

Agreed. We have been meaning to do something like this, but just haven't taken the time to carefully consider how we should choose a better criterion value.

Ryan

17

So I now have a model_driver_QUIP that gives the same forces, basically
consistent numerical and analytical, for every NBC method in
vc_forces_numer_deriv, but I have 2 related problems/questions:

Great.

I've also confirmed that it's consistent with the same structure evaluated
directly by QUIP, so not only is it internally consistent, it's also
validated against the original implementation.

Very good.

I'm not exactly sure what you mean by this. If you are referring to the
"pass/fail" that the verification check prints at the end of each line,
then it is a known issue that in most cases even correct Models have
nearly every line listed as "fail". Ellad Tadmor implemented this routine
based on the discussion in Numerical Recipes. However, it seems that the
criterion used to determine the pass/fail state is too strict in practice.

Yes, exactly. So I won't worry about the "fail".

OK

That is a good point. We don't have a better "all around stress-test"
type of check. However, in our examples are some individual Tests that
use these features:

ex_test_Ar_FCCcohesive_NEIGH_PURE: computes the cohesive energy of a
         periodic crystal using ghost-atoms

ex_test_Ar_FCCcohesive_NEIGH_RVEC: computes the cohesive energy of a
         periodic crystal using a one-atom configuration and many
         image-atoms

I'll try those. One thing that might be useful for others, which I've
implemented in a particular and not necessarily general way is to print out
(I used QUIP's format, but perhaps you could do it in some standard format
if you can find one) the input geometries and output quantities that are
being evaluated by these tests. That's very useful for validation, i.e.
comparing to previous or independent implementations.

Yeah, that is a good idea. I'll put that on my list.

Ryan

18

My test suggests that this might not be easy. The force difference in my implementation is always less than about 1e-7, and the ratio of the estimated error and force difference is almost always small (mostly < 10, the rest of the time < 100). However, every so often the numerical error estimate is tiny, 1e-14, and then the ratio between the force difference and estimated error becomes huge (> 10^7).

Noam