using lammps as library with other c++ codes

Dear all,

A research project requires me to manipulate coordinates and forces
step by step which is better done by using lammps as library, coupling
with some of my simple c++ codes. I am looking at the given
simple.cpp, especially the 2 lines:

double *x = new double[3*natoms];
lammps_get_coords(lmp,x);

I would like to do similar things for force and potential energy, like:

double *force = new double[3*natoms];
double potE;
lammps_get_force(lmp,force);
lammps_get_potE(lmp,potE);

I have two questions.

(1)
If the potential energy can be extracted out by using the existing
"lammps_extract_variable" function, should it be used like ?
double potE = lammps_extract_variable(lmp, pe, group1);
(suppose I have defined "group1" to be the atoms that I want the pe of)
If the potential energy cannot be extracted out by using the
"lammps_extract_variable" function, what should I do ?

(2)
For the force I think I have to define a new function
"lammps_get_force" in library.cpp, I suppose it should be similar to
the definition of lammps_get_coords (copied below), but I really don't
know how and where I should change or not change to get what I want.
Could any expert please tell me how ? appreciate a lot!

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

void lammps_get_coords(void *ptr, double *coords)
{
  LAMMPS *lmp = (LAMMPS *) ptr;

  // error if tags are not defined or not consecutive

  if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) return;
  if (lmp->atom->natoms > MAXSMALLINT) return;

  int natoms = static_cast<int> (lmp->atom->natoms);
  double *copy = new double[3*natoms];
  for (int i = 0; i < 3*natoms; i++) copy[i] = 0.0;

  double **x = lmp->atom->x;
  int *tag = lmp->atom->tag;
  int nlocal = lmp->atom->nlocal;

  int id,offset;
  for (int i = 0; i < nlocal; i++) {
    id = tag[i];
    offset = 3*(id-1);
    copy[offset+0] = x[i][0];
    copy[offset+1] = x[i][1];
    copy[offset+2] = x[i][2];
  }

  MPI_Allreduce(copy,coords,3*natoms,MPI_DOUBLE,MPI_SUM,lmp->world);
  delete [] copy;
}

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

Thanks.

L.Zhang
UF Physics

Dear all,

A research project requires me to manipulate coordinates and forces
step by step which is better done by using lammps as library, coupling
with some of my simple c++ codes. I am looking at the given
simple.cpp, especially the 2 lines:

double *x = new double[3*natoms];
lammps_get_coords(lmp,x);

I would like to do similar things for force and potential energy, like:

double *force = new double[3*natoms];
double potE;
lammps_get_force(lmp,force);
lammps_get_potE(lmp,potE);

I have two questions.

(1)
If the potential energy can be extracted out by using the existing
"lammps_extract_variable" function, should it be used like ?
double potE = lammps_extract_variable(lmp, pe, group1);
(suppose I have defined "group1" to be the atoms that I want the pe of)
If the potential energy cannot be extracted out by using the
"lammps_extract_variable" function, what should I do ?

first, you should look at src/library.cpp
then, you have to define your variable in lammps script.
if you want the global potential energy, you don't need
to define a group, just an equal style variable.

overall, before playing with the library interface,
you should familiarize yourself with the features
that you want to use in regular lammps scripts
and runs.

(2)
For the force I think I have to define a new function
"lammps_get_force" in library.cpp, I suppose it should be similar to
the definition of lammps_get_coords (copied below), but I really don't
know how and where I should change or not change to get what I want.
Could any expert please tell me how ? appreciate a lot!

i don't know what to tell you. you already found the location
and have a prototype code. that means you are already 90%
done. if you don't know what to do with this at this point,
you should not be programming. just do it, man!

axel.

For the force I think I have to define a new function

no, the lib method lammps_extract_atom() already does
this, if you give it name = "f"

If the potential energy can be extracted out by using the existing
"lammps_extract_variable" function, should it be used like ?
double potE = lammps_extract_variable(lmp, pe, group1);

As Axel said, you just define a variable like

variable foo equal pe

same as you would in a LAMMPS input script

Then you can use lammps_extract_variable on foo to
get the PE of the system. But only at appropriate
times, e.g. after a run, when LAMMPS has computed
the current potential energy.

Steve

For the force I think I have to define a new function

no, the lib method lammps_extract_atom() already does
this, if you give it name = "f"

hmmm... but that would only be a pointer the force
array of the local and ghost atoms of one processor, right?

the lammps_get_coords() function collects all
coordinates from all processes with an
MPI_Allreduce() that lammps_extract_atom()
doesn't do.

axel.

Depends on whether you want the forces gathered from
all procs, or just want to manipulate them in place in
a distributed manner. I agree that you need to
add a get_forces() kind of method to the lib interface
if you want to gather them onto one proc.

Steve

Just added more general methods to the LAMMPS library interface
so that any atom property can be gathered across all procs or scattered
cak to all procs. It should work with forces. See the Aug 18 patch.

Steve