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

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

From: Julien Guénolé <julien.guenole@...4839...>
Date: Fri, 28 Jul 2017 11:40:05 +0200

On 30/06/17 17:02, Mordechai Kornbluth wrote:
> 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).

I started to work with this potential (in an hybrid scheme) and I
confirm that the 3b interaction is effectively taken into account even
if atom j and k are the same.
I'm however not sure about many other aspect, like the units or how the
cutoff is considered... I need to look more carefully at the code.

Do you have any news on your side on the discrepancies between the doc
and the code? (except the comments from Axel on the 30.06)

> 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.
> [1]
> [2]
> [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.
> [4]
> [5] lines 117-130 here:
Dr. Julien GUÉNOLÉ
Postdoctoral researcher