C++ interface: unable to delete array after casting property/atom via lammps_extract_compute

LAMMPS LAMMPS Mailing List Mirror
1

Hello,
I have been developing a C++ code using the LAMMPS interface. I have recently noticed memory leakage issues and started cleaning the code up (various arrays created with “new double” and then not deleted []…).

I noticed that I am unable to delete [] any array to which I casted property/atom via lammps_extract_compute.

I have downloaded the latest version of LAMMPS (7Aug19) and edited the example/COUPLE/simple/simple.cpp to reproduce the error. I have copied the edited main below, highlighting between ************ any lines that I have added. Everything works fine and I can read and print atom IDs using the compute, but when I get to
delete [] aID
at the very bottom, the code crashes with Abort trap: 6 (signal 6)

I am concerned because I am using this approach extensively in my code, creating such aIDs (and others for radii, type, energy, stress per atom) in a function that is called in a loop for thousands of times. If I do not delete them, I am afraid this could cause memory leakage. Notice that any other array that I create and manipulate (as long as not casting from lammps) can be deleted without any problem, as expected.

Last note, unrelated but since I’m here… If you can suggest a better way to evaluate a property/atom compute rather than my approach of creating, using, and then overwriting a temporary dump every time, I’d be very interested (let me know if I should open a separate thread for this).

Thank you,
Enrico

int main(int narg, char **arg)
{
// setup MPI and various communicators
// driver runs on all procs in MPI_COMM_WORLD
// comm_lammps only has 1st P procs (could be all or any subset)

MPI_Init(&narg,&arg);

if (narg != 3) {
printf(“Syntax: simpleCC P in.lammps\n”);
exit(1);
}

int me,nprocs;
MPI_Comm_rank(MPI_COMM_WORLD,&me);
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);

int nprocs_lammps = atoi(arg[1]);
if (nprocs_lammps > nprocs) {
if (me == 0)
printf(“ERROR: LAMMPS cannot use more procs than available\n”);
MPI_Abort(MPI_COMM_WORLD,1);
}

int lammps;
if (me < nprocs_lammps) lammps = 1;
else lammps = MPI_UNDEFINED;
MPI_Comm comm_lammps;
MPI_Comm_split(MPI_COMM_WORLD,lammps,0,&comm_lammps);

// open LAMMPS input script

FILE *fp;
if (me == 0) {
fp = fopen(arg[2],“r”);
if (fp == NULL) {
printf(“ERROR: Could not open LAMMPS input script\n”);
MPI_Abort(MPI_COMM_WORLD,1);
}
}

// run the input script thru LAMMPS one line at a time until end-of-file
// driver proc 0 reads a line, Bcasts it to all procs
// (could just send it to proc 0 of comm_lammps and let it Bcast)
// all LAMMPS procs call input->one() on the line

LAMMPS *lmp = NULL;
if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps);

int n;
char line[1024];
while (1) {
if (me == 0) {
if (fgets(line,1024,fp) == NULL) n = 0;
else n = strlen(line) + 1;
if (n == 0) fclose(fp);
}
MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
if (n == 0) break;
MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
if (lammps == 1) lammps_command(lmp,line);
}

// run 10 more steps
// get coords from LAMMPS
// change coords of 1st atom
// put coords back into LAMMPS
// run a single step with changed coords

double *x = NULL;
double *v = NULL;

/************************
/************************
/************************
double aID = NULL;
/***********************
/************************
/************************

if (lammps == 1) {
lmp->input->one(“run 10”);

int natoms = static_cast (lmp->atom->natoms);
x = new double[3natoms];
v = new double[3natoms];
lammps_gather_atoms(lmp,(char *) “x”,1,3,x);
lammps_gather_atoms(lmp,(char *) “v”,1,3,v);
double epsilon = 0.1;
x[0] += epsilon;
lammps_scatter_atoms(lmp,(char *) “x”,1,3,x);

// these 2 lines are the same

// lammps_command(lmp,“run 1”);

/************************

/************************
/************************
lmp->input->one(“compute tempID all property/atom id”);
lmp->input->one(“dump tdID all custom 1 dump.temp c_tempID”);
aID = new double[natoms]; // array of unsorted IDs
aID = ((double *) lammps_extract_compute(lmp,(char *) “tempID”,1,1));

lmp->input->one(“run 1”);

lmp->input->one(“undump tdID”);
lmp->input->one(“uncompute tempID”);

for (int i=0; i<natoms; i++) printf("%d\n",(int)aID[i]);
/************************
/************************
/************************

}

// extract force on single atom two different ways

if (lammps == 1) {
double **f = (double **) lammps_extract_atom(lmp,(char *) “f”);
printf(“Force on 1 atom via extract_atom: %g\n”,f[0][0]);

double *fx = (double *)
lammps_extract_variable(lmp,(char *) “fx”,(char *) “all”);
printf(“Force on 1 atom via extract_variable: %g\n”,fx[0]);
}

// use commands_string() and commands_list() to invoke more commands

char *strtwo = (char *) “run 10\nrun 20”;
if (lammps == 1) lammps_commands_string(lmp,strtwo);

char *cmds[2];
cmds[0] = (char *) “run 10”;
cmds[1] = (char *) “run 20”;
if (lammps == 1) lammps_commands_list(lmp,2,cmds);

// delete all atoms
// create_atoms() to create new ones with old coords, vels
// initial thermo should be same as step 20

int *type = NULL;

if (lammps == 1) {
int natoms = static_cast (lmp->atom->natoms);
type = new int[natoms];
for (int i = 0; i < natoms; i++) type[i] = 1;

lmp->input->one(“delete_atoms group all”);
lammps_create_atoms(lmp,natoms,NULL,type,x,v,NULL,0);
lmp->input->one(“run 10”);
}

delete [] x;
delete [] v;
delete [] type;

/************************
/************************
/************************
delete [] aID;
/************************
/************************
/************************

// close down LAMMPS

delete lmp;

// close down MPI

if (lammps == 1) MPI_Comm_free(&comm_lammps);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
}

2

Hello,
I have been developing a C++ code using the LAMMPS interface. I have recently noticed memory leakage issues and started cleaning the code up (various arrays created with “new double” and then not deleted []…).

I noticed that I am unable to delete [] any array to which I casted property/atom via lammps_extract_compute.

you must not delete those. those are pointers to storage managed by LAMMPS.

I have downloaded the latest version of LAMMPS (7Aug19) and edited the example/COUPLE/simple/simple.cpp to reproduce the error. I have copied the edited main below, highlighting between ************ any lines that I have added. Everything works fine and I can read and print atom IDs using the compute, but when I get to
delete [] aID
at the very bottom, the code crashes with Abort trap: 6 (signal 6)

for a good reason. you should not use delete[] on this data. in fact, your code accessing the data is faulty. you must not use “uncompute” until you are done with the data to the pointer you have extracted from the compute. it only works in the case of atom IDs, because this is not data managed by compute property/atom directly but elsewhere and compute property/atom merely forwards the pointer to it.

I am concerned because I am using this approach extensively in my code, creating such aIDs (and others for radii, type, energy, stress per atom) in a function that is called in a loop for thousands of times. If I do not delete them, I am afraid this could cause memory leakage. Notice that any other array that I create and manipulate (as long as not casting from lammps) can be deleted without any problem, as expected.

Last note, unrelated but since I’m here… If you can suggest a better way to evaluate a property/atom compute rather than my approach of creating, using, and then overwriting a temporary dump every time, I’d be very interested (let me know if I should open a separate thread for this).

for atom ids, and plenty other per-atom properties you could use lammps_extract_atom() which doesn’t need something to trigger the compute.

axel.

3

Thank you Axel,

so, if I understand correctly, because the memory is managed by LAMMPS, there would be no memory leakage if I were to:

  1. extract_compute my array aID and use its content
  2. then set aiD = NULL
  3. change the number of atoms in LAMMPS (delete_atoms or fix_deposit)
  4. allocate my array aID again using new double[natoms], and then back to point 1 above in a loop?

Of course, assuming I a not making other mistakes, i.e. properly managing computes (unlike my previous example) and updating natoms etc.

Thanks for the hint about extract_atoms. Months ago I started with that indeed, but then I needed extract_compute to get energy per atom, and I was unsure whether the ordering of ids from extract_atoms would have matched the ordering of energies from extract_compute. So, since I had to compute energies anyhow, I opted for extract_compute-ing also the ids, types, etc.

Enrico

4

Thank you Axel,

so, if I understand correctly, because the memory is managed by LAMMPS, there would be no memory leakage if I were to:

  1. extract_compute my array aID and use its content
  2. then set aiD = NULL
  3. change the number of atoms in LAMMPS (delete_atoms or fix_deposit)
  4. allocate my array aID again using new double[natoms], and then back to point 1 above in a loop?

no. you don’t allocate aiD anywhere. lammps_extract_atom() and lammps_extract_compute() with compute property/atom would return the same direct pointer to the internal storage managed by LAMMPS.
if you do allocate aID and then use it with lammps_extract_atom() or lammps_extract_compute(), you will leak memory, since that allocated data is lost.

Of course, assuming I a not making other mistakes, i.e. properly managing computes (unlike my previous example) and updating natoms etc.

Thanks for the hint about extract_atoms. Months ago I started with that indeed, but then I needed extract_compute to get energy per atom, and I was unsure whether the ordering of ids from extract_atoms would have matched the ordering of energies from extract_compute. So, since I had to compute energies anyhow, I opted for extract_compute-ing also the ids, types, etc.

the order for domain distributed per-atom data is always the same and stretches from 0 to atom->nlocal-1 (plus from atom->nlocal to atom->nlocal+atom->nghost-1, if that property is present for ghost atoms as well).
compute property/atom does not compute anything, it just provides access to internal data, that is not otherwise always directly available inside of LAMMPS scripts. before this compute existed, all fixes had to be programmed to access data directly and it was becoming a maintenance problem. when using the library interface, lammps_extract_atom() is the direct route and preferred. that said, you will have to renew the pointers at every timestep, since those arrays may have changed during re-neighboring and redistributing atoms into sub-domains in LAMMPS.

axel.

5

Thanks, now I have understood how to use it (without ever allocating memory to my array):

double *aID;
aID = lammps_extract_compute(…)
do stuff with it
change number of atoms or else with lammps as needed
aID = lammps_extract_compute(…)
etc…

I have written the above just for future reference for whom might look this up in the mailing list.

Thank you also for the explanation of how the various extract_ have been conceived and some of the logics behind them. I will consider again using extract_atoms when possible.

Enrico