# [lammps-users] how to debug self written command

Pieter, Thanks for sharing your experience.

Thanks, Guillaume. Actually, I want to write a fix_force command to
relax the velocity of atoms to a specific value. so i based my
fix_couple_force.cpp on fix_ave_force.cpp, and in the function
<post_force>, I averaged the atoms' velocity on a 2*2 grid, then i
want to change the force of every atom in the grids to let it relax to
some specific grid value. but i am confused by the result of the
following code, which want to average the number of density in a 2*2
grid:

nx=2;nz=2;
delta_x=(domain->boxhi[0]-domain->boxlo[0])/nx;
invdelta_x=1./delta_x;
invdelta_z=1./delta_z;

double *vectorv,*vectorf;
int nstride = 3;
vectorv = &atom->v[0][0];

m = 0;
for (i = 0; i < nlocal; i++) {
if(x[i][2]<=ztop){
iz= static_cast<int> ((ztop-x[i][2])*invdelta_z);
if(iz<nz){
ix=static_cast<int>((x[i][0]-domain->boxlo[0])*invdelta_x);
igrid=iz*nx+ix;
count_one[igrid]+=1.0;
v_grid_one[igrid][0]+=vectorv[m];
}
}

}
m += nstride;
}

after some try-error process, i know LAMMPS will hang on because ix
will get a value of 2 sometimes, i don't know why this would happen
considering the chosen delta_x, am i wrong with the domain->boxlo[0]??
if i insert the following code, the bug is gone:
if(ix<0) ix=0;
if (ix>1) ix=1;

Another question: what's the execution order of the following fix
subroutines? is it going like this:
init - setup - initial_integrate - (calculate force) - post_force -
(verlet integration) - final_integrate - end_of_step

Thanks!

Yin-chun

Yin-chun,

At cetain points during the timestep, particles can move beyond the
box bounds. Looks to me like that's why you're getting a value of 2
for ix sometimes. You might try inserting your fix at the
"pre_neighbor" stage rather than at the "post_force" stage.

And I agree with the order of execution of the statements you have:

init - setup - initial_integrate - (calculate force) - post_force -
(verlet integration) - final_integrate - end_of_step

But I'd leave out "verlet integration" because LAMMPS treats the
verlet integration like a fix (see fix nve) --- it happens in
intial_integrate and final_integrate. Also note that "init" and
"setup" do not happen every timestep --- they occur before entering
the time stepping loop (see Verlet::iterate).

Paul