Hi steve and developers,
Recently, I add a new fix command via modifying the fix_addforce.cpp. The function of my command is that when I add a force on i atom, I want to add a counterforce on a neighbor atom via randomly choosing an atom in i’s neighbor list. To achieve the function, I use a full neighbor list in init() subfunction via adding the following red codes:
/* ---------------------------------------------------------------------- */
void FixActDPD::init()
{
// need a full neighbor list, built whenever re-neighboring occurs
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}
Then, I modify code of post_force(int vflag) subfunction:
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
if (varflag == CONSTANT) {
double unwrap[3];
for (int ii = 0; ii < nlocal; ii++){
i = ilist[ii];
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
domain->unmap(x[i],image[i],unwrap);
////////////////////////////////////////////////////////////////////
////calcuting force in dipole
///d_mu/dt = omega cross mu
///renormalize mu to unit dipole
if(mu[i][3] > 1.0e-30){
//printf(“tag=%d\n”,tag[i]);
ka =0,k=0;
jlist = firstneigh[i];
jnum = numneigh[i];
for(int jj = 0; jj < jnum; jj++){
j = jlist[jj];
j &= NEIGHMASK;
if( mask[j] & groupbit ) {
k++;
if( mu[j][3] < 1.0e-30 ) ka++;
}
}
if(k==0) continue;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
xvalue = pvalue * mu[i][0]/mu[i][3];
yvalue = pvalue * mu[i][1]/mu[i][3];
zvalue = pvalue * mu[i][2]/mu[i][3];
///////////////////////////////////////
foriginal[0] -= xvalueunwrap[0] + yvalueunwrap[1] + zvalue*unwrap[2];
foriginal[1] += f[i][0];
foriginal[2] += f[i][1];
foriginal[3] += f[i][2];
f[i][0] += xvalue;
f[i][1] += yvalue;
f[i][2] += zvalue;
////add conterforce on water atoms
if(ka!=0){
while(1){
int jj = (int)floor(jnumrandom->uniform());
j = jlist[jj];
j &= NEIGHMASK;
//random choosing atom
if( mu[j][3] < 1.0e-30 ){
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = sqrt(delxdelx + delydely + delzdelz);
if(rsq < 1.0e-30) error->all(FLERR,“overlap”);
delx /= rsq;
dely /= rsq;
delz /= rsq;
if((delxmu[i][0] + delymu[i][1] + delzmu[i][2])/mu[i][3] <= -cosangle) break;
}
}
}
else{
while(1){
int jj = (int)floor(jnumrandom->uniform());
j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = sqrt(delxdelx + delydely + delzdelz);
if(rsq < 1.0e-30) error->all(FLERR,“overlap”);
delx /= rsq;
dely /= rsq;
delz /= rsq;
if((delxmu[i][0] + delymu[i][1] + delzmu[i][2])/mu[i][3] <= -cosangle) break;
}
}
f[j][0] -= xvalue;
f[j][1] -= yvalue;
f[j][2] -= zvalue;
}
}
However, I can not get the result what I want to obtain when using the code.
The questions I want to ask are: 1) Do it need rebuild the neighbor list when using full neigborlist? 2) When I add force on i atom and couterforce on one of its neighbor atoms, Do it need add code for communication in subfunction?
The fix command is attached.
Any suggestion is appreciated. Thanks in advance
fix_actdpd2.0.cpp (12.6 KB)