MPI support for SGP construction

.gitignore

@@ -10,3 +10,4 @@ _C_flare*
 **/xml
 **/dist
 **egg-info
+*.DS_Store

CMakeLists.txt

@@ -9,7 +9,7 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/)
 
 include(ExternalProject)
 
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
@@ -111,6 +111,22 @@ endif()
 
 # Greatly reduces the code bloat
 target_compile_options(pybind11 INTERFACE "-fvisibility=hidden")
+
+# Distmatrix
+###############################################################################
+FetchContent_Declare(
+    distmatrix
+    GIT_REPOSITORY https://github.com/mir-group/distmatrix.git
+)
+FetchContent_MakeAvailable(distmatrix)
+
+link_directories(${distmatrix_BINARY_DIR})
+link_directories(${distmatrix_BINARY_DIR}/scalapack_build/lib)
+
+include_directories(${distmatrix_SOURCE_DIR}/include)
+
+
+#target_link_libraries(flare PUBLIC distmatrix)
 ###############################################################################
 
 # Specify source files.
@@ -119,14 +135,19 @@ set(FLARE_SOURCES
     src/flare_pp/radial.cpp
     src/flare_pp/cutoffs.cpp
     src/flare_pp/structure.cpp
+    src/flare_pp/utils.cpp
     src/flare_pp/bffs/sparse_gp.cpp
+    src/flare_pp/bffs/parallel_sgp.cpp
     src/flare_pp/bffs/gp.cpp
     src/flare_pp/descriptors/descriptor.cpp
+    src/flare_pp/descriptors/b1.cpp
     src/flare_pp/descriptors/b2.cpp
     src/flare_pp/descriptors/b2_norm.cpp
     src/flare_pp/descriptors/b2_simple.cpp
     src/flare_pp/descriptors/b3.cpp
-    src/flare_pp/descriptors/wigner3j.cpp
+    src/flare_pp/descriptors/bk.cpp
+    src/flare_pp/descriptors/indices.cpp
+    src/flare_pp/descriptors/coeffs.cpp
     src/flare_pp/descriptors/two_body.cpp
     src/flare_pp/descriptors/three_body.cpp
     src/flare_pp/descriptors/three_body_wide.cpp
@@ -152,6 +173,14 @@ include_directories(
 add_library(flare STATIC ${FLARE_SOURCES})
 set_target_properties(flare PROPERTIES POSITION_INDEPENDENT_CODE ON)
 add_dependencies(flare eigen_project)
+add_dependencies(flare distmatrix)
+
+# Link to MPI
+find_package(MPI REQUIRED)
+target_link_libraries(flare PUBLIC MPI::MPI_CXX)
+
+# Link to distmatrix
+target_link_libraries(flare PUBLIC distmatrix -lgfortran)
 
 # Link to json.
 target_link_libraries(flare PUBLIC nlohmann_json::nlohmann_json)
@@ -175,6 +204,8 @@ if (DEFINED ENV{MKL_INCLUDE} AND DEFINED ENV{MKL_LIBS})
   include_directories($ENV{MKL_INCLUDE})
   target_link_libraries(flare PUBLIC $ENV{MKL_LIBS})
   target_compile_definitions(flare PUBLIC EIGEN_USE_MKL_ALL=1)
+  target_include_directories(flare PUBLIC $ENV{MKL_MPI})
+  message(MKL_MPI=$ENV{MKL_MPI})
 # Give the option to do a "bare-bones" build without Lapack/Lapacke.
 elseif(DEFINED ENV{NO_LAPACK})
   message(STATUS "Building without Lapack.")

README.md

@@ -13,3 +13,36 @@ If you're installing on Harvard's compute cluster, load the following modules fi
 ```
 module load cmake/3.17.3-fasrc01 python/3.6.3-fasrc01 gcc/9.3.0-fasrc01
 ```
+
+### MPI compilation
+To install the MPI version of sparse GP, on Harvard's compute cluster, load the following modules first:
+```
+module load cmake/3.17.3-fasrc01 python/3.6.3-fasrc01 gcc/9.3.0-fasrc01
+module load intel-mkl/2017.2.174-fasrc01 openmpi/4.0.5-fasrc01
+```
+
+Then clone the repository
+```
+git clone https://github.com/mir-group/flare_pp.git
+git checkout mpi_distmat
+```
+
+Create a directory for building the library
+```
+mkdir build
+cd build
+```
+
+Compile with `cmake` and `make`. Here we need to specify the MPI compiler with 
+`CC=mpicc CXX=mpic++ FC=mpif90`. We also set the option `-DSCALAPACK_LIB=NOTFOUND`,
+such that we download and compile our own static Scalapack library to support 
+python binding with MPI parallelized sparse GP
+```
+CC=mpicc CXX=mpic++ FC=mpif90 cmake .. -DSCALAPACK_LIB=NOTFOUND
+make -j
+```
+
+Copy the python binding library to the folder of `flare_pp`
+```
+cp _C_flare.cpython-36m-x86_64-linux-gnu.so ../flare_pp
+```

ctests/CMakeLists.txt

@@ -1,3 +1,8 @@
+if (NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/dft_data.xyz)
+    file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/dft_data.xyz
+                DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+endif ()
+
 # make executable
 add_executable(
   tests
@@ -6,8 +11,10 @@ add_executable(
   test_structure.cpp
   test_n_body.cpp
   test_sparse_gp.cpp
+  test_parallel_sgp.cpp
   test_descriptor.cpp
   test_json.cpp
+  test_utils.cpp
 )
 
 include_directories(../src/flare_pp)

ctests/dft_data.xyz

@@ -0,0 +1,17 @@
+2
+Lattice="1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" Properties=species:S:1:pos:R:3:momenta:R:3:forces:R:3 sparse_indices= energy=0.9133510291442271 stress="0.4314439934824156 0.7292230215429331 0.38721540888657746 0.7292230215429331 0.326605318720777 0.25504993852246405 0.38721540888657746 0.25504993852246405 0.3185010782409763" pbc="T T T"
+He       0.15713841       0.99710398       0.64092684       0.19861472       0.12959016       0.46411878       0.48801624       0.78256914       0.52750114
+He       0.69338569       0.90346464       0.14626089       0.92338744       0.37879257       0.79694737       0.92320467       0.71863087       0.54575622
+4
+Lattice="1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" Properties=species:S:1:pos:R:3:momenta:R:3:forces:R:3 sparse_indices=2 energy=0.5845339958967781 stress="0.8832643728374101 0.4736589482749042 0.6125575360846077 0.4736589482749042 0.04982809557310608 0.9519383879956919 0.6125575360846077 0.9519383879956919 0.22321884060680475" pbc="T T T"
+H        0.98699548       0.27109152       0.82840818       0.39723635       0.69230317       0.71742402       0.97501125       0.69790992       0.42736135
+He       0.95542741       0.62792707       0.71255717       0.20846441       0.71391509       0.56988269       0.73802987       0.55611127       0.20978028
+H        0.95551793       0.79583919       0.01626457       0.50770084       0.31544704       0.09250715       0.97886919       0.15580748       0.02892780
+He       0.34877543       0.15749443       0.15293379       0.64543458       0.63466428       0.15787333       0.67504974       0.30204244       0.41673029
+5
+Lattice="1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" Properties=species:S:1:pos:R:3:momenta:R:3:forces:R:3 sparse_indices="0 2 4" energy=0.26796364806180983 stress="0.5706389111071576 0.3483258644776196 0.06089244456604703 0.3483258644776196 0.5296282371894802 0.21476826115341807 0.06089244456604703 0.21476826115341807 0.008875719863243536" pbc="T T T"
+H        0.35071505       0.73963362       0.02696195       0.50938856       0.68674285       0.31711543       0.21822063       0.99866869       0.23949716
+H        0.78288668       0.67299943       0.47254490       0.03722967       0.36556295       0.28723150       0.36477796       0.28274931       0.78633052
+H        0.43430905       0.05355055       0.52035306       0.40526494       0.75298417       0.54189491       0.15815742       0.15396677       0.20679348
+He       0.50335965       0.87700387       0.47565834       0.59881725       0.71577512       0.87647597       0.44898467       0.56835792       0.32339863
+He       0.28826927       0.14879537       0.18758152       0.00333732       0.24138565       0.27410589       0.43239731       0.43696276       0.94641792

ctests/test_descriptor.cpp

@@ -1,52 +1,132 @@
+#include "bk.h"
+#include "b1.h"
+#include "b2.h"
 #include "b3.h"
 #include "descriptor.h"
 #include "test_structure.h"
 #include "gtest/gtest.h"
 #include <Eigen/Dense>
 #include <cmath>
+#include <list>
 #include <iostream>
 
-TEST_F(StructureTest, RotationTest) {
+// Test different types B1, B2, B3 to match with Bk
+template <typename T>
+class DescTest : public StructureTest {
+public:
+  using List = std::list<T>;
+};
+
+using DescTypes = ::testing::Types<B1, B2, B3>;
+TYPED_TEST_SUITE(DescTest, DescTypes);
+
+TYPED_TEST(DescTest, TestBk) {
+  // Set up B1/2/3 descriptors
+  std::vector<int> descriptor_settings_1{this->n_species, this->N, this->L};
+  std::vector<Descriptor *> descriptors1;
+
+  TypeParam desc1(this->radial_string, this->cutoff_string, this->radial_hyps, 
+                  this->cutoff_hyps, descriptor_settings_1);
+  descriptors1.push_back(&desc1);
+
+  // Set up Bk descriptors
+  int K = desc1.K;
+  std::vector<int> descriptor_settings_2{this->n_species, K, this->N, this->L};
+  Bk desc2(this->radial_string, this->cutoff_string, this->radial_hyps, 
+           this->cutoff_hyps, descriptor_settings_2);
+  std::vector<Descriptor *> descriptors2;
+  descriptors2.push_back(&desc2);
+
+  Structure struc1 = Structure(this->cell, this->species, this->positions, this->cutoff, descriptors1);
+  Structure struc2 = Structure(this->cell, this->species, this->positions, this->cutoff, descriptors2);
+
+  // Check the descriptor dimensions
+  std::vector<int> last_index = desc2.nu[desc2.nu.size()-1];
+  int n_d = last_index[last_index.size()-1] + 1; // the size of list nu
+  int n_d1 = struc1.descriptors[0].n_descriptors;
+  int n_d2 = struc2.descriptors[0].n_descriptors;
+  EXPECT_EQ(n_d, n_d1);
+  EXPECT_EQ(n_d1, n_d2);
+
+  // Check that Bk and B3 give the same descriptors.
+  double d1, d2;
+  int nu_ind;
+  for (int i = 0; i < struc1.descriptors.size(); i++) {
+    for (int j = 0; j < struc1.descriptors[i].descriptors.size(); j++) {
+      for (int k = 0; k < struc1.descriptors[i].descriptors[j].rows(); k++) {
+        for (int l = 0; l < struc1.descriptors[i].descriptors[j].cols(); l++) {
+          d1 = struc1.descriptors[i].descriptors[j](k, l);
+          d2 = struc2.descriptors[i].descriptors[j](k, l);
+          EXPECT_NEAR(d1, d2, 1e-8);
+        }
+      }
+    }
+  }
+
+  for (int i = 0; i < struc1.descriptors.size(); i++) {
+    for (int j = 0; j < struc1.descriptors[i].descriptor_force_dervs.size(); j++) {
+      for (int k = 0; k < struc1.descriptors[i].descriptor_force_dervs[j].rows(); k++) {
+        for (int l = 0; l < struc1.descriptors[i].descriptor_force_dervs[j].cols(); l++) {
+          d1 = struc1.descriptors[i].descriptor_force_dervs[j](k, l);
+          d2 = struc2.descriptors[i].descriptor_force_dervs[j](k, l);
+          EXPECT_NEAR(d1, d2, 1e-8);
+        }
+      }
+    }
+  }
+
+}
 
+// Test Bk with K=1,2,3 for rotational invariance 
+class DescRotTest : public StructureTest,
+                    public testing::WithParamInterface<int> {
+public:
   // Choose arbitrary rotation angles.
   double xrot = 1.28;
   double yrot = -3.21;
   double zrot = 0.42;
 
+  Eigen::MatrixXd rotated_pos;
+  Eigen::MatrixXd rotated_cell;
+
   // Define rotation matrices.
   Eigen::MatrixXd Rx{3, 3}, Ry{3, 3}, Rz{3, 3}, R{3, 3};
-  Rx << 1, 0, 0, 0, cos(xrot), -sin(xrot), 0, sin(xrot), cos(xrot);
-  Ry << cos(yrot), 0, sin(yrot), 0, 1, 0, -sin(yrot), 0, cos(yrot);
-  Rz << cos(zrot), -sin(zrot), 0, sin(zrot), cos(zrot), 0, 0, 0, 1;
-  R = Rx * Ry * Rz;
-
-  Eigen::MatrixXd rotated_pos = positions * R.transpose();
-  Eigen::MatrixXd rotated_cell = cell * R.transpose();
+  DescRotTest() {
+    Rx << 1, 0, 0, 0, cos(xrot), -sin(xrot), 0, sin(xrot), cos(xrot);
+    Ry << cos(yrot), 0, sin(yrot), 0, 1, 0, -sin(yrot), 0, cos(yrot);
+    Rz << cos(zrot), -sin(zrot), 0, sin(zrot), cos(zrot), 0, 0, 0, 1;
+    R = Rx * Ry * Rz;
+    rotated_pos = positions * R.transpose();
+    rotated_cell = cell * R.transpose();
+  }
+};
 
-  // Define descriptors.
-  descriptor_settings[2] = 2;
-  B3 descriptor = B3(radial_string, cutoff_string, radial_hyps, cutoff_hyps,
-                     descriptor_settings);
+TEST_P(DescRotTest, RotationTest) {
+  int K = GetParam();
+  std::vector<int> descriptor_settings{n_species, K, N, L};
+  Bk desc(radial_string, cutoff_string, radial_hyps, cutoff_hyps,
+          descriptor_settings);
 
   std::vector<Descriptor *> descriptors;
-  descriptors.push_back(&descriptor);
+  descriptors.push_back(&desc);
 
   Structure struc1 = Structure(cell, species, positions, cutoff, descriptors);
   Structure struc2 =
-      Structure(rotated_cell, species, rotated_pos, cutoff, descriptors);
+      Structure(this->rotated_cell, species, this->rotated_pos, cutoff, descriptors);
 
   // Check that B1 is rotationally invariant.
-  double d1, d2, diff;
-  double tol = 1e-10;
+  double d1, d2;
 
+  std::cout << "n_descriptors=" << struc1.descriptors[0].n_descriptors << std::endl;
   for (int n = 0; n < struc1.descriptors[0].n_descriptors; n++) {
     d1 = struc1.descriptors[0].descriptors[0](0, n);
     d2 = struc2.descriptors[0].descriptors[0](0, n);
-    diff = d1 - d2;
-    EXPECT_LE(abs(diff), tol);
+    EXPECT_NEAR(d1, d2, 1e-10);
   }
 }
 
+INSTANTIATE_TEST_SUITE_P(DescBodies, DescRotTest, testing::Values(1, 2, 3));
+
   // TEST_F(DescriptorTest, SingleBond) {
   //   // Check that B1 descriptors match the corresponding elements of the
   //   // single bond vector.