Using fix external gives different results from using addforce (per atom variable)

LAMMPS LAMMPS Development
1

Hello!
I am trying to apply an active force to particles. I tried to do this using fix external and filling the fexternal array but in this way I do not get the correct results. I tried another way, still using the fix external, where I, instead of using the fexternal array, gathered the force vector from all particles and simply added the active force directly and then scattered. This works!
Here are the two implementations. I have no idea why they yield different results. i checked and the fix is called at post_force. That’s why directly messing with the force array works, it is not cleared before the forces are applied.
Any ideas?

void active_force_fails(void *ptr, long int timestep, int nlocal, int *ids, double **x, double **fexternal)
{
double *theta = (double )malloc(nlocalsizeof(double));
lammps_gather_atoms(ptr,“d_theta”,1,1,theta);

double *dt_ptr = (double *)lammps_extract_global(ptr,"dt");
double dt = *dt_ptr;

double *drot_ptr = (double *) lammps_extract_variable(ptr,"Drot", NULL);
double drot = *drot_ptr;

double *factive_ptr = (double *) lammps_extract_variable(ptr,"factive", NULL);
double factive = *factive_ptr;

double upd;
for(int i=0;i<nlocal;i++) {
	fexternal[i][0] = factive*cos(theta[i]);
	fexternal[i][1] = factive*sin(theta[i]);
            fexternal[i][2] = 0.0;
	upd=sqrt(2*drot*dt)*gaussian_random(0.0,1.0);
	theta[i] += upd;
}

lammps_scatter_atoms(ptr,"d_theta",1,1,theta);

free(theta);

}

void active_force_works(void *ptr, long int timestep, int nlocal, int *ids, double **x, double **fexternal)
{
double *theta = (double )malloc(nlocalsizeof(double));
lammps_gather_atoms(ptr,“d_theta”,1,1,theta);

double *dt_ptr = (double *)lammps_extract_global(ptr,"dt");
double dt = *dt_ptr;

double *drot_ptr = (double *) lammps_extract_variable(ptr,"Drot", NULL);
double drot = *drot_ptr;

double *factive_ptr = (double *) lammps_extract_variable(ptr,"factive", NULL);
double factive = *factive_ptr;
double *f = (double *)malloc(3*nlocal*sizeof(double));
lammps_gather_atoms(ptr,(char *) "f",1,3,f);

double upd,sqr2dt;
sqr2dt = sqrt(2*drot*dt);
for(int i=0;i<nlocal;i++) {
	f[i*3]   += factive*cos(theta[i]);
	f[3*i+1] += factive*sin(theta[i]);
	upd=sqr2dt*gaussian_random(0.0,1.0);
	theta[i] += upd;
}
lammps_scatter_atoms(ptr,"d_theta",1,1,theta);
lammps_scatter_atoms(ptr,(char *) "f",1,3,f);

free(f);
free(theta);

}

2

The code you quote does not look correct to me in both cases. It is difficult to tell, since the rest (input and how it is hooked up with LAMMPS) is missing.

For starters, you are allocating arrays of the size “nlocal” but then are using gather and scatter operations on them. That makes no sense, since those arrays should be of the size “natoms”.
This should crash unless you have only one MPI rank and thus nlocal == natoms.

Then you are mixing gathered data (which is sorted by atom ID) with local arrays which are indexed differently. No surprise that things don’t add up.

The next unexpected thing is that your computation seems entirely based on data that is available as variables? Why try to use fix external in the first place? Why not just do the computation with equal and atom style variables and then use fix addforce.

I think that you get the expected result (which you call “correct”) only because of calling the scatter function on the force array, which is definitely incorrect and should lead to inconsistent data (except for the added force).

3

Dear Axel,
I am very thankful for your comment.
Indeed I was wrongly dealing with indexing…your comment made me dive into how indexing works. Now using the atom->tag property i have correctly made the correspondence between my local arrays and the particles and my fix is working properly.
Sorry for the delay in replying but your comment allowed me to finish my fix which i have tested for the last weeks.
Once again thanks!

4

Thanks for reporting back. This will be useful for future readers of this discussion.