[WIP] Adding Ewald+SPME longrange force solvers

CMakeLists.txt

@@ -59,6 +59,11 @@ if(NOT DEFINED KOKKOS_DIR)
   set(KOKKOS_DIR ~/install/kokkos/)
 endif()
 
+option(CabanaMD_ENABLE_LongRange "Build CabanaMD with longrange forces (Ewald)" OFF)
+if( CabanaMD_ENABLE_LongRange )
+  add_definitions(-DCabanaMD_ENABLE_LongRange=ON)
+endif()
+
 ##---------------------------------------------------------------------------##
 ## Print the revision number to stdout
 ##---------------------------------------------------------------------------##
@@ -78,8 +83,13 @@ MESSAGE(STATUS "CabanaMD Revision = '${CabanaMD_GIT_COMMIT_HASH}'")
 ##---------------------------------------------------------------------------##
 ## Build CabanaMD
 ##---------------------------------------------------------------------------##
-include_directories(src src/force_types)
-file(GLOB SOURCES src/*.cpp src/force_types/*.cpp)
+if( CabanaMD_ENABLE_LongRange )
+  include_directories(src src/force_types src/longrange_force_types)
+  file(GLOB SOURCES src/*.cpp src/force_types/*.cpp src/longrange_force_types/*.cpp)
+else()
+  include_directories(src src/force_types)
+  file(GLOB SOURCES src/*.cpp src/force_types/*.cpp)
+endif()
 
 include_directories(SYSTEM ${KOKKOS_DIR}/include)
 link_directories(${KOKKOS_DIR}/lib)
@@ -89,4 +99,5 @@ link_directories(${CABANA_DIR}/lib)
 
 add_executable(CabanaMD ${SOURCES})
 target_link_libraries(CabanaMD cabanacore kokkos dl MPI::MPI_CXX)
+
 install(TARGETS CabanaMD DESTINATION bin)

src/cabanamd.cpp

@@ -98,6 +98,7 @@ void CabanaMD::init(int argc, char* argv[]) {
     force->init_coeff(neigh_cutoff,
                       input->input_data.words[input->force_coeff_lines(line)]);
   }
+  //TODO: make sure longrange force can also be grabbed from input deck
 
   // Create Communication Submodule
   comm = new Comm(system, neigh_cutoff);
@@ -107,7 +108,10 @@ void CabanaMD::init(int argc, char* argv[]) {
 
   // system->print_particles();
   if(system->do_print) {
-    printf("Using: %s %s %s %s\n",force->name(),comm->name(),binning->name(),integrator->name());
+    if(input->longrange)
+      printf("Using: %s %s %s %s %s\n",force->name(),lrforce->name(),comm->name(),binning->name(),integrator->name());
+    else
+      printf("Using: %s %s %s %s\n",force->name(),comm->name(),binning->name(),integrator->name());
   }
 
   // Ok lets go ahead and create the particles if that didn't happen yet
@@ -132,6 +136,8 @@ void CabanaMD::init(int argc, char* argv[]) {
   auto f = Cabana::slice<Forces>(system->xvf);
   Cabana::deep_copy(f, 0.0);
   force->compute(system);
+  if(input->longrange)
+    lrforce->compute(system);
 
   if(input->comm_newton) {
     // Reverse Communicate Force Update on Halo
@@ -146,6 +152,7 @@ void CabanaMD::init(int argc, char* argv[]) {
     KinE kine(comm);
     T_FLOAT T = temp.compute(system);
     T_FLOAT PE = pote.compute(system,force)/system->N;
+    //TODO: add longrange contribution to PE
     T_FLOAT KE = kine.compute(system)/system->N;
     if(system->do_print) {
       if (!system->print_lammps) {
@@ -176,13 +183,14 @@ void CabanaMD::run(int nsteps) {
   KinE kine(comm);
 
   double force_time = 0;
+  double lrforce_time = 0;
   double comm_time  = 0;
   double neigh_time = 0;
   double integrate_time = 0;
   double other_time = 0;
 
   double last_time = 0;
-  Kokkos::Timer timer,force_timer,comm_timer,neigh_timer,integrate_timer,other_timer;
+  Kokkos::Timer timer,force_timer,lrforce_timer,comm_timer,neigh_timer,integrate_timer,other_timer;
 
   // Timestep Loop
   for(int step = 1; step <= nsteps; step++ ) {
@@ -228,7 +236,14 @@ void CabanaMD::run(int nsteps) {
     force->compute(system);
     force_time += force_timer.seconds();
 
-    // This is where Bonds, Angles and KSpace should go eventually 
+    //Compute Long Range Force
+    if(input->longrange) {
+      lrforce_timer.reset();
+      lrforce->compute(system);
+      lrforce_time += lrforce_timer.seconds();
+    }
+
+    // This is where Bonds, Angles, etc. should go eventually
 
     // Reverse Communicate Force Update on Halo
     if(input->comm_newton) {

src/cabanamd.h

@@ -51,18 +51,19 @@
 #include <Cabana_Core.hpp>
 
 #include <types.h>
-#include<system.h>
-#include<integrator_nve.h>
-#include<force.h>
-#include<comm_mpi.h>
-#include<input.h>
-#include<binning_cabana.h>
+#include <system.h>
+#include <integrator_nve.h>
+#include <force.h>
+#include <comm_mpi.h>
+#include <input.h>
+#include <binning_cabana.h>
 
 class CabanaMD {
   public:
     System* system;
     Integrator* integrator;
     Force* force;
+    Force* lrforce;
     Comm* comm;
     Input* input;
     Binning* binning;

src/force.h

@@ -56,6 +56,7 @@ class Force {
 public:
   bool half_neigh, comm_newton;
   Force(System* system, bool half_neigh_);
+  Force(System* system, bool half_neigh_, MPI_Comm comm);
 
   virtual void init_coeff(T_X_FLOAT neigh_cut, char** args);
   virtual void create_neigh_list(System* system);

src/input.cpp

@@ -130,12 +130,19 @@ void ItemizedFile::add_line(const char* const line) {
 
 Input::Input(System* p):system(p),input_data(ItemizedFile()),data_file_data(ItemizedFile()) {
 
+  binning_type = BINNING_LINKEDCELL;
   comm_type = COMM_MPI;
   integrator_type = INTEGRATOR_NVE;
   neighbor_type = NEIGH_2D;
   force_type = FORCE_LJ;
   force_iteration_type = FORCE_ITER_NEIGH_FULL;
-  binning_type = BINNING_LINKEDCELL;
+  lrforce_type = FORCE_NONE;
+  lrforce_iteration_type = force_iteration_type;
+
+  longrange = false;
+  #ifdef CabanaMD_ENABLE_LongRange
+  longrange = true;
+  #endif
 
   // set defaults (matches ExaMiniMD LJ example)
 

src/input.h

@@ -182,7 +182,10 @@ class Input {
   int comm_newton;
 
   int force_type;
+  int lrforce_type;
   int force_iteration_type;
+  int lrforce_iteration_type;
+  bool longrange;
 
   T_F_FLOAT force_cutoff;
   int force_line;

src/longrange_force_types/force_ewald_cabana_neigh.cpp

@@ -0,0 +1,18 @@
+/****************************************************************************
+ * Copyright (c) 2018-2019 by the Cabana authors                            *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the Cabana library. Cabana is distributed under a   *
+ * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
+ * the top-level directory.                                                 *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+#include<force_ewald_cabana_neigh_impl.h>
+
+
+template class ForceEwald<t_verletlist_half_2D>;
+template class ForceEwald<t_verletlist_full_2D>;
+template class ForceEwald<t_verletlist_half_CSR>;
+template class ForceEwald<t_verletlist_full_CSR>;

src/longrange_force_types/force_ewald_cabana_neigh.h

@@ -0,0 +1,105 @@
+/****************************************************************************
+ * Copyright (c) 2018-2019 by the Cabana authors                            *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the Cabana library. Cabana is distributed under a   *
+ * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
+ * the top-level directory.                                                 *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+
+//TODO: I don't really understand this first part
+#ifdef MODULES_OPTION_CHECK
+  if( (strcmp(argv[i], "--force-iteration") == 0) ) {
+    if( (strcmp(argv[i+1], "NEIGH_FULL") == 0) )
+      lrforce_iteration_type = FORCE_ITER_NEIGH_FULL;
+    if( (strcmp(argv[i+1], "NEIGH_HALF") == 0) )
+      lrforce_iteration_type = FORCE_ITER_NEIGH_HALF;
+  }
+  if( (strcmp(argv[i], "--neigh-type") == 0) ) {
+    if( (strcmp(argv[i+1], "NEIGH_2D") == 0) )
+      neighbor_type = NEIGH_2D;
+    if( (strcmp(argv[i+1], "NEIGH_CSR") == 0) )
+      neighbor_type = NEIGH_CSR;
+  }
+#endif
+#ifdef FORCE_MODULES_INSTANTIATION
+    else if (input->lrforce_type == FORCE_EWALD) {
+      bool half_neigh = input->lrforce_iteration_type == FORCE_ITER_NEIGH_HALF;
+      if (input->neighbor_type == NEIGH_2D) {
+        if (half_neigh)
+          lrforce = new ForceEwald<t_verletlist_half_2D>(system,half_neigh);
+        else
+          lrforce = new ForceEwald<t_verletlist_full_2D>(system,half_neigh);
+      }
+      else if (input->neighbor_type == NEIGH_CSR) {
+	if (half_neigh)
+          lrforce = new ForceEwald<t_verletlist_half_CSR>(system,half_neigh);
+        else
+          lrforce = new ForceEwald<t_verletlist_full_CSR>(system,half_neigh);
+      }
+      #undef FORCETYPE_ALLOCATION_MACRO
+    }
+#endif
+
+
+#if !defined(MODULES_OPTION_CHECK) && \
+    !defined(FORCE_MODULES_INSTANTIATION)
+
+#ifndef FORCE_EWALD_CABANA_NEIGH_H
+#define FORCE_EWALD_CABANA_NEIGH_H
+#include <Cabana_Core.hpp>
+#include<force.h>
+#include<types.h>
+#include<system.h>
+#include <comm_mpi.h>
+
+template<class t_neighbor>
+class ForceEwald: public Force {
+private:
+  int N_local,ntypes;
+  typename AoSoA::member_slice_type<Positions> x;
+  typename AoSoA::member_slice_type<Forces> f;
+  typename AoSoA::member_slice_type<Forces>::atomic_access_slice f_a;
+  typename AoSoA::member_slice_type<IDs> id;
+  typename AoSoA::member_slice_type<Types> type;
+  typename AoSoA::member_slice_type<Charges> q;
+  typename AoSoA::member_slice_type<Potentials> p;
+
+  double _alpha;
+  double _r_max;
+  double _k_max;
+
+  //dielectric constant
+  double _eps_r = 1.0; //Assume 1 for now (vacuum)
+
+  double *EwaldUk_coeffs;
+
+  //Kokkos::View<double *, MemorySpace> domain_width;
+
+  MPI_Comm comm;
+
+public:
+
+  bool half_neigh;
+  T_X_FLOAT neigh_cut;
+
+  t_neighbor neigh_list;
+
+  ForceEwald(System* system, bool half_neigh);
+  ForceEwald(System* system, bool half_neigh, Comm comm);
+
+  void init_coeff(char** args);
+
+  //void create_neigh_list(System* system);
+
+  void compute(System* system);
+  T_F_FLOAT compute_energy(System* system);
+
+  const char* name();
+};
+
+#endif
+#endif