Making a new compute which uses a custom per-atom vector

Dear LAMMPS users.

Hi,

I am trying to make a new compute, this compute calculates a global scalar value based on a custom per-atom vector which is created by ‘fix property/atom’.

I believe that my code works well regarding finding and accessing a custom per-atom vector (by using atom->find_custom, atom->dvector) but there is something wrong about updating values.

Below is my test run result and ‘c_elecThermalEnergy’ is my new compute. The value should be ~ 28.61 but sometimes the value is suddenly changed (completely wrong), even NaN.

In addition to that, sometimes the value is just wrong at the start.

I am not sure what is the source of such misbehavior. Because fundamentally my compute just did a similar thing, I compared my new codes with compute_ke.cpp & compute_ke.h, but looks like that my code is just fine.

I am not fully understanding every details of LAMMPS modification, but I believe that I missed something about ‘flag’ or update routine.

Below is my new compute constructor,

ComputeElecThermalEnergy::ComputeElecThermalEnergy(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{

int iarg = 0;
int flag = 0;
int index_custom = 0;


// !! ARGUMENT !! compute ID group-ID eThermalEnergy d_ID lambda ${lambda} V ${volumePerAtom}
// narg = # of arguments (counting from ID to the end)
// arg = arguments (strings); start from 0 index (ID) to

if (narg != 8) error->all(FLERR,"Illegal compute eThermalEnergy command");  // !!ERROR!! not enough arguments

// force->numeric does interpert a string as a floating number by using 'atof'.
// strcmp; compare strings.
/*
int iarg = 4;
iarg++;
*/

iarg = 3;
if (strstr(arg[iarg],"d_") == arg[iarg]) { // check the name of this argument started with 'd_',
    
    index_custom = atom->find_custom(&arg[iarg][2],flag);
    if (index_custom < 0 || flag != 1) error->all(FLERR,"Illegal compute eThermalEnergy command; the electric temperature per-atom vector does not exist.");
    
    temp = atom->dvector[index_custom];
}else{
    error->all(FLERR,"Illegal compute eThermalEnergy command; the electric temperature per-atom vector should be float.");
}

        
if (strcmp(arg[4],"lambda") == 0){
    lambda = force->numeric(FLERR,arg[5]);
    if (lambda <= 0){ error->all(FLERR,"Illegal compute eThermalEnergy command; lambda should be greater than 0"); }
}else{ error->all(FLERR,"Illegal compute eThermalEnergy command"); }

if (strcmp(arg[6],"V") == 0){
    V = force->numeric(FLERR,arg[7]);
    if (V <= 0){ error->all(FLERR,"Illegal compute eThermalEnergy command; V should be greater than 0"); }
}else{ error->all(FLERR,"Illegal compute eThermalEnergy command"); }


// Flags
scalar_flag = 1;
extscalar = 1;
timeflag = 1;
//comm_reverse = 1;

}

Below is a compute_scalar method of my new compute,

double ComputeElecThermalEnergy::compute_scalar()
{
int i;
int nlocal = atom->nlocal;
int *mask = atom->mask;
double domainSum = 0.0;

/* Probably these are unnecessary parts.
int *mask = atom->mask;
bigint natoms = atom->natoms;
 */

// Don't know why, these below lines should be to make this code work.
invoked_scalar = update->ntimestep;
if (update->eflag_global != invoked_scalar)
    error->all(FLERR,"~ was not tallied on needed timestep");

 /* Probably these are unnecessary parts.
// grow local energy array if necessary
// needs to be atom->nmax in length

if (atom->nmax > nmax) {
memory->destroy(energy);
nmax = atom->nmax;
memory->create(energy,nmax,"ee:energy");
vector_atom = energy;
}
*/

// Calculate 'domainSum'.
for (i = 0; i < nlocal; i++){
    if (mask[i] & groupbit){
        domainSum += (temp[i]*temp[i]);
    }
}

// Calculate 'globalSum'.
MPI_Allreduce(&domainSum,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
scalar *= 0.5*lambda*V;

return scalar;
}

If you have any comment, it would be really appreciated.

Youhwan

How is this property initialized?

It would be easier to have the entire code and a suitable example input deck.
Without I just have to guess. It looks like you are looking up a pointer in the constructor and are not updating it. However, atoms will migrate between processors and thus also the per-atom arrays will change. Thus you would need to look those pointer up every time.

It is good practice to check memory accesses with valgrind. Your code should be valgrind clean (outside of issues triggered by OpenMPI or OpenMP, thus it is recommended to do this with MPICH or a serial executable and without OpenMP enabled (you can still enable the OPENMP package).

Thank you for your comment.

The property is initialized in my LAMMPS script. For example,

fix 1 all property/atom d_elecTemp
compute 1 all property/atom d_elecTemp
set group all d_elecTemp ${initialValue}

And my new compute identifies a custom property which it should use, by getting a property name as an argument.

compute elecThermalEnergy all eThermalEnergy d_elecTemp lambda v_lambda V v_volumePerAtom

I will definitely check them with valgrind.

When quoting input, you should wrap those text blocks into “```” to have them properly typeset.

That is not addressing my concern. You do look it up correctly, but do you look it up every time you access it? As I mentioned the memory addresses change on every reneighbor when atoms can migrate between processors (or from one end of a box to the other due to PBC).

Thank you so much.

I misunderstood your comment, I never thought about that the pointer address can be changed. The problem is now solved.

Indeed, in my code, the pointer is not changed, just fixed after the constructor is called.

I put a single line in my code thus for every time compute_scalar() is called, renews the pointer.

double ComputeElecThermalEnergy::compute_scalar()
{
int i;
int nlocal = atom->nlocal;
int *mask = atom->mask;
double domainSum = 0.0;

/* Probably these are unnecessary parts.
int *mask = atom->mask;
bigint natoms = atom->natoms;
 */

invoked_peratom = update->ntimestep;
// Don't know why, these below lines should be to make this code work.
invoked_scalar = update->ntimestep;
if (update->eflag_global != invoked_scalar)
    error->all(FLERR,"~ was not tallied on needed timestep");

 /* Probably these are unnecessary parts.
// grow local energy array if necessary
// needs to be atom->nmax in length

if (atom->nmax > nmax) {
memory->destroy(energy);
nmax = atom->nmax;
memory->create(energy,nmax,"ee:energy");
vector_atom = energy;
}
*/

//
temp = atom->dvector[index_custom];


// Calculate 'domainSum'.
for (i = 0; i < nlocal; i++){
    if (mask[i] & groupbit){
        domainSum += (temp[i]*temp[i]);
    }
}

// Calculate 'globalSum'.
MPI_Allreduce(&domainSum,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
scalar *= 0.5*lambda*V;

return scalar;
}

I do have following questions though, actually, I also made a new fix. This new fix also uses the same custom per-atom vector. In this fix, the pointer indicating the custom vector is not changed at all just like the case of my previous compute code.

Nonetheless, I did not find any problem in that fix. Could you give me a short clue why my fix is just OK without redefining the pointer?

This fix using post_force() method to do necessary calculations.

Thank you
Youhwan

Sorry, but I cannot read minds, so without knowing what exact the the fix does, how and when, I cannot say anything. It may still be wrong, but just not causing a problem this far.

Again, valgrind is the tool of choice to validate whether memory accesses are valid.

1 Like