Hi everyone,
My atom style is hybrid bond and dipole, and I need average the values of atom->mu of atoms in each molecule. And then I use these average value as input values to solve the question dmu_i/dt = -k for each atom. My idea is that I use atom->omega to store the input value of dipole atoms.
Then I try to modify the atom_vec_dipole.h/cpp, bond_harmonic.h/cpp, fix nve.h/cpp to achieve my goal.
I tried to open the atom_omega_flag in the atom_vec_dipole.cpp/h files, like the atom_vec_line.h/cpp. It does work.
Then I try to modify the bond_harmonic.h/cpp (to bond_cell.h/cpp) and calculate the average by traversing the bond list and then assign the to each atom on the same molecule via storing them in atom->omega( the key code is given below).
The final step is that I try to use the atom->omega in the modified fix nve.h/cpp files(fix_nve_cell.h/cpp) via **omega = atom->omega;
However, unexpectedly, the atom->omega values used in fix_nve_cell.cpp are not changed by bond_cell.h/cpp.
All output omegas are zero in fix_nve_cell.cpp.
I do not know why it happens. I try to use com_modify vel yes to communicate the omega. It does not work.
I think the atom->omega is like atom->f ( bond style could change f ) or atom->v.
Any idea is appreciated to help me deal with the problem.
Best wishes.
/* ---------------------------------------------------------------------- */
//
void BondCell::compute(int eflag, int vflag)
{
int i1,icell,n,type;
tagint *molecule = atom->molecule ;
double **omega = atom->omega;
double **mu = atom->mu;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
// calculating the average for each molecule
for (icell = 0; icell < ncell; icell ++)
avomega[icell][0] = avomega[icell][1] = avomega[icell][2] = 0.0;
for (n = 0; n < nbondlist; n++) {
i1 = bondlist[n][0];
avomega[mol2cell[molecule[i1]]][0] += mu[i1][0];
avomega[mol2cell[molecule[i1]]][1] += mu[i1][1];
avomega[mol2cell[molecule[i1]]][2] += mu[i1][2];
}
// average omega
for (icell = 0; icell < ncell; icell++) {
avomega[icell][0] /= nCELL[icell];
avomega[icell][1] /= nCELL[icell];
avomega[icell][2] /= nCELL[icell];
//normalization
if (domain->dimension == 2){
avomega[icell][2] = 0.0;
}
double msq = sqrt(avomega[icell][0]*avomega[icell][0] + avomega[icell][1]*avomega[icell][1] + avomega[icell][2]*avomega[icell][2]);
if(msq > 0.0){
avomega[icell][0] /= msq;avomega[icell][1] /= msq;avomega[icell][2] /= msq;
}
}
for (n = 0; n < nbondlist; n++) {
i1 = bondlist[n][0];
type = bondlist[n][2];
omega[i1][0] = -alpha[type] * avomega[ mol2cell[molecule[i1]] ][0];
omega[i1][1] = -alpha[type] * avomega[ mol2cell[molecule[i1]] ][1];
omega[i1][2] = -alpha[type] * avomega[ mol2cell[molecule[i1]] ][2];
if (screen) fprintf(screen,“omega_cell: %f %f \n”,omega[i1][0],omega[i1][1]); //output for checking the values
}
//if (screen) fprintf(screen,“index of local atom %d \n”,2);
}