Fortran API / lammps_extract_global

LAMMPS LAMMPS Development
1

I have been working to add features to /fortran/lammps.f90 the last couple of days, and I have implementations all the way through lammps_extract_global (working through library.h in order) working and tested.

I thought I’d run this by Axel, at the very least, before I proceed much further. In the /examples/COUPLE/fortran2 API, I made lammps_extract_global and similar functions into subroutines whose first argument is the variable being returned. This isn’t identical to the C interface, but it’s close, and saves the trouble of having an array of extra constants around.

I came up with an alternative for /fortran/lammps.f90, which is to resolve what kind of data lammps_extract_global (from the C API) is returning similar to the way Python does it, namely to have another parameter called “dtype” whose sole duty is to state to the compiler that the function is returning an integer, a 64-bit integer, a 64-bit real number, a string, etc. For example,

   integer (C_int) :: nlocal
   real (C_double) :: dt
   ! more code
   nlocal = lmp%extract_global('nlocal', LMP_INT)
   dt = lmp%extract_global('dt', LMP_DOUBLE)

where LMP_INT and LMP_DOUBLE are (public) constants defined within /fortran/lammps.f90.

First, would doing it this way be confusing to users and/or a maintenance issue (or otherwise icky)?

This led me to thinking along the lines of how to deal with LAMMPS_BIGINT and the LAMMPS_SMALLBIG/LAMMPS_BIGBIG/etc. settings. Since those settings are all set by preprocessor directives, Fortran won’t see them at all, so there isn’t an easy solution (such as defining LMP_BIGINT=LMP_INT if bigints are ints and LMP_BIGINT=LMP_INT64 if bigints are int64_ts). We can pull the size of bigint, tagint, and imageint out of extract_setting once LAMMPS is running, but that doesn’t let the user declare something like

   integer (C_bigint) :: ntimestep

at compile time (which is what they would do with the C or C++ APIs). Right now, I would lean toward having the Fortran API be agnostic to such things and rely on the user to “know” which integer kind to pass, but it would be nice to be “clever” about it (similar to the Python interface’s “dtype=LAMMPS_AUTODETECT” and subsequent type resolution, which of course is not possible in strongly-typed languages). Any thoughts?

2

If you need a polymorph return type in C where you decide the type at runtime, you can do it with:

typedef struct {
        int type;
        union {
           void *p;
           int i;
           long l;
           float f;
           double d;
       } data;
} polytype_t;

And then have constants for which you can set the member type to and that indicates which of the union items has been used e.g. data.i or data.d. Now Fortran doesn’t have a union, but you could use a derived type instead (i.e. the same as another struct) and use a little more memory but then have a self describing return value. This can be extended to arrays and multi-dimensional arrays and then you could also store the array dimensions as members etc.
This would work also with the various integer types, since you can look their sizes up at runtime and the Fortran interface would then use INTEGER(kind=4) and INTEGER(kind=8).

3

That’s an idea: create a derived type with overloaded assignment, equality, etc. operators that will assign to whatever is appropriate based on (private) member identities.

I like that for the purposes of internal implementation—it would still be on the user to know what to put on the left-hand-side of the assignment, but it might eliminate the proliferation of redundantly coded generic functions. It also makes it possible to solve the “LAMMPS_BIGBIG” problem, at least in principle.

Thanks.

4

I messed around with this for a while, and it turns out the derived type plus overloaded assignment operator idea works with minimal code duplication if one is OK with making copies of the data from LAMMPS. If we want to return pointers to LAMMPS data with extract_global and friends, then it’s harder.

For example,

ntimestep = lmp%extract_global("ntimestep")

should in principle allow “ntimestep” to update as LAMMPS runs. This is possible, but it’s complicated by the fact that the “=>” operator is not overload-able combined with Fortran’s rules about polymorphic dummy arguments with the pointer attribute. It’s possible to do what I want to do, but not without code duplication, or at least I haven’t thought of a way to do it cleanly yet. I’ll keep thinking about this issue…

5

You don’t want to change these global properties directly. That can create all kinds of trouble, since in most cases there would need to be some additional changes of derived or related data and sometimes even a full re-initialization.

I think in python this is the same.

6

I completely agree that users should not change values returned from extract_global, but they might want to read them again later (e.g., the time step will update as the run continues). We should definitely warn the user not to modify the results from extract_global, but there is nothing physically preventing them from doing it through the C API (and there is a big, orange-colored warning message in the documentation telling them not to modify the data, which I will replicate in the Fortran documentation as well). I did notice that Python made copies, but I assumed that was at least in part because of the incompatibilities between Python and C/C++ (i.e., Python doesn’t have pointers, or at least cannot map C arrays directly onto Python lists).

My original implementation in examples/COUPLE/fortran2 (fall 2012 or earlier) made copies of the data returned from LAMMPS, without the possibility of sending anything back, but almost immediately we encountered a use case (for extract_atom) that required direct access to the data without copies being made. See this post from almost 10 years ago (which makes me feel really old…): questions about library.cpp - #2 by sjplimp.

I actually managed to get something working, though it required more subroutines than I was hoping for (but probably fewer lines of code in the end). Now my biggest question is about error messages: the mechanism most of LAMMPS uses requires the preprocessor, and if I recall correctly you wanted to avoid doing that. This rules out using the __LINE__, __FILE__ mechanism used in the Error class. For now, I’m just dumping a message to standard error (ERROR_UNIT in Fortran) and issuing STOP if the user passes an incompatible pointer to extract_global, but I don’t expect that to be a long-term solution.

7

I figure this thread is as good a place as any to continue this conversation with @akohlmey.

As I’m certain you already saw, I pushed Fortran implementations of all functions in library.h up through lammps_extract_global (i.e., get_thermo, extract_box, reset_box, memory_usage, get_mpi_comm, extract_setting, and extract_global). I also wrote unit tests for those functions, which is what prompted my bug report about the Fortran tests not building earlier. I’m not in any way trying to push anyone look at them quickly, just making sure you know what’s in there. If that implementation, particularly the overloaded assignment operator I implemented to make extract_global work, is not desirable in some way, please let me know and I’ll try to come up with an alternative. I do not plan to push further changes to that particular pull request unless someone asks me to; any further changes I suggest for the Fortran API will go through another branch.

I want to ask for preferences on what to do with the “type” and “style” parameters passed to extract_compute (which will also affect extract_variable, extract_fix, etc.) The C interface defines them as global constants, such as LMP_STYLE_GLOBAL, LMP_STYLE_ATOM, and LMP_STYLE_LOCAL. The Python interface (almost accidentally) defines them as lammps.LMP_STYLE_GLOBAL, lammps.LMP_STYLE_ATOM, and so forth, though using “from lammps import *” would put them in the global namespace, as with the C API.

For the Fortran API, I can see three ways to do this: (1) as public parameters in the module (like the C library), (2) as type-bound “constants” like the Python module, and (3) as a “nested” derived type. The difference in the calling signature would be something like

   ! option (1)
   compute = lmp%extract_compute('1', LMP_STYLE_LOCAL, LMP_TYPE_SCALAR)

versus

   ! option (2)
   compute = lmp%extract_compute('1', lmp%STYLE_LOCAL, lmp%TYPE_SCALAR)

versus

   ! option (3)
   compute = lmp%extract_compute('1', lmp%style%local, lmp%type%scalar)

Style (2) could also use lmp%LMP_STYLE_LOCAL and such if desired.

Thoughts? I lean toward option (3) as the most “elegant” solution, mainly because it only requires the user to use type(lammps) instead of that plus several parameters. Option (2) would function similarly. Option (1), however, is closest to what C and Python do currently.

8

As it so happens, I’ve spend the last few hours, reviewing and testing most of those changes and applying some updates and (minor) corrections as well as a “property” example similar to the one for the C library.

Absolutely, I am about to finish the existing pull request as it looks pretty good already and works much better than I had expected. You have done a great job there.

I have been thinking about error handling and I think we just need to add a function to the C-library interface that can trigger the error calls.

Since this is new, we are free to choose. The main requirement is that the words and their order remain the same and that is fulfilled by all 3 options. While we try to keep all three interfaces similar, there are also differences when the expected way of how something works in a given language is different.

9

The only thing there is how to handle line numbers without resorting to the preprocessor.

Honestly, it worked better than I expected, too. I wanted to overload the => operator rather than the = operator, but that evidently is not allowed. Fortunately, this worked. (It was largely due to your suggestion that I thought of doing it that way, so…thanks!) It also gives compile-time errors if you try to assign a non-pointer object to a pointer from extract_global, which I’ve already done by accident twice when testing. Explaining that error message to users might be a tweak I’ll apply to the next PR.

I’ll proceed with (3) for now, but it shouldn’t be immensely difficult to change to (1) or (2) later if we choose to do that.

Thank you for your guidance and help!

10

We have to explicitly provide data. But I am thinking for simplicity to make it line number -1 and file name ‘(unknown)’ and set this already in the C-library interface. So that the C-library signature would be:
void lammps_error(void *handle, int error_code, const char *error_text);

for the error code there are (additive) enumerated constants:
LAMMPS_ERROR_WARNING=0, LAMMPS_ERROR_ONE=1, LAMMPS_ERROR_ALL=2, LAMMPS_ERROR_WORLD=4, LAMMPS_ERROR_UNIVERSE=8
so that it can branch out to the different Error class calls. Unknown error codes would just trigger a warning about them.

11

The file name is known, as everything is in lammps.f90. Setting the line number to some (negative) constant and telling the Error class that this constant means “error from Fortran interface” might be a good solution.

12

We could also do some trickery and have either sed or cpp generate lammps.f90 from another file (say, lammps.f90.in) when LAMMPS is built, replacing __LINE__ with the line number, sort of like how doc/utils/sphinx-config/conf.py is constructed.

13

I would like this to be a generic function and not just a utility to be called from the Fortran interface.
I have just pushed my changes to your PR branch. Please check it out and add your comments/suggestions (or implement them :wink: ).

14

I probably won’t do it soon, but when it comes time to implement the remaining 30 functions, what is the status of lammps_gather and lammps_scatter? They seem to do the same thing as lammps_gather_atoms and lammps_scatter_atoms except for a couple of modifications (e.g., to gather/scatter per-atom fix or compute data). Will they eventually replace lammps_gather_atoms and lammps_scatter_atoms? Do we need both in the Fortran API?

Thank you for your help in fixing the bugs revealed during the automated checks.

15

These functions have not been reviewed and properly tested and documented. There were some changes incorporated to accommodate some external package, but since then there has not been sufficient free time to complete the review, documentation, testing, and refactoring of the library interface.

It is an ongoing work in progress, but the incentive to complete is not very large as we are moving from frequently used to more obscure functions.

At the moment the various ongoing internal refactoring projects in LAMMPS have priority and we have far too few people interested in collaborating on LAMMPS maintenance. You just need to look at who does what to see how urgent we need motivated and competent help.

16

@akohlmey: I was thinking about the frequent logical-to-int or int-to-logical conversions that occur in the Fortran interface, and it occurred to me that we could eliminate some redundant programming (e.g., all of the “lmp_config_has_XXX” routines) by having the C functions return a _Bool value (or bool if using <stdbool.h>) instead of an int. Fortran can call such C functions using the LOGICAL(c_Bool) type and kind (c_Bool is part of ISO_C_BINDING).

For example, lammps_config_package_count is called directly by the Fortran interface, with no “wrapper”: it returns an integer and doesn’t need a pointer to LAMMPS, so it “just works” directly. The rest of those “_config” functions would work similarly if the “int-to-logical” conversion became a “_Bool-to-logical(c_Bool)” non-conversion.

Caveats with this idea: (1) it would mean changing the return values of certain things in library.h and library.cpp to be _Bool instead of int; (2) users would get a warning if they pass an ordinary logical (not of kind(c_bool)) to it, though passing the constants .TRUE. and .FALSE. should be fine, and there should also be no warning for assigning return values from logical(c_Bool) to logical of default kind or using logical(c_Bool) in if or while statements; (3) there might be ramifications for the Python API, too.

What is your take on this sort of conversion? Is using _Bool at all worth thinking about? For that matter, there are other places in LAMMPS that use int when bool could be used; are there plans to convert those, or is that something best left alone?

17

No. We want to keep the header requirements on library.h as minimal as possible.
Even mpi.h is not included (and the corresponding function exposed) unless an explicit define is set.

Not at the moment. We don’t object to people using bool, but for anything that requires MPI communication, it must be an int, since we stick with the C interface to MPI.

18

_Bool itself does not require any headers (only bool does, if using it in C). I get your point, though.

19

I don’t know that this is true: in the MPI 4 standard (and probably the MPI 3 standard, too, but I don’t have that in front of me), there is a table and accompanying statement that says, “If there is an accompanying C++ compiler, then the datatypes in table 3.4 are also supported in C and Fortran”; Table 3.4 contains MPI_CXX_BOOL (for bool) as well as three more supported C++'s std::complex<kind> types, where type is float, double, or long double. There is also MPI_C_BOOL for _Bool in C.

That said, I won’t expend any time on the _Bool idea in the library. Thanks.

20

@akohlmey The C library interface currently has several functions that accept a string and its length; an example is lammps_get_os_info, which has the prototype

void lammps_get_os_info(char *buffer, int buf_size)

However, malloc and friends require the length to be type size_t. This causes some issues on the Fortran side, where integer kinds in function calls are checked more thoroughly by the compiler. I see two possible solutions: (1) Update the C interface to have size_t buf_size and similar edits elsewhere, or (2) use INT(buf_size, KIND=c_size_t) as the argument to lammps_malloc or INT(buf_size, KIND=c_int) as arguments to C library functions on the Fortran side.

Neither should cause significant disruption to the user interface, as the pass-by-value nature of C should promote int to size_t when necessary in function calls.

Which do think is preferable?