LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] pair_style nb3b/harmonic implementation vs documentation
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] pair_style nb3b/harmonic implementation vs documentation

From: Mordechai Kornbluth <mkornbl@...24...>
Date: Fri, 30 Jun 2017 11:02:08 -0400

I apologize if this is not the correct place to post this.

The nb3b/harmonic pair style doesn’t seem to behave as expected or documented with respect to the input-file interpretation. [0] I would greatly appreciate it if someone would enlighten me whether I am misreading the code & documentation, or if this should be submitted as an issue to be fixed.

The documentation [1] and a subsequent mailing-list discussion [2] indicate a few major features of this pair_style syntax:
1. It is OK if not all triplets (i,j,k) are listed in the potential file.
2. The value of K corresponds to the energy of a displacement of the angle of 1 radian. [3]
3. If atoms j and k are the same, the cutoff is used, and the physical parameters (K, theta_0) are totally ignored.
4. If atoms j and k are different, the physical parameters (K, theta_0) are used and the cutoff is totally ignored

However, the source code [4] indicates the following:
1. In setup_params, line 425, an error is raised if the potential file is missing an entry, against #1 above.
    - This was added in the commit of 30 Jul 2016, and seems to make the documentation now out-of-date.
2. The i-j-k loop iterates over all triplets, meaning that the (i,j,k) angle will be counted twice: Once as (i,j,k) and once as (i,k,j). If so, the energy of a displacement of 1 radian would correspond to 2K due to this double-counting.
3. Even if atoms j and k are the same,
    - The K, theta_0 parameters are read (in read_file, lines 386-387)
    - The i,j,k triplet is part of the loop (in compute, lines 103-125, over all atoms)
    - The parameters (K,theta_0) are applied (in threebody, lines 482-483)
    - This seems to be an undocumented effect of the commit on 15 Oct 2015, a mere 3 days after the mailing list discussion [2].
    - Yet the commit message indicates only the changing of > to >=, not the deletion of the ~25 lines of code from pair_nb3b_harmonic.cpp [4].
    - Also, the deleted lines [4] don’t seem to be doing anything; perhaps in an earlier form they ensured the angle was not double-counted (point 2 above).
4. The implementation is in agreement with the documentation. However, perhaps the developers should think about implementing a cutoff for the (i,j,k) triplet; it entails merely adding the following two lines before the call to threebody on line 138, and of course modifying the documentation appropriately:
    - if (rsq1 > params[ijkparam].cutsq) continue;
    - if (rsq2 > params[ijkparam].cutsq) continue;

Thanks in advance!

[0] I am assuming a convention where the input “i j k” identifies an angle where “i” is the central atom, i.e. the j-k-i angle. This is the convention of the nb3b.harmonic file.
[3] Although theta_0 is given in degrees, it is apparent from the source code [3] l.386,440,483 that K is in energy/radian^2.
[5] lines 117-130 here: