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 (mask[i] & groupbit) {
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