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: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 30 Jun 2017 11:18:06 -0400

On Fri, Jun 30, 2017 at 11:02 AM, Mordechai Kornbluth <mkornbl@...24...> wrote:
> I apologize if this is not the correct place to post this.

it is OK to post your observations here.
an alternative would be to file an issue at
https://github.com/lammps/lammps/issues
since the latter is quite new and we are still refining our guidelines
which should be sent where.
thus for the time being we accept both.

i am copying the author of the code so you have an authoritative
source for the items, that i don't know about.

> 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.

you are correct. i observed an issue, where out-of-bounds array
accesses would happen and suggested two options to resolve this. the
chosen resolution was to require a complete set of entries in the
potential file by adding dummy entries. this should indeed by
documented. i will add this.

i have to defer to answers to the rest of your questions to others.

axel.

> 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.
> [1] http://lammps.sandia.gov/doc/pair_nb3b_harmonic.html
> [2] http://lammps.sandia.gov/threads/msg56857.html
> [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]
> https://github.com/lammps/lammps/blob/master/src/MANYBODY/pair_nb3b_harmonic.cpp
> [5] lines 117-130 here:
> https://github.com/lammps/lammps/blob/aa6624c029fd5eb32aeec35ab41e79574ddf6a29/src/MANYBODY/pair_nb3b_harmonic.cpp
>
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> lammps-users mailing list
> lammps-users@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/lammps-users
>



-- 
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.