Zero output from fix ave/atom

LAMMPS LAMMPS Development
1

Hi everyone,

I am using LAMMPS to evaluate the decay of NMR signal in magnetic field gradient. I have written a “compute” that evaluates the accumulated phase of each particle after each time step. I dump the phase after a set number of timesteps. While this works with no problems on one node, and for shorter simulations, at longer evolutions the output phase for some of the particles becomes zero.

Here is how I prepare my output:

compute 	NMRsig all NMRsig ${omegax} ${omegay} ${omegaz} 
fix             signalf all ave/atom 1 1 1 c_NMRsig[*]	
dump		NMRsig_out all custom 10000000 nmrsig200_1.txt id f_signalf[*]
dump_modify     NMRsig_out sort id
#number of steps to run
run_style	verlet
run		1000000000

And here is a set of output

ITEM: TIMESTEP
10000000
ITEM: NUMBER OF ATOMS
1000
ITEM: BOX BOUNDS pp pp pp
-1.7000000000000000e+02 1.7000000000000000e+02
-1.7000000000000000e+02 1.7000000000000000e+02
-1.7000000000000000e+02 1.7000000000000000e+02
ITEM: ATOMS id f_signalf[1] f_signalf[2] f_signalf[3] f_signalf[4] f_signalf[5] f_signalf[6]
…
550 1 -3.55613e-05 1 1.51378e-05 1 -3.99799e-06
551 1 2.68159e-06 1 -2.86324e-05 1 -7.39508e-06
552 1 3.20163e-05 1 -2.91998e-05 1 3.03028e-05
553 1 1.48324e-05 1 -3.1044e-05 1 1.81926e-05
554 1 1.89969e-05 1 -1.43878e-05 1 1.79598e-05
555 1 -2.38021e-05 1 6.45593e-05 1 -3.14228e-05
556 0 0 0 0 0 0
557 1 -2.5495e-05 1 -3.21013e-05 1 -8.20196e-06
558 1 -1.02374e-05 1 -1.53342e-05 1 4.43111e-05
559 1 -4.53915e-06 1 -3.78215e-07 1 1.11562e-05
560 1 -1.44737e-05 1 -1.52298e-05 1 -4.23223e-05
561 1 3.33534e-05 1 5.21771e-05 1 6.90638e-06
562 1 -8.80949e-06 1 -9.83268e-06 1 2.66596e-06
563 1 -2.12596e-05 1 -1.92553e-05 1 -2.27869e-05
564 1 -5.033e-06 1 2.88323e-05 1 1.61323e-05
565 1 -3.56757e-06 1 2.81538e-05 1 1.95783e-06
566 1 -3.42753e-05 1 -1.31754e-05 1 1.56249e-05
567 1 4.85273e-05 1 -7.47894e-05 1 -8.89502e-06
568 0 0 0 0 0 0
569 1 -2.22802e-05 1 1.06457e-05 1 7.42265e-06
570 1 3.90405e-05 1 -1.16549e-05 1 -3.92964e-05
571 1 -2.32027e-05 1 2.17756e-05 1 -4.24659e-06
572 1 2.69602e-05 1 -6.13531e-06 1 4.6356e-05
573 1 -3.37474e-05 1 -2.0135e-05 1 -8.65166e-06
574 1 -4.48482e-05 1 0.000130097 1 -1.80315e-05
575 1 -2.99853e-05 1 -2.1133e-05 1 1.50947e-05
576 1 1.3288e-05 1 -9.63368e-06 1 2.16961e-05
577 1 3.22323e-07 1 3.24864e-05 1 -2.25511e-05
578 1 -1.12512e-05 1 2.51472e-05 1 2.46447e-05
579 1 2.82556e-05 1 -1.92971e-05 1 1.15178e-05
580 1 3.33253e-05 1 -1.36951e-05 1 2.15541e-05
581 1 4.26993e-06 1 -4.16841e-05 1 5.44215e-05
582 1 -8.09298e-06 1 -5.64074e-05 1 -1.69614e-05
583 1 -2.15552e-05 1 1.30414e-06 1 4.55563e-05
584 1 5.98638e-06 1 -2.47841e-05 1 -3.21058e-05
585 1 5.45646e-05 1 1.00233e-05 1 3.5604e-06
586 0 0 0 0 0 0
587 1 -8.23416e-06 1 2.84374e-05 1 -1.42745e-05
588 1 2.23211e-05 1 -2.18053e-05 1 -3.46998e-06
589 1 -7.45714e-06 1 -9.49005e-06 1 1.7271e-05
590 1 -1.46878e-05 1 4.53119e-06 1 -7.63254e-07

The location of zero output is different at each time step.
Any idea what could be the cause and how to avoid this?

2

Please always report which LAMMPS version you are working with.

Well, the most obvious source of possible errors is your new code. Without carefully reviewing that code, it is difficult to make any statements about what the problem could be.

It is not clear why to pass this through fix ave/atom. Why not pass the compute reference to the dump directly?

3

Thanks for getting back to me.
The LAMMPS version is:
LAMMPS (29 Sep 2021 - Update 1)

The reason for using the fix/dump combination is that I need the compute to be called at every time step, but the dump file is written once every 10^7 time steps. Is there a way to do this using only a dump?

Here is the relevant part of the code from mycompute.cpp. I unwrap the coordinates, calculate phase increments along x,y,z and multiply them to the accumulated phase from previous time steps. The six components are real and imaginary parts of the three phases.

void ComputeNMRSig::compute_peratom()
{
invoked_peratom = update->ntimestep;

// grow local accphase array if necessary

if (atom->nmax > nmax) {
memory->destroy(accphase);
nmax = atom->nmax;
memory->create(accphase,nmax,6,“nmrsig/atom:accphase”);
array_atom = accphase;
}

// xt,yt,zt = unwrapped position of atom
// original unwrapped position is stored by fix
// for triclinic, need to unwrap current atom coord via h matrix

const std::complex I(0, 1);

double **x = atom->x;
int *mask = atom->mask;
imageint *image = atom->image;
int nlocal = atom->nlocal;

double *h = domain->h;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;

int xbox,ybox,zbox;
double xt,yt,zt;

if (0 == update->ntimestep % 1000000000 )
{ for (int i = 0; i < nlocal; i++){
accphase[i][0] = 1.0;
accphase[i][1] = 0.0;
accphase[i][2] = 1.0;
accphase[i][3] = 0.0;
accphase[i][4] = 1.0;
accphase[i][5] = 0.0;
}
}

if (domain->triclinic == 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xt = x[i][0] + xboxxprd;
yt = x[i][1] + yboxyprd;
zt = x[i][2] + zbox*zprd;
std::complex phx(accphase[i][0] ,accphase[i][1]);
std::complex phy(accphase[i][2] ,accphase[i][3]);
std::complex phz(accphase[i][4] ,accphase[i][5]);
accphase[i][0] = std::real(std::exp(I * omegax * xt) * phx);
accphase[i][1] = std::imag(std::exp(I * omegax * xt) * phx);
accphase[i][2] = std::real(std::exp(I * omegay * yt) * phy);
accphase[i][3] = std::imag(std::exp(I * omegay * yt) * phy);
accphase[i][4] = std::real(std::exp(I * omegaz * zt) * phz);
accphase[i][5] = std::imag(std::exp(I * omegaz * zt) * phz);
}
else {accphase[i][0] = accphase[i][2] = accphase[i][4] = 1.0;
accphase[i][1] = accphase[i][3] = accphase[i][5] = 0.0;
}

} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xt = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
yt = x[i][1] + h[1]*ybox + h[3]*zbox;
zt = x[i][2] + h[2]*zbox;
std::complex phx(accphase[i][0] ,accphase[i][1]);
std::complex phy(accphase[i][2] ,accphase[i][3]);
std::complex phz(accphase[i][4] ,accphase[i][5]);
accphase[i][0] = std::real(std::exp(I * omegax * xt) * phx);
accphase[i][1] = std::imag(std::exp(I * omegax * xt) * phx);
accphase[i][2] = std::real(std::exp(I * omegay * yt) * phy);
accphase[i][3] = std::imag(std::exp(I * omegay * yt) * phy);
accphase[i][4] = std::real(std::exp(I * omegaz * zt) * phz);
accphase[i][5] = std::imag(std::exp(I * omegaz * zt) * phz);
}
else {accphase[i][0] = accphase[i][2] = accphase[i][4] = 1.0;
accphase[i][1] = accphase[i][3] = accphase[i][5] = 0.0;
}
}
}

4

Please note that there has been a new stable release a couple of months ago. And for developing new features, we generally recommend to work of the develop branch (and stay up-to-date with it) since we have ongoing refactoring efforts of internal code and thus porting from older versions can become very tedious while staying current makes those issues incremental and it will be easier to get your code integrated.

If this needs to be called every step (or otherwise regularly), then it should be a fix and not a compute.
Compute styles are “on demand”, Fix styles decide on their own how often and when in the time step loop they get executed.

Fixes can also migrate per-atom data with the atoms when they move between subdomains.
A compute has no way to do that. Computes that need to keep per-atom data around - like compute msd - delegate that work to a(n internal) fix STORE/PERATOM instance.