ALJ potential: Inconsistencies in documentation/code/paper implementations

[navidpanchi]

Hi,

I was reviewing the documentation for the ALJ potential and noticed that the formula there doesn’t seem to match the one in the original paper or the implementation.

The paper describes $\epsilon_s$ (the shape-modified epsilon, referred to as $\epsilon_c$ in the documentation) multiplied with the central interaction term, but the documentation shows it being multiplied with the closest point interaction. From the implementation, it looks like the paper’s version is actually followed.

Also, there may be a bug in the implementation. The parameter $\alpha$ affects $\epsilon_s$. The shape encoding appears to be active only when $\alpha$ is 0 or 2. In other cases, it seems ignored, and the potential reduces to two (isotropic) LJ terms applied to center-center and closest point distances (with or without averaging). Please refer to the code below, from EvaluatePairALJ.h, lines 609–630:

          Scalar k1 = sqrt(dot(a, a));
          Scalar k2 = sqrt(dot(dr + b, dr + b));
          Scalar rho = sigma12
                       / (r - (k1 - Scalar(0.5) * _params.sigma_i)
                          - (k2 - Scalar(0.5) * _params.sigma_j));
          Scalar invr_rsq = rho * rho;
          Scalar invr_6 = invr_rsq * invr_rsq * invr_rsq;
          Scalar numer
              = (invr_6 * invr_6 - invr_6) - (_params.alpha % 2 == 0 ? SHIFT_RHO_DIFF : 0);

          invr_rsq = sigma12_sq / rsq;
          invr_6 = invr_rsq * invr_rsq * invr_rsq;
          Scalar invr_12 = invr_6 * invr_6;
          Scalar denom = invr_12 - invr_6 - (_params.alpha % 2 == 0 ? SHIFT_RHO_DIFF : 0);
          Scalar scale_factor = denom != 0 ? (numer / denom) : 1;

          Scalar four_epsilon = Scalar(4.0) * _params.epsilon;
          Scalar four_scaled_epsilon = four_epsilon;
          if (_params.alpha % 2 == 0)
              {
              four_scaled_epsilon = four_epsilon * scale_factor;
              }

After this point the scale factor is never used again. If this is intended behavior, it's not very clear from the documentation.

Thanks!
Navid Panchi

[joaander]

I was reviewing the documentation for the ALJ potential and noticed that the formula there doesn’t seem to match the one in the original paper or the implementation.

The original author of ALJ graduated several years ago. He may see this and reply, but it is not likely. I have limited knowledge of the model.

The paper describes ϵ s (the shape-modified epsilon, referred to as ϵ c in the documentation) multiplied with the central interaction term, but the documentation shows it being multiplied with the closest point interaction. From the implementation, it looks like the paper’s version is actually followed.

Feel free to submit a pull request that corrects the documentation.

Also, there may be a bug in the implementation. The parameter α affects ϵ s . The shape encoding appears to be active only when α is 0 or 2. In other cases, it seems ignored, and the potential reduces to two (isotropic) LJ terms applied to center-center and closest point distances (with or without averaging). Please refer to the code below, from EvaluatePairALJ.h, lines 609–630:

After this point the scale factor is never used again. If this is intended behavior, it's not very clear from the documentation.

As you said above, the scale factor only applies to the center-center term as stated in the paper (and the documentation is incorrect). In that case, you wouldn't expect the scale factor to continue to be used after the center-center term is computed. I do not know why the scaling is only applied when alpha is 0 or 2.

Most users employ ALJ with alpha=0 to resolve the volume exclusion effects. There are a very small number of cases (specifically polymers adhered to the surface of a shape as the original authors envisioned) where the attractive terms make sense. In the general case of one polyhedron interacting with another polyhedron, the attractive terms in ALJ don't model any realistic physical behavior. It is based only on the closest distance, so there is no energetic preference for a face to face contact over a vertex to vertex contact.

[navidpanchi]

Thanks! I will create a pull request for the documentation update.