Question about omega

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);
}

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 <mu> for each atom. My idea is that I use atom->omega to store the input value <mu> 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.

neither of the two atom styles (bond or dipole) maintains the omega
property. for that you would need to add "sphere" as well.

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.

unlikely. please also note that the implementation of atom styles has
changed recently in the development branch (master) and thus you are
advised to use at least the latest unstable branch to base your work
on. or else you will likely be cut off from any future new features,
bug fixes and improvements.

Then I try to modify the bond_harmonic.h/cpp (to bond_cell.h/cpp) and calculate the average <mu> by traversing the bond list and then assign the <mu> to each atom on the same molecule via storing them in atom->omega( the key code is given below).

that sounds like something that could be done using the chunk feature,
compute fragment/atom will allow to detect molecules by bond topology,
the you can use compute chunk/atom and compute ave/chunk to compute
the average (possibly via an atom style variable, if you need to do
some extra work. you can use fix property/atom and compute
chunk/spread, if you want to redistribute the per-chunk value to the
constituent atoms.

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.

this seems all to be due to the omission in the definition of the atom
style. also there appears to be much less work needed to do explicit
programming by utilizing the existing facilities in LAMMPS (which also
are likely to be more correct, as what you are doing is somewhat
complex and thus has a potential for mistakes).

axel.

Hi Axel,

Thanks for your response. I try to explain what I want to do more clearly below. The LAMMPS version is Mar2020.

The function I want to achieve is that the orientation of each atom on one molecule is updated in dependence of the average orientation of all atoms on the molecule.

I modify all *.cpp/hs with aims to add a new integration command, fix nve/cell to update the orientation of dipole, not with aims to calculate the .
That is to say, I’d like to solve the equation d(mu)_i/dt = -k(), here mu_i is for each atom of one molecule, is the average mu for all atoms on the molecule. I just use the omega quality to restore the for each atom. I use the bond_cell.h/cpp (because I want to traverse the bond list) to calculate the . Then I add a Euler integration method in the fix_nve_cell.h/cpp to solve equation.

I have add the omega property via modifying the atom_vec_dipole.h/cpp, like atom_vec_sphere.h/cpp(I have compared the difference between atom_vec_dipole and actom_vec_sphere, and copy all sub-functions related to omega in the atom_vec_dipole).

The question is that, when, I use **omega = atom->omega to access the omega qualities, it is zero, not the value I assigned. For example,
In the bond_cell.cpp file, the code below

atom_vec_dipole.h (2.62 KB)

atom_vec_dipole.cpp (28.4 KB)

bond_cell.cpp (8.97 KB)

fix_nve_cell.cpp (4.87 KB)

fix_nve_cell.h (1.61 KB)

bond_cell.h (1.64 KB)

Sorry, that is not something that I have time for to discuss in more
detail and assist you with.

as the saying goes, you have already moved too far off the reservation
and made changes that will not be tolerated to be merged into the main
branch of LAMMPS. also, as I already mentioned, we have recently
refactored the way custom atom styles are defined and managed, so any
effort to sort out details will be a waste of my time.

i have already given you an outline of how I would approach something
like this using existing functionality. this is not the only way, but
in general, if you need (internal) per atom storage, you can use fix
STORE, if this has to be visible externally (and may need to be
entered via a section in the data file, or need to be available in
other fixes or computes or variables) you can use fix property/atom.
those are clean, tested, and are used throughout LAMMPS in multiple
places.

what you do is "a hack(tm)" by comparison and has mutliple bad
choices. that doesn't mean it won't work. it just means, that if you
decide to go that way (and LAMMPS is open source, so have every right
to do it), then you are - at least as far as I am concerned - on your
own. I prefer to spend my time on code that has a chance to be
integrated into LAMMPS so that the entire LAMMPS community can
benefit, not just you. This perspective is not given here. As far as
you have the right to choose how and what you work on, so have I.

Sorry,
     Axel.