LAMMPS library: Gather specific atoms and insert small molecules

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,
dear Steve and Axel,

I'm following up on this old post on the list.
Trying to get all coordinates for a given
type, e.g. 2, I put this in my LAMMPS script:

group MYGRP type 2
variable MYCOORDS atom x
variable MYCOUNT equal count(MYGRP)

Using the functions from library.cpp by

double* pnatoms = (double*) lammps_extract_variable(lmp, "MYCOUNT", NULL);
                        int natoms = (int)(*pnatoms);
                                           double* coords = (double*)
lammps_extract_variable(lmp, "MYCOORDS", "MYGRP");

However, coords contains only zeros, e.g. 0.0, 0.0, 0.0.
Apparently I'm doing something wrong, but I'm not sure what.

Maybe somebody can point me in the right direction.

Best regards,
Stephan

2

I'm apologizing for the multiple posts,
I had a problem with my mail client.

Sorry folks.

Hope somebody is still interested in my question.

Best regards,
Stephan

3

Dear LAMMPS users,
dear Steve and Axel,

I'm following up on this old post on the list.
Trying to get all coordinates for a given
type, e.g. 2, I put this in my LAMMPS script:

group MYGRP type 2
variable MYCOORDS atom x
variable MYCOUNT equal count(MYGRP)

Using the functions from library.cpp by

double* pnatoms = (double*) lammps_extract_variable(lmp, "MYCOUNT", NULL);

                      int natoms = (int)(*pnatoms);

          double* coords = (double*)
lammps_extract_variable(lmp, "MYCOORDS", "MYGRP");

However, coords contains only zeros, e.g. 0.0, 0.0, 0.0.
Apparently I'm doing something wrong, but I'm not sure what.

it ​works for me. here is a complete(!) test example based on the "crack"
example.

axel.

​#include <stdio.h>
#include <inttypes.h>
#include "library.h"

int main (int argc, char **argv)

{
    void *lmp;
    int i,gatoms;
    int64_t natoms;

    lammps_open_no_mpi(0,NULL,&lmp);
    lammps_file(lmp,"in.crack");
    lammps_command(lmp,"group MYGRP type 2");
    lammps_command(lmp,"variable MYCOORDS atom x");
    lammps_command(lmp,"variable MYCOUNT equal count(MYGRP)");

    natoms = *(int64_t *)lammps_extract_global(lmp,"natoms");
    gatoms = (int) *(double *)lammps_extract_variable(lmp, "MYCOUNT", NULL);
    printf("total atoms = %ld atoms in group = %d\n",natoms,gatoms);
    double *coords = (double
*)lammps_extract_variable(lmp,"MYCOORDS","MYGRP");
    for (i=0; i < natoms; ++i) {
        printf("x[%03d]=%10.4f\n",i,coords[i]);
    }
    lammps_free((void *)coords);
    lammps_close(lmp);
    return 0;
}
​

4

Dear Axel,

thanks for the example.
I'm using the lammps stable branch.

When I run your code with my LAMMPS
input script, again, I'm getting only zeros.

I suspect maybe that my input script is wrong,
thus I'm attaching it here:

units real
neigh_modify delay 2 every 1

atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic

pair_style lj/charmm/coul/long 8 10
pair_modify mix arithmetic
kspace_style pppm 1e-4

read_data protein.data

special_bonds charmm
fix 1 all nve
velocity all create 0.0 1504301625 dist uniform

thermo 1
thermo_style multi
timestep 0.005

dump 1 all atom 10 protein.dump
dump_modify 1 image yes scale yes
group NTGRP type 2
variable NT atom x
variable NTCOUNT equal count(NTGRP)

Is there anything fishy here?
Thanks in advance.

Best regards,
Stephan

5

Dear Axel,

thanks for the example.
I'm using the lammps stable branch.

When I run your code with my LAMMPS
input script, again, I'm getting only zeros.

​are you sure you don't miss the real data. how many atoms are in the
selected group?
perhaps you should prefix the test print​, so it prints out only the
non-zero numbers, e.g. like this:

    for (i=0; i < natoms; ++i) {
        if (coords[i] != 0.0)
            printf("x[%03d]=%10.4f\n",i,coords[i]);
    }

i've taken your input and then used the data file from the LAMMPS peptide
example, and i get - with the modification from above - the following
output.

LAMMPS (17 Aug 2017)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread.
(../comm.cpp:90)
  using 1 OpenMP thread(s) per MPI task
Reading data file ...
  orthogonal box = (36.8402 41.0137 29.7681) to (64.2116 68.3851 57.1395)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  2004 atoms
  reading velocities ...
  2004 velocities
  scanning bonds ...
  3 = max bonds/atom
  scanning angles ...
  6 = max angles/atom
  scanning dihedrals ...
  14 = max dihedrals/atom
  scanning impropers ...
  1 = max impropers/atom
  reading bonds ...
  1365 bonds
  reading angles ...
  786 angles
  reading dihedrals ...
  207 dihedrals
  reading impropers ...
  12 impropers
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj: 0 0 0
  special bond factors coul: 0 0 0
  4 = max # of 1-2 neighbors
  7 = max # of 1-3 neighbors
  14 = max # of 1-4 neighbors
  18 = max # of special neighbors
3 atoms in group NTGRP
3 atoms in group MYGRP
total atoms = 2004 atoms in group = 3
x[001]= 45.1039
x[068]= 47.7078
x[079]= 49.3779
Total wall time: 0:00:00

this looks perfectly alright to me. 3 atoms in the group -> 3 non-zero
x-coordinates in the variable.

axel.

6

Dear Axel,

thanks for your excellent help!

I used gatoms in the for loop, not natoms.
Now the question is, how could I get only
the x coord of the three atoms in the group
without checking for inequality with 0.0?

Best regards,
Stephan

7

Dear Axel,

thanks for your excellent help!

I used gatoms in the for loop, not natoms.
Now the question is, how could I get only
the x coord of the three atoms in the group
without checking for inequality with 0.0?

​you could pull in another atom style variable that contains the output of
​gmask(GROUP), which is 1 for group members and 0 for other atoms, and then
use this variable as selector.

axel.

8

Dear Axel,

thanks, this sounds good.

What I'm not understanding so far is,
that I thought by using something like

double *coords = (double*)lammps_extract_variable(lmp,"MYCOORDS","MYGRP");

I would get only the coordinates for the atoms in MYGRP
and thus would not have to use gmask as selector?

A last question (hopefully) concerning this:
If I run that code now with MPI (more than one proc),
I get some garbage output, like:
x[1005]=7199496020021154586587594...

Is this to be expected without using the selector?

Thanks in advance.

Best regards,
Stephan

9

Dear Axel,

thanks, this sounds good.

What I'm not understanding so far is,
that I thought by using something like

double *coords = (double*)lammps_extract_variable(lmp,"MYCOORDS","MYGRP");

I would get only the coordinates for the atoms in MYGRP
and thus would not have to use gmask as selector?

​no. atom style variables are always returned as​ an array of the size of
the number of local atoms.

A last question (hopefully) concerning this:
If I run that code now with MPI (more than one proc),
I get some garbage output, like:
x[1005]=7199496020021154586587594...

Is this to be expected without using the selector?

​no. this is because the demo code i wrote only works for 1 processor. for
1 processor nlocal == natoms.
atom style variables are stored distributed across the domains, so the
length of the array is nlocal and not natoms for more than one MPI rank.

sorry about that. so the more correct library interface test code looks
like this:

#include <mpi.h>
#include <stdio.h>
#include <inttypes.h>
#include "library.h"

int main (int argc, char **argv)

{
    void *lmp;
    int i,gatoms;
    int nlocal;

    MPI_Init(&argc,&argv);
    lammps_open(0,NULL,MPI_COMM_WORLD,&lmp);
    lammps_file(lmp,"in.protein");
    lammps_command(lmp,"group MYGRP type 2");
    lammps_command(lmp,"variable MYCOORDS atom x");
    lammps_command(lmp,"variable MYIDS atom id");
    lammps_command(lmp,"variable MYCOUNT equal count(MYGRP)");

    nlocal = *(int *)lammps_extract_global(lmp,"nlocal");
    gatoms = (int) *(double *)lammps_extract_variable(lmp,"MYCOUNT",NULL);
    printf("local atoms = %d atoms in group = %d\n",nlocal,gatoms);
    double *coords = (double
*)lammps_extract_variable(lmp,"MYCOORDS","MYGRP");
    double *ids = (double *)lammps_extract_variable(lmp,"MYIDS","MYGRP");
    for (i=0; i < nlocal; ++i) {
        if (ids[i] != 0.0)
            printf("x[%03d]=%10.4f\n",(int)ids[i],coords[i]);
    }
    lammps_free((void *)coords);
    lammps_close(lmp);
    MPI_Finalize();
    return 0;
}

​axel.​

10

Dear Axel,

okay this makes sense. Thank you!
(And thus works now)

    Dear Axel,

    thanks, this sounds good.

    What I'm not understanding so far is,
    that I thought by using something like

    double *coords =
    (double*)lammps_extract_variable(lmp,"MYCOORDS","MYGRP");

    I would get only the coordinates for the atoms in MYGRP
    and thus would not have to use gmask as selector?

​no. atom style variables are always returned as​ an array of the size
of the number of local atoms.

Ok.

    A last question (hopefully) concerning this:
    If I run that code now with MPI (more than one proc),
    I get some garbage output, like:
    x[1005]=7199496020021154586587594...

    Is this to be expected without using the selector?

​no. this is because the demo code i wrote only works for 1 processor.
for 1 processor nlocal == natoms.
atom style variables are stored distributed across the domains, so the
length of the array is nlocal and not natoms for more than one MPI rank.

Great, thank you.

I've to be sorry because I probably could have found it
in Section 6.19 of the docu. Is the documentation on
http://lammps.sandia.gov/doc/Section_howto.html#library-interface-to-lammps
the most recent and is there more documentation I miss?

sorry about that. so the more correct library interface test code looks
like this:

#include <mpi.h>
#include <stdio.h>
#include <inttypes.h>
#include "library.h"

int main (int argc, char **argv)

{
void *lmp;
int i,gatoms;
int nlocal;

MPI\_Init\(&amp;argc,&amp;argv\);
lammps\_open\(0,NULL,MPI\_COMM\_WORLD,&amp;lmp\);
lammps\_file\(lmp,&quot;in\.protein&quot;\);
lammps\_command\(lmp,&quot;group MYGRP type 2&quot;\);
lammps\_command\(lmp,&quot;variable MYCOORDS atom x&quot;\);
lammps\_command\(lmp,&quot;variable MYIDS atom id&quot;\);
lammps\_command\(lmp,&quot;variable MYCOUNT equal count\(MYGRP\)&quot;\);

nlocal = \*\(int \*\)lammps\_extract\_global\(lmp,&quot;nlocal&quot;\);
gatoms = \(int\) \*\(double \*\)lammps\_extract\_variable\(lmp,&quot;MYCOUNT&quot;,NULL\);
printf\(&quot;local atoms = %d  atoms in group = %d\\n&quot;,nlocal,gatoms\);
double \*coords = \(double

*)lammps_extract_variable(lmp,"MYCOORDS","MYGRP");
double *ids = (double *)lammps_extract_variable(lmp,"MYIDS","MYGRP");
for (i=0; i < nlocal; ++i) {
if (ids[i] != 0.0)
printf("x[%03d]=%10.4f\n",(int)ids[i],coords[i]);
}
lammps_free((void *)coords);
lammps_close(lmp);
MPI_Finalize();
return 0;
}

​axel.​

Best,
Stephan

11

Dear Axel,

okay this makes sense. Thank you!
(And thus works now)

​[...]
​

I've to be sorry because I probably could have found it
in Section 6.19 of the docu. Is the documentation on
http://lammps.sandia.gov/doc/Section_howto.html#library-
interface-to-lammps
the most recent and is there more documentation I miss?

​the online documentation should always be the most up to date.
beyond that, the most complete docs are always the source code. library.cpp
has plenty of comments about the implementation.

​axel.​

12

Thanks!
So now I'm getting the local atoms in my group on each proc.
However, I would need all atoms in the group on each proc.

Would you suggest to just do a plain MPI_Allreduce on
the double* coordsX? Or is there a LAMMPS library
function I've overseen?

Thanks in advance.

Best regards,
Stephan