support for tiled decompositions in PPPM

/* -*- c++ -*- ----------------------------------------------------------

   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.

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

#ifndef LMP_GRIDCOMM2_H

#define LMP_GRIDCOMM2_H

#include "pointers.h"

namespace LAMMPS_NS {

class GridComm2 : protected Pointers {

 public:

  GridComm2(class LAMMPS *, MPI_Comm, int, int, int,

	    int, int, int, int, int, int,

	    int, int, int, int, int, int);

  GridComm2(class LAMMPS *, MPI_Comm, int, int, int,

	    int, int, int, int, int, int,

	    int, int, int, int, int, int,

	    int, int, int, int, int, int);

  ~GridComm2();

  void setup(int &, int &);

  int ghost_adjacent();

  void forward_comm(class KSpace *, int, int, int, void *, void *, MPI_Datatype);

  void reverse_comm(class KSpace *, int, int, int, void *, void *, MPI_Datatype);

 private:

  int me,nprocs;

  int layout;                 // REGULAR or TILED

  MPI_Comm gridcomm;

  // inputs from caller via constructor

  int nx,ny,nz;               // size of global grid in all 3 dims

  int inxlo,inxhi;            // inclusive extent of my grid chunk

  int inylo,inyhi;            //   0 <= in <= N-1

  int inzlo,inzhi;   

  int outxlo,outxhi;          // inclusive extent of my grid chunk plus

  int outylo,outyhi;          //   ghost cells in all 6 directions

  int outzlo,outzhi;          // lo indices can be < 0, hi indices can be >= N

  int outxlo_max,outxhi_max;  // ??

  int outylo_max,outyhi_max;

  int outzlo_max,outzhi_max;

  // -------------------------------------------

  // internal variables for REGULAR layout

  // -------------------------------------------

  int procxlo,procxhi;     // 6 neighbor procs that adjoin me

  int procylo,procyhi;     //   not used for comm_style = tiled

  int proczlo,proczhi;

  int ghostxlo,ghostxhi;   // # of my owned grid planes needed

  int ghostylo,ghostyhi;   //   by neighobr procs in each dir as their ghost planes

  int ghostzlo,ghostzhi;

  // swap = exchange of owned and ghost grid cells between 2 procs, including self

  struct Swap {

    int sendproc;       // proc to send to for forward comm

    int recvproc;       // proc to recv from for forward comm

    int npack;          // # of datums to pack

    int nunpack;        // # of datums to unpack

    int *packlist;      // 3d array offsets to pack

    int *unpacklist;    // 3d array offsets to unpack

  };

  int nswap,maxswap;

  Swap *swap;

  // -------------------------------------------

  // internal variables for TILED layout

  // -------------------------------------------

  int *overlap_procs;

  MPI_Request *requests;

  // RCB tree of cut info

  // each proc contributes one value, except proc 0

  struct RCBinfo {

    int dim;        // 0,1,2 = which dim the cut is in

    int cut;        // grid index of lowest cell in upper half of cut

  };

  RCBinfo *rcbinfo;

  // overlap = a proc whose owned cells overlap with my extended ghost box

  // includes overlaps across periodic boundaries, can also be self

  struct Overlap {

    int proc;            // proc whose owned cells overlap my ghost cells

    int box[6];          // box that overlaps otherproc's owned cells

                         // this box is wholly contained within global grid

    int pbc[3];          // PBC offsets to convert box to a portion of my ghost box

                         // my ghost box may extend beyond global grid

  };

  int noverlap,maxoverlap;

  Overlap *overlap;

  // request = sent to each proc whose owned cells overlap my ghost cells

  struct Request {

    int sender;          // sending proc

    int index;           // index of overlap on sender

    int box[6];          // box that overlaps receiver's owned cells

                         // wholly contained within global grid

  };

  Request *srequest,*rrequest;

  // response = reply from each proc whose owned cells overlap my ghost cells

  struct Response {

    int index;           // index of my overlap for the initial request

    int box[6];          // box that overlaps responder's owned cells

                         // wholly contained within global grid

                         // has to unwrapped by PBC to map to my ghost cells

  };

  Response *sresponse,*rresponse;

  // send = proc to send a subset of my owned cells to, for forward comm

  // for reverse comm, proc I receive ghost overlaps with my owned cells from

  // offset used in reverse comm to recv a message in middle of a large buffer

  struct Send {

    int proc;

    int npack;

    int *packlist;

    int offset;

  };

  // recv = proc to recv a subset of my ghost cells from, for forward comm

  // for reverse comm, proc I send a subset of my ghost cells to

  // offset used in forward comm to recv a message in middle of a large buffer

  struct Recv {

    int proc;

    int nunpack;

    int *unpacklist;

    int offset;

  };

  int adjacent;      // 0 on a proc who receives ghosts from a non-neighbor proc

  // copy = subset of my owned cells to copy into subset of my ghost cells

  // that describes forward comm, for reverse comm it is the opposite

  struct Copy {

    int npack;

    int nunpack;

    int *packlist;

    int *unpacklist;

  };

  int nsend,nrecv,ncopy;

  Send *send;

  Recv *recv;

  Copy *copy;

  // -------------------------------------------

  // internal methods

  // -------------------------------------------

  void setup_regular(int &, int &);

  void setup_tiled(int &, int &);

  void ghost_box_drop(int *, int *);

  void box_drop_grid(int *, int, int, int &, int *);

  int ghost_adjacent_regular();

  int ghost_adjacent_tiled();

  void forward_comm_regular(class KSpace *, int, int, int,

			    void *, void *, MPI_Datatype);

  void forward_comm_tiled(class KSpace *, int, int, int,

			  void *, void *, MPI_Datatype);

  void reverse_comm_regular(class KSpace *, int, int, int,

			    void *, void *, MPI_Datatype);

  void reverse_comm_tiled(class KSpace *, int, int, int,

			  void *, void *, MPI_Datatype);

  void grow_swap();

  void grow_overlap();

  int indices(int *&, int, int, int, int, int, int);

};

}

#endif

/* -*- c++ -*- ----------------------------------------------------------

   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.

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

#ifdef KSPACE_CLASS

KSpaceStyle(pppm2,PPPM2)

#else

#ifndef LMP_PPPM2_H

#define LMP_PPPM2_H

#include "kspace.h"

#if defined(FFT_FFTW3)

#define LMP_FFT_LIB "FFTW3"

#elif defined(FFT_MKL)

#define LMP_FFT_LIB "MKL FFT"

#elif defined(FFT_CUFFT)

#define LMP_FFT_LIB "cuFFT"

#else

#define LMP_FFT_LIB "KISS FFT"

#endif

#ifdef FFT_SINGLE

typedef float FFT_SCALAR;

#define LMP_FFT_PREC "single"

#define MPI_FFT_SCALAR MPI_FLOAT

#else

typedef double FFT_SCALAR;

#define LMP_FFT_PREC "double"

#define MPI_FFT_SCALAR MPI_DOUBLE

#endif

namespace LAMMPS_NS {

class PPPM2 : public KSpace {

 public:

  PPPM2(class LAMMPS *);

  virtual ~PPPM2();

  virtual void settings(int, char **);

  virtual void init();

  virtual void setup();

  virtual void setup_grid();

  virtual void compute(int, int);

  virtual int timing_1d(int, double &);

  virtual int timing_3d(int, double &);

  virtual double memory_usage();

  virtual void compute_group_group(int, int, int);

 protected:

  int me,nprocs;

  int nfactors;

  int *factors;

  double cutoff;

  double volume;

  double delxinv,delyinv,delzinv,delvolinv;

  double h_x,h_y,h_z;

  double shift,shiftone;

  int peratom_allocate_flag;

  int nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in;

  int nxlo_out,nylo_out,nzlo_out,nxhi_out,nyhi_out,nzhi_out;

  int nxlo_ghost,nxhi_ghost,nylo_ghost,nyhi_ghost,nzlo_ghost,nzhi_ghost;

  int nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft;

  int nlower,nupper;

  int ngrid,nfft,nfft_both;

  FFT_SCALAR ***density_brick;

  FFT_SCALAR ***vdx_brick,***vdy_brick,***vdz_brick;

  FFT_SCALAR ***u_brick;

  FFT_SCALAR ***v0_brick,***v1_brick,***v2_brick;

  FFT_SCALAR ***v3_brick,***v4_brick,***v5_brick;

  double *greensfn;

  double **vg;

  double *fkx,*fky,*fkz;

  FFT_SCALAR *density_fft;

  FFT_SCALAR *work1,*work2;

  double *gf_b;

  FFT_SCALAR **rho1d,**rho_coeff,**drho1d,**drho_coeff;

  double *sf_precoeff1, *sf_precoeff2, *sf_precoeff3;

  double *sf_precoeff4, *sf_precoeff5, *sf_precoeff6;

  double sf_coeff[6];          // coefficients for calculating ad self-forces

  double **acons;

  // FFTs and grid communication

  class FFT3d *fft1,*fft2;

  class Remap *remap;

  class GridComm2 *gc;

  FFT_SCALAR *gc_buf1,*gc_buf2;

  int ngc_buf1,ngc_buf2,npergrid;

  // group-group interactions

  int group_allocate_flag;

  FFT_SCALAR ***density_A_brick,***density_B_brick;

  FFT_SCALAR *density_A_fft,*density_B_fft;

  int **part2grid;             // storage for particle -> grid mapping

  int nmax;

  double *boxlo;

                               // TIP4P settings

  int typeH,typeO;             // atom types of TIP4P water H and O atoms

  double qdist;                // distance from O site to negative charge

  double alpha;                // geometric factor

  virtual void set_grid_global();

  void set_grid_local();

  void adjust_gewald();

  virtual double newton_raphson_f();

  double derivf();

  double final_accuracy();

  virtual void allocate();

  virtual void allocate_peratom();

  virtual void deallocate();

  virtual void deallocate_peratom();

  int factorable(int);

  double compute_df_kspace();

  double estimate_ik_error(double, double, bigint);

  virtual double compute_qopt();

  virtual void compute_gf_denom();

  virtual void compute_gf_ik();

  virtual void compute_gf_ad();

  void compute_sf_precoeff();

  virtual void particle_map();

  virtual void make_rho();

  virtual void brick2fft();

  virtual void poisson();

  virtual void poisson_ik();

  virtual void poisson_ad();

  virtual void fieldforce();

  virtual void fieldforce_ik();

  virtual void fieldforce_ad();

  virtual void poisson_peratom();

  virtual void fieldforce_peratom();

  void procs2grid2d(int,int,int,int *, int*);

  void compute_rho1d(const FFT_SCALAR &, const FFT_SCALAR &,

                     const FFT_SCALAR &);

  void compute_drho1d(const FFT_SCALAR &, const FFT_SCALAR &,

                     const FFT_SCALAR &);

  void compute_rho_coeff();

  virtual void slabcorr();

  // grid communication

  virtual void pack_forward2(int, void *, int, int *);

  virtual void unpack_forward2(int, void *, int, int *);

  virtual void pack_reverse2(int, void *, int, int *);

  virtual void unpack_reverse2(int, void *, int, int *);

  // triclinic

  int triclinic;               // domain settings, orthog or triclinic

  void setup_triclinic();

  void compute_gf_ik_triclinic();

  void poisson_ik_triclinic();

  void poisson_groups_triclinic();

  // group-group interactions

  virtual void allocate_groups();

  virtual void deallocate_groups();

  virtual void make_rho_groups(int, int, int);

  virtual void poisson_groups(int);

  virtual void slabcorr_groups(int,int,int);

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

   denominator for Hockney-Eastwood Green's function

     of x,y,z = sin(kx*deltax/2), etc

            inf                 n-1

   S(n,k) = Sum  W(k+pi*j)**2 = Sum b(l)*(z*z)**l

           j=-inf               l=0

          = -(z*z)**n /(2n-1)! * (d/dx)**(2n-1) cot(x)  at z = sin(x)

   gf_b = denominator expansion coeffs

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

  inline double gf_denom(const double &x, const double &y,

                         const double &z) const {

    double sx,sy,sz;

    sz = sy = sx = 0.0;

    for (int l = order-1; l >= 0; l--) {

      sx = gf_b[l] + sx*x;

      sy = gf_b[l] + sy*y;

      sz = gf_b[l] + sz*z;

    }

    double s = sx*sy*sz;

    return s*s;

  };

};

}

#endif

#endif

/* ERROR/WARNING messages:

E: Illegal ... command

Self-explanatory.  Check the input script syntax and compare to the

documentation for the command.  You can use -echo screen as a

command-line option when running LAMMPS to see the offending line.

E: Must redefine kspace_style after changing to triclinic box

UNDOCUMENTED

E: Cannot (yet) use PPPM with triclinic box and kspace_modify diff ad

This feature is not yet supported.

E: Cannot (yet) use PPPM with triclinic box and slab correction

This feature is not yet supported.

E: Cannot use PPPM with 2d simulation

The kspace style pppm cannot be used in 2d simulations.  You can use

2d PPPM in a 3d simulation; see the kspace_modify command.

E: PPPM can only currently be used with comm_style brick

This is a current restriction in LAMMPS.

E: Kspace style requires atom attribute q

The atom style defined does not have these attributes.

E: Cannot use non-periodic boundaries with PPPM

For kspace style pppm, all 3 dimensions must have periodic boundaries

unless you use the kspace_modify command to define a 2d slab with a

non-periodic z dimension.

E: Incorrect boundaries with slab PPPM

Must have periodic x,y dimensions and non-periodic z dimension to use

2d slab option with PPPM.

E: PPPM order cannot be < 2 or > than %d

This is a limitation of the PPPM implementation in LAMMPS.

E: KSpace style is incompatible with Pair style

Setting a kspace style requires that a pair style with matching

long-range Coulombic or dispersion components be used.

E: Pair style is incompatible with TIP4P KSpace style

The pair style does not have the requires TIP4P settings.

E: Bond and angle potentials must be defined for TIP4P

Cannot use TIP4P pair potential unless bond and angle potentials

are defined.

E: Bad TIP4P angle type for PPPM/TIP4P

Specified angle type is not valid.

E: Bad TIP4P bond type for PPPM/TIP4P

Specified bond type is not valid.

W: Reducing PPPM order b/c stencil extends beyond nearest neighbor processor

This may lead to a larger grid than desired.  See the kspace_modify overlap

command to prevent changing of the PPPM order.

E: PPPM order < minimum allowed order

The default minimum order is 2.  This can be reset by the

kspace_modify minorder command.

E: PPPM grid stencil extends beyond nearest neighbor processor

This is not allowed if the kspace_modify overlap setting is no.

E: KSpace accuracy must be > 0

The kspace accuracy designated in the input must be greater than zero.

E: Must use kspace_modify gewald for uncharged system

UNDOCUMENTED

E: Could not compute grid size

The code is unable to compute a grid size consistent with the desired

accuracy.  This error should not occur for typical problems.  Please

send an email to the developers.

E: PPPM grid is too large

The global PPPM grid is larger than OFFSET in one or more dimensions.

OFFSET is currently set to 4096.  You likely need to decrease the

requested accuracy.

E: Could not compute g_ewald

The Newton-Raphson solver failed to converge to a good value for

g_ewald.  This error should not occur for typical problems.  Please

send an email to the developers.

E: Non-numeric box dimensions - simulation unstable

The box size has apparently blown up.

E: Out of range atoms - cannot compute PPPM

One or more atoms are attempting to map their charge to a PPPM grid

point that is not owned by a processor.  This is likely for one of two

reasons, both of them bad.  First, it may mean that an atom near the

boundary of a processor's sub-domain has moved more than 1/2 the

"neighbor skin distance"_neighbor.html without neighbor lists being

rebuilt and atoms being migrated to new processors.  This also means

you may be missing pairwise interactions that need to be computed.

The solution is to change the re-neighboring criteria via the

"neigh_modify"_neigh_modify command.  The safest settings are "delay 0

every 1 check yes".  Second, it may mean that an atom has moved far

outside a processor's sub-domain or even the entire simulation box.

This indicates bad physics, e.g. due to highly overlapping atoms, too

large a timestep, etc.

E: Cannot (yet) use K-space slab correction with compute group/group for triclinic systems

This option is not yet supported.

E: Cannot (yet) use kspace_modify diff ad with compute group/group

This option is not yet supported.

U: Cannot (yet) use PPPM with triclinic box and TIP4P

This feature is not yet supported.

*/

  kspace = new_kspace(style,trysuffix,sflag);

  store_style(kspace_style,style,sflag);

  if (comm->style == 1 && !kspace_match("ewald",0))

    error->all(FLERR,

               "Cannot yet use KSpace solver with grid with comm style tiled");

  //if (comm->style == 1 && !kspace_match("ewald",0))

  // error->all(FLERR,

  //             "Cannot yet use KSpace solver with grid with comm style tiled");

}

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