Unverified Commit 8c52032b authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #1902 from stanmoore1/kk_compute_coord_atom 

Add Kokkos version of compute coord/atom
parents 53ac67f5 0cf56360

doc/src/Commands_compute.rst

+1 −1
Original line number Diff line number Diff line
@@ -42,7 +42,7 @@ KOKKOS, o = USER-OMP, t = OPT. 
   * :doc:`com <compute_com>`

   * :doc:`com/chunk <compute_com_chunk>`

   * :doc:`contact/atom <compute_contact_atom>`

   * :doc:`coord/atom <compute_coord_atom>`

   * :doc:`coord/atom (k) <compute_coord_atom>`

   * :doc:`damage/atom <compute_damage_atom>`

   * :doc:`dihedral <compute_dihedral>`

   * :doc:`dihedral/local <compute_dihedral_local>`

doc/src/compute_coord_atom.rst

+27 −0
Original line number Diff line number Diff line
@@ -3,6 +3,9 @@ 
compute coord/atom command

==========================



compute coord/atom/kk command

===================================



Syntax

""""""
@@ -109,6 +112,30 @@ too frequently. 
   :doc:`special_bonds <special_bonds>` command that includes all pairs in

   the neighbor list.



----------





Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are

functionally the same as the corresponding style without the suffix.

They have been optimized to run faster, depending on your available

hardware, as discussed on the :doc:`Speed packages <Speed_packages>` doc

page.  The accelerated styles take the same arguments and should

produce the same results, except for round-off and precision issues.



These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,

USER-OMP and OPT packages, respectively.  They are only enabled if

LAMMPS was built with those packages.  See the :doc:`Build package <Build_package>` doc page for more info.



You can specify the accelerated styles explicitly in your input script

by including their suffix, or you can use the :doc:`-suffix command-line switch <Run_options>` when you invoke LAMMPS, or you can use the

:doc:`suffix <suffix>` command in your input script.



See the :doc:`Speed packages <Speed_packages>` doc page for more

instructions on how to use the accelerated styles effectively.





----------



**Output info:**



For *cstyle* cutoff, this compute can calculate a per-atom vector or

doc/src/compute_orientorder_atom.rst

+11 −2
Original line number Diff line number Diff line
@@ -19,13 +19,14 @@ Syntax 


  .. parsed-literal::



     keyword = *cutoff* or *nnn* or *degrees* or *components*

     keyword = *cutoff* or *nnn* or *degrees* or *components* or *chunksize*

       *cutoff* value = distance cutoff

       *nnn* value = number of nearest neighbors

       *degrees* values = nlvalues, l1, l2,...

       *wl* value = yes or no

       *wl/hat* value = yes or no

       *components* value = ldegree

       *chunksize* value = number of atoms in each pass 



Examples

""""""""
@@ -107,6 +108,14 @@ in conjunction with :doc:`compute coord_atom <compute_coord_atom>` to 
calculate the ten Wolde's criterion to identify crystal-like

particles, as discussed in :ref:`ten Wolde <tenWolde2>`.



The optional keyword *chunksize* is only applicable when using the

the KOKKOS package and is ignored otherwise. This keyword controls

the number of atoms in each pass used to compute the bond-orientational

order parameters and is used to avoid running out of memory. For example

if there are 4000 atoms in the simulation and the *chunksize*

is set to 2000, the parameter calculation will be broken up

into two passes.



The value of :math:`Q_l` is set to zero for atoms not in the

specified compute group, as well as for atoms that have less than

*nnn* neighbors within the distance cutoff, unless *nnn* is NULL.
@@ -192,7 +201,7 @@ Default 


The option defaults are *cutoff* = pair style cutoff, *nnn* = 12,

*degrees* = 5 4 6 8 10 12 i.e. :math:`Q_4`, :math:`Q_6`, :math:`Q_8`, :math:`Q_{10}`, and :math:`Q_{12}`,

*wl* = no, *wl/hat* = no, and *components* off

*wl* = no, *wl/hat* = no, *components* off, and *chunksize* = 2000



----------

src/KOKKOS/Install.sh

+2 −0
Original line number Diff line number Diff line
@@ -85,6 +85,8 @@ action comm_kokkos.cpp 
action comm_kokkos.h

action comm_tiled_kokkos.cpp

action comm_tiled_kokkos.h

action compute_coord_atom_kokkos.cpp

action compute_coord_atom_kokkos.h

action compute_orientorder_atom_kokkos.cpp

action compute_orientorder_atom_kokkos.h

action compute_temp_kokkos.cpp

src/KOKKOS/compute_coord_atom_kokkos.cpp

0 → 100644
+249 −0
Original line number Diff line number Diff line 
/* ----------------------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov



   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.



   See the README file in the top-level LAMMPS directory.

------------------------------------------------------------------------- */



#include "compute_coord_atom_kokkos.h"

#include <cmath>

#include <cstring>

#include "compute_orientorder_atom_kokkos.h"

#include "atom_kokkos.h"

#include "update.h"

#include "modify.h"

#include "neighbor_kokkos.h"

#include "neigh_list.h"

#include "neigh_request.h"

#include "force.h"

#include "pair.h"

#include "comm.h"

#include "group.h"

#include "memory_kokkos.h"

#include "error.h"

#include "atom_masks.h"



using namespace LAMMPS_NS;



#define INVOKED_PERATOM 8



/* ---------------------------------------------------------------------- */



template<class DeviceType>

ComputeCoordAtomKokkos<DeviceType>::ComputeCoordAtomKokkos(LAMMPS *lmp, int narg, char **arg) :

  ComputeCoordAtom(lmp, narg, arg)

{

  kokkosable = 1;

  atomKK = (AtomKokkos *) atom;

  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;

  datamask_read = EMPTY_MASK;

  datamask_modify = EMPTY_MASK;



  d_typelo = typename AT::t_int_1d("coord/atom:typelo",ncol);

  d_typehi = typename AT::t_int_1d("coord/atom:typehi",ncol);



  auto h_typelo = Kokkos::create_mirror_view(d_typelo);

  auto h_typehi = Kokkos::create_mirror_view(d_typehi);



  for (int i = 0; i < ncol; i++) {

    h_typelo(i) = typelo[i];

    h_typehi(i) = typehi[i];

  }



  Kokkos::deep_copy(d_typelo,h_typelo);

  Kokkos::deep_copy(d_typehi,h_typehi);

}



/* ---------------------------------------------------------------------- */



template<class DeviceType>

ComputeCoordAtomKokkos<DeviceType>::~ComputeCoordAtomKokkos<DeviceType>()

{

  if (copymode) return;



  memoryKK->destroy_kokkos(k_cvec,cvec);

  memoryKK->destroy_kokkos(k_carray,carray);

}



/* ---------------------------------------------------------------------- */



template<class DeviceType>

void ComputeCoordAtomKokkos<DeviceType>::init()

{

  ComputeCoordAtom::init();



  // need an occasional full neighbor list



  // irequest = neigh request made by parent class



  int irequest = neighbor->nrequest - 1;



  neighbor->requests[irequest]->

    kokkos_host = std::is_same<DeviceType,LMPHostType>::value &&

    !std::is_same<DeviceType,LMPDeviceType>::value;

  neighbor->requests[irequest]->

    kokkos_device = std::is_same<DeviceType,LMPDeviceType>::value;

}



/* ---------------------------------------------------------------------- */



template<class DeviceType>

void ComputeCoordAtomKokkos<DeviceType>::compute_peratom()

{

  invoked_peratom = update->ntimestep;



  // grow coordination array if necessary



  if (atom->nmax > nmax) {

    if (ncol == 1) {

      memoryKK->destroy_kokkos(k_cvec,cvec);

      nmax = atom->nmax;

      memoryKK->create_kokkos(k_cvec,cvec,nmax,"coord/atom:cvec");

      vector_atom = cvec;

      d_cvec = k_cvec.template view<DeviceType>();

    } else {

      memoryKK->destroy_kokkos(k_carray,carray);

      nmax = atom->nmax;

      memoryKK->create_kokkos(k_carray,carray,nmax,ncol,"coord/atom:carray");

      array_atom = carray;

      d_carray = k_carray.template view<DeviceType>();

    }

  }



  if (cstyle == ORIENT) {

    if (!(c_orientorder->invoked_flag & INVOKED_PERATOM)) {

      c_orientorder->compute_peratom();

      c_orientorder->invoked_flag |= INVOKED_PERATOM;

    }

    nqlist = c_orientorder->nqlist;

    normv = c_orientorder->array_atom;

    comm->forward_comm_compute(this);



    if (!c_orientorder->kokkosable)

      error->all(FLERR,"Must use compute orientorder/atom/kk with compute coord/atom/kk");



    if (c_orientorder->execution_space == Host) {

      ComputeOrientOrderAtomKokkos<LMPHostType>* c_orientorder_kk;

      c_orientorder_kk = (ComputeOrientOrderAtomKokkos<LMPHostType>*) c_orientorder;

      c_orientorder_kk->k_qnarray.modify<LMPHostType>();

      c_orientorder_kk->k_qnarray.sync<DeviceType>();

      d_normv = c_orientorder_kk->k_qnarray.view<DeviceType>();

    } else {

      ComputeOrientOrderAtomKokkos<LMPDeviceType>* c_orientorder_kk;

      c_orientorder_kk = (ComputeOrientOrderAtomKokkos<LMPDeviceType>*) c_orientorder;

      c_orientorder_kk->k_qnarray.modify<LMPHostType>();

      c_orientorder_kk->k_qnarray.sync<DeviceType>();

      d_normv = c_orientorder_kk->k_qnarray.view<DeviceType>();

    }

  }



  // invoke full neighbor list (will copy or build if necessary)



  neighbor->build_one(list);



  inum = list->inum;

  NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);

  d_numneigh = k_list->d_numneigh;

  d_neighbors = k_list->d_neighbors;

  d_ilist = k_list->d_ilist;



  // compute coordination number(s) for each atom in group

  // use full neighbor list to count atoms less than cutoff



  atomKK->sync(execution_space,X_MASK|TYPE_MASK|MASK_MASK);

  x = atomKK->k_x.view<DeviceType>();

  type = atomKK->k_type.view<DeviceType>();

  mask = atomKK->k_mask.view<DeviceType>();



  copymode = 1;

  if (cstyle == CUTOFF) {

    if (ncol == 1) {

      typename Kokkos::RangePolicy<DeviceType, TagComputeCoordAtom<CUTOFF,1> > policy(0,inum);

      Kokkos::parallel_for("ComputeCoordAtom",policy,*this);

    } else {

      typename Kokkos::RangePolicy<DeviceType, TagComputeCoordAtom<CUTOFF,0> > policy(0,inum);

      Kokkos::parallel_for("ComputeCoordAtom",policy,*this);

    }

  } else if (cstyle == ORIENT) {

    typename Kokkos::RangePolicy<DeviceType, TagComputeCoordAtom<ORIENT,1> > policy(0,inum);

    Kokkos::parallel_for("ComputeCoordAtom",policy,*this);

  }

  copymode = 0;



  if (ncol == 1 || cstyle == ORIENT) {

    k_cvec.modify<DeviceType>();

    k_cvec.sync<LMPHostType>();

  } else {

    k_carray.modify<DeviceType>();

    k_carray.sync<LMPHostType>();

  }



}



template<class DeviceType>

template<int CSTYLE, int NCOL>

KOKKOS_INLINE_FUNCTION

void ComputeCoordAtomKokkos<DeviceType>::operator()(TagComputeCoordAtom<CSTYLE,NCOL>, const int &ii) const

{

  const int i = d_ilist[ii];

  if (NCOL == 1 || CSTYLE == ORIENT)

    d_cvec(i) = 0.0;

  else

    for (int m = 0; m < ncol; m++) d_carray(i,m) = 0.0;

  if (mask[i] & groupbit) {

    const X_FLOAT xtmp = x(i,0);

    const X_FLOAT ytmp = x(i,1);

    const X_FLOAT ztmp = x(i,2);

    const int jnum = d_numneigh[i];



    int n = 0;

    for (int jj = 0; jj < jnum; jj++) {

      int j = d_neighbors(i,jj);

      j &= NEIGHMASK;



      if (NCOL == 1)

        if (!(mask[j] & jgroupbit)) continue;



      const int jtype = type[j];

      const F_FLOAT delx = x(j,0) - xtmp;

      const F_FLOAT dely = x(j,1) - ytmp;

      const F_FLOAT delz = x(j,2) - ztmp;

      const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;

      if (rsq < cutsq) {

        if (CSTYLE == CUTOFF) {

          if (NCOL == 1) {

            if (jtype >= d_typelo[0] && jtype <= d_typehi[0])

              n++;

          } else {

              for (int m = 0; m < ncol; m++)

                if (jtype >= d_typelo[m] && jtype <= d_typehi[m])

                  d_carray(i,m) += 1.0;

          }

        } else if (CSTYLE == ORIENT) {

            double dot_product = 0.0;

            for (int m=0; m < 2*(2*l+1); m++) {

              dot_product += d_normv(i,nqlist+m)*d_normv(j,nqlist+m);

            }

            if (dot_product > threshold) n++;

        }

      }

    }



    if (NCOL == 1 || CSTYLE == ORIENT)

      d_cvec[i] = n;

  }



}



namespace LAMMPS_NS {

template class ComputeCoordAtomKokkos<LMPDeviceType>;

#ifdef KOKKOS_ENABLE_CUDA

template class ComputeCoordAtomKokkos<LMPHostType>;

#endif

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE