Questions and improvements about the pair_ylz.cpp file

Dear, Dr Axel,

I download the LAMMPS source with git, release branch, on ubuntu 20.04, compile lammps by make.

I am one of the participants of src/ASPHERE/pair_ylz.cpp. When submitting to github, you added a piece of code for ‘check that all atoms are ellipsoids’ at lines 263-273. At that time, you said that if you didn’t add check, there would be an error.

But we have tested the examples/ASPHERE examples of flat_membrane and vesicle. There is no error if you do not add this check code.

In addition, we also have a red blood cell calculation example that needs to use pair_ylz, in which atom_style is hybrid, ellipsoid and molecular, then an error will be reported if check is carried out at this time. But once I comment it out, the code runs smoothly.

Therefore, I would like to ask why this code should be added, or can it be commented out or deleted?

I am a new user, so I cannot put my input file of my red blood cell model here. I have sent my input file to you by email.

We are looking forward to hearing from you.

Thank you!

Yes, there would be an error with legitimate data files, where not all atoms are flagged as ellipsoids. In the statement marked below you unconditionally access the “atomvec->bonus” array for all atoms the pair style is applied to, and thus assume that all elements of the ellipsoid array are >=0. Those will be -1 for non-ellipsoid atoms.

  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
  int *ellipsoid = atom->ellipsoid;

[...]

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    itype = type[i];

    iquat = bonus[ellipsoid[i]].quat; // HERE is the problem
    MathExtra::quat_to_mat_trans(iquat, a1);

Here is the equivalent block form pair style gayberne:

  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
  int *ellipsoid = atom->ellipsoid;

[...]

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    itype = type[i];

    if (form[itype][itype] == ELLIPSE_ELLIPSE) {
      iquat = bonus[ellipsoid[i]].quat;
      MathExtra::quat_to_mat_trans(iquat,a1);
      MathExtra::diag_times3(well[itype],a1,temp);
      MathExtra::transpose_times3(a1,temp,b1);
      MathExtra::diag_times3(shape2[itype],a1,temp);
      MathExtra::transpose_times3(a1,temp,g1);
    }

Here the bonus data is only accessed when the pair has been identified as an interaction between ellipsoids.

Because in those inputs all atoms are flagged as ellipsoids and have the corresponding ellipsoid data. If that was not the case, your pair style would cause a segmentation fault.

Unless you modify your pair style to have well defined interactions that include non-ellipsoid atoms, there must be a check. To be compatible with a hybrid pair style a more complex check is needed that uses the neighbor list, which will only include pairs of atoms for which the specific pair style is chosen.

The following commit implements that more elaborate check and will be included in the next feature release of LAMMPS (which we expect to be ready in a few days).