questions about library.cpp

Hi,

I am trying to add some functions to library.cpp. I am wondering why in library.cpp you guys always put two dimension pointers from lammps into one dimension pointers. Like when coordinates are got from lammps in the function below, x is a two dimension like x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],… But one dimension copy is used to store x by expand the two dimension pointers x. I think it should be easier to just use two dimension maybe **copy do the same job. I am not really familiar with c++/c. The command line “ lmp->memory->create(copy,countnatoms,“lib/gather:copy”);” is creating vector copy if I think it right. I also find a template **create which is able to create a 2-d array in memory.h. So why not creat a 2-d copy instead of 1-d copy? Is there any reason?

Thank you.

Below is one function in library.cpp.

/* ----------------------------------------------------------------------

gather the named atom-based entity across all processors

name = desired quantity, e.g. x or charge

type = 0 for integer values, 1 for double values

count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f

return atom-based values in data, ordered by count, then by atom ID

e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],…

data must be pre-allocated by caller to correct length

------------------------------------------------------------------------- */

void lammps_gather_atoms(void *ptr, char *name,

int type, int count, void *data)

{

LAMMPS *lmp = (LAMMPS *) ptr;

// error if tags are not defined or not consecutive

int flag = 0;

if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;

if (lmp->atom->natoms > MAXSMALLINT) flag = 1;

if (flag && lmp->comm->me == 0) {

lmp->error->warning(FLERR,“Library error in lammps_gather_atoms”);

return;

}

int natoms = static_cast (lmp->atom->natoms);

int i,j,offset;

void *vptr = lmp->atom->extract(name);

// copy = Natom length vector of per-atom values

// use atom ID to insert each atom’s values into copy

// MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID

if (type == 0) {

int *vector = NULL;

int **array = NULL;

if (count == 1) vector = (int *) vptr;

else array = (int **) vptr;

int *copy;

lmp->memory->create(copy,count*natoms,“lib/gather:copy”);

for (i = 0; i < count*natoms; i++) copy[i] = 0;

int *tag = lmp->atom->tag;

int nlocal = lmp->atom->nlocal;

if (count == 1)

for (i = 0; i < nlocal; i++)

copy[tag[i]-1] = vector[i];

else

for (i = 0; i < nlocal; i++) {

offset = count*(tag[i]-1);

for (j = 0; j < count; j++)

copy[offset++] = array[i][j];

}

MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);

lmp->memory->destroy(copy);

} else {

double *vector = NULL;

double **array = NULL;

if (count == 1) vector = (double *) vptr;

else array = (double **) vptr;

double *copy;

lmp->memory->create(copy,count*natoms,“lib/gather:copy”);

for (i = 0; i < count*natoms; i++) copy[i] = 0.0;

int *tag = lmp->atom->tag;

int nlocal = lmp->atom->nlocal;

if (count == 1) {

for (i = 0; i < nlocal; i++)

copy[tag[i]-1] = vector[i];

} else {

for (i = 0; i < nlocal; i++) {

offset = count*(tag[i]-1);

for (j = 0; j < count; j++)

copy[offset++] = array[i][j];

}

}

MPI_Allreduce(copy,data,count*natoms,MPI_DOUBLE,MPI_SUM,lmp->world);

lmp->memory->destroy(copy);

}

}

So why not creat a 2-d copy instead of 1-d copy? Is there any reason?

It's because we wanted the library interface to be usable by languages
other than C/C++, such as Fortran. Fortran wouldn't know how
to interpret a double **array with an extra layer of pointers. But all
languages can take a pointer to a chunk of memory that contains
the data (i.e. a 1d chunk of linear memory).

That said, some of the newer functions like lammps_extract_compute
do return double ** pointers, so we're inconsistent.

Karl - is there some consistent way to have Fortran treat C pointers
as return args - like a double vs double * vs double **? Does your
Fortran wrapper work correctly with the current lammps_extract_atom
for a double **, like x for coords? Or does it have to make extra
copies of x for Fortran to use it, which kind of defeats the purpose
of that function's interface?

Steve

To answer Yumeng's question: the Fortran C bindings (from the
ISO_C_BINDING intrinsic module) declare a derived type, called "C_PTR,"
that mimics C's void* type. There is no distinction made between pointers
to different types, as they are actually treated the same way by C in
memory---the only difference is how they are dereferenced.

As I have written it now, the wrapper will copy the C array into another
Fortran array that is allocated. It sounds like I should think about
modifying it so that it translates the C pointer into a Fortran pointer,
preserving the original memory and the possibility of modifying it.
This is trivial for 1-dimensional arrays, but I'm not 100% certain
whether it will work for 2- and higher-dimensional arrays. If the array
is contiguous (e.g., was allocated with LAMMPS's memory->create), then
it should work pretty well, I would think, but I haven't tried.

If this type of functionality is something users want, I will be happy to
explore it---I haven't had all that much experience with trying to pass
arrays of pointers to doubles, say, between languages. The
lammps_extract_XXXX routines (where XXXX is not "variable") were almost an
afterthought when writing the module.

It would be helpful to know which functions are intended to be used in
this manner. I would anticipate that lammps_extract_global, for example,
would be a read-only situation. Extract_atom, on the other hand, might
not be.

-Karl

Hi, Steve,

Thank you very much for your answers. Does it mean if I do not intend to
couple Fortran code with LAMMPS I am able to use 2-D copy to get all the
data I wanted?

Best,
Yumeng Li

Hi, Karl,

Thank you so much for your reply.

I read through the file you created for coupling LAMMPS with Fortran code.
As my understanding, the Fortran interface does not able to write data back
to the actual LAMMPS data at least for now, which is what I am trying to do.
I am not pretty familiar with C. I am able to not able to write the
interface code. Because I need some new functions to library.cpp. I think I
will stay with the sample C code to couple with LAMMPS and do not need to
worry about the extra interface for Fortran. I want to let the C code to
deal with calling LAMMPS and MPI. I am wondering if I can call Fortran
codes (I am more familiar with )by this C code to do some calculations and
return some arrays and vectors to the C code. Based on some research, it
looks like it is not possible. For most of examples for calling Fortran by C
I saw, they are just passing data to Fortran not the opposite way. You are
familiar with both C and Fortran. I will appreciate if you can give me some
advises.

Best,
Yumeng Li

Fortran can call C and vice versa. What Karl has set up allows
the Fortran to get a ptr into the C data structures within LAMMPS,
so that you can do whatever you want with those values on the
Fortran side. E.g. examine them, make your own copy, change them.

Steve

I didn't understand your first message (you said "if I do not intend to
couple Fortran code with LAMMPS I am able..."; did you mean if you DO
intend?), but I can reply to the second one.

As Steve said, the module in the examples/COUPLE/fortran2 directory just
makes it easier to communicate between LAMMPS and Fortran-encoded
programs. The primary goal is to automate the processes of Fortran/C
string conversions and Fortran/C pointer conversions, which would
otherwise impede the use of the library interface. If you are trying to
add new functions (written in C or C++), I encourage you to look at
LAMMPS-wrapper.cpp and the corresponding C interfaces in LAMMPS.F90
(starting on line 68), which provide the Fortran interfaces to those
functions. You can then write Fortran wrappers, if desired, that
provide the actual Fortran functions you'll be calling.

If you give me a specific example of the data you're wishing to manipulate
from a Fortran program, I can either tell you how to do it with the
existing interface or tell you in a more specific way what you need to
do to extend the existing one to cover your needs.

Karl D. Hammond
University of Tennessee, Knoxville