LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Accessing number of bonds in python operated LAMMPS or length of internal LAMMPS data array
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Accessing number of bonds in python operated LAMMPS or length of internal LAMMPS data array


From: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 13 Apr 2018 11:04:15 -0400

On Fri, Apr 13, 2018 at 10:42 AM, Robert Meißner
<robert.meissner@...4427...> wrote:
> Thanks for the prompt reply, Axel!
>
> Assuming I want to extract the topology in a gather_atoms() fashion, in
> which variable is this info stored? Atom info is in lmp->atom->… but where
> can I find the bonds? Probably I need furthermore an nlocal counter for the
> number of bonds, since the number of nlocal does not match the number of
> local bonds. Maybe you know a function to which I can refer as a blueprint?

persistent bond information is stored with atoms in the atom class in
lmp->atom->num_bond
lmp->atom->bond_type
lmp->atom->bond_atom

those are all per-atom arrays and may be NULL, if the atom style does
not support bonds.

compute_bond_local() uses that information to compile the local bond lists.
thus my suggestion that the approach requiring the least amount of
changes, would be to add a scalar compute function to compute
bond/local, that triggers the local compute, but also stores the
number of local bonds in a scalar number. that would not require any
changes to the library interface. you just query the compute twice
(and the result of the local compute would be cached, so it should
come at no additional cpu cost).

> Since we're on the subject, how does the command extract_box() command work?
> I did several approaches but was not able to extract the periodicity or
> anything using this command.

are you talking about the C interface or the python wrapper?
the C interface should be obvious from the prototype and a look at library.cpp
for python, you would have to create suitably dimensioned and typed
ctype arrays and pass them as arguments.

in python/lammps.py you have:

   self.lib.lammps_extract_box.argtypes = \
      [c_void_p,POINTER(c_double),POINTER(c_double),
       POINTER(c_double),POINTER(c_double),POINTER(c_double),
       POINTER(c_int),POINTER(c_int)]

so you have pointers to an array of 3 doubles for boxlo, boxhi, then
pointers to 1 double for xy, xy, yx, then 3 ints for periodicity and a
pointer to a single int for box_change.

> Furthermore, thanks for the hints when changing the bond topology. I plan to
> use create_bonds or delete_bonds (or maybe fix bonds/create) directly from a
> LAMMPS input script to build new topology and avoid the problems you’ve
> mentioned. The python function should just identify atom IDs and store them
> in variables for later use.

if you are working on something like this, it is probably advisable to
work off of the unstable branch from github (or bitbucket) and keep
that up-to-date with each patch release.
python support is under active development, and if you stick with an
old version, you may later have a really bad time to port it to the
latest version. doing this incrementally is usually much less effort.

cheers,
    axel.


> Best
> R
>
>
> Am 13.04.2018 um 14:53 schrieb Axel Kohlmeyer <akohlmey@...24...>:
>
> On Fri, Apr 13, 2018 at 7:51 AM, Robert Meißner <robert.meissner@...12...4427...>
> wrote:
>
> Dear LAMMPS Developers,
>
> since I'm facing the same problem right now, has there been any progress
> achieved on this matter? I want to extract the bond topology to
> create/delete some distinct bonds by identifying atoms from their bonded
> neighbors. If I'm running LAMMPS from python (thus in serial) the following
>
>
> FWIW, you can also launch the python interpreter in parallel (via
> pypar or mpi4py).
>
>
> code will give me something which works:
>
> ---
>
> import numpy as np
>
> bonds = np.array([],dtype=int).reshape(-1,3);
>
> def bond_topology(lmpptr):
>  from lammps import lammps
>  lmp = lammps(ptr=lmpptr)
>
>  global bonds
>
>  lmp.command("compute bonds all property/local btype batom1 batom2")
>  lmp.command("run 0")
>  nbonds = lmp.extract_global("nbonds",0)
>  b = lmp.extract_compute("bonds",2,2)
>
>  for i in range(nbonds):
>    if b[i][0]: bonds = np.vstack([bonds,[b[i][0],b[i][1],b[i][2]]])
>  lmp.command("uncompute bonds")
>
>
> this just works by accident and not by construction. it looks to me,
> the straightforward way to address this would be to add a scalar
> context to compute bond/local (and other local computes) which then
> allows to obtain the desired (and currently inaccessible) information
> of the local array lengths.
>
> steve, do you agree? or am i overlooking something?
>
> ---
>
> However, when running python from LAMMPS obviously every MPI proc starts his
> own python interpreter and has only his local bond topology making this code
> fragment pretty useless since I need the full topology. Maybe someone has
> already made a routine like lmp.gather_bonds(...) or something similar to
> get the topology from python?
>
>
> what you see (in the source) is what you get. the philosophy is that
> library.cpp is deliberately incomplete and that only then functions
> are added to library.cpp, when people have a demonstrated need that
> cannot be satisfied differently.
> yours certainly fits the bill, however, i am not aware of any effort
> to make the topology information available.
>
> in general, please be very cautious when making changes to the bond
> topology in LAMMPS. it is quite tricky to do so when running in
> parallel. for performance reasons, many parts of LAMMPS assume that
> the bond topology does not change, so data structures are initialized
> at the beginning of a run. so those need to be updated when changes
> are made, and in particularly updated consistently across all parallel
> ranks. the (rather recently added) ability to access the python
> interface from within a run and thus perform modifications, that were
> previously only possible outside a run, has opened the door to a lot
> of ways to make the internal state of LAMMPS inconsistent. if you are
> using the very latest version of LAMMPS, multiple problematic classes
> of interface calls are now prohibited.
>
> axel.
>
>
>
> Best
>
> Robert
>
>
>
> --
> Prof Dr.-Ing Robert Horst Meißner
>
> Institute of Polymer and Composites
> Denickestr. 15
> D-21073 Hamburg
> Tel.:  +49 (0)40 42878 2580
> E-mail: robert.meissner@...4427...
> Homepage: www.tuhh.de/alt/kvweb
>
>
> ------------------------------------------------------------------------------
> 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.
>
>
>
> ------------------------------------------------------------------------------
> 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.