From 85b8ededdf85650e52924fb7711e45f13d80bdbd Mon Sep 17 00:00:00 2001
From: Carl Poelking <cp605@cam.ac.uk>
Date: Tue, 26 Apr 2016 17:48:43 +0100
Subject: [PATCH] Build system mods, cxx_tests, NablaYlm.

---
 CMakeLists.txt                           |   3 +-
 cxx_tests/CMakeLists.txt                 |  61 ++++
 cxx_tests/build.sh                       |   5 +
 cxx_tests/cmake/FindFFTW3.cmake          |  22 ++
 cxx_tests/cmake/FindGSL.cmake            |  59 ++++
 cxx_tests/cmake/FindGTEST.cmake          |  11 +
 cxx_tests/cmake/FindMKL.cmake            | 339 +++++++++++++++++++++++
 cxx_tests/cmake/FindSOAP.cmake           |  11 +
 cxx_tests/gtest.sh                       |   4 +
 cxx_tests/src/CMakeLists.txt             |  17 ++
 cxx_tests/src/gtest_dylm.cpp             | 190 +++++++++++++
 cxx_tests/src/test_dylm.cpp              |  52 ++++
 cxx_tests/src/test_dylm.py               | 116 ++++++++
 cxx_tests/src/test_functions.cpp         |  26 ++
 cxx_tests/src/test_functions.py          |  10 +
 src/functions.cpp                        |  95 -------
 src/{ => soap}/CMakeLists.txt            |  10 +-
 src/{ => soap}/SOAPRC.in                 |   4 +
 src/{ => soap}/__init__.py               |   0
 src/{ => soap}/angularbasis.cpp          |   2 +-
 src/{ => soap}/angularbasis.hpp          |   6 +-
 src/{ => soap}/atomicspectrum.cpp        |   2 +-
 src/{ => soap}/atomicspectrum.hpp        |  14 +-
 src/{ => soap}/base/exceptions.hpp       |   0
 src/{ => soap}/base/logger.hpp           |   0
 src/{ => soap}/base/objectfactory.hpp    |   0
 src/{ => soap}/basis.cpp                 |   4 +-
 src/{ => soap}/basis.hpp                 |  14 +-
 src/{ => soap}/bindings.cpp              |   0
 src/{ => soap}/bindings.hpp              |   6 +-
 src/{ => soap}/boundary.cpp              |   2 +-
 src/{ => soap}/boundary.hpp              |   2 +-
 src/{ => soap}/cutoff.cpp                |   2 +-
 src/{ => soap}/cutoff.hpp                |   8 +-
 src/soap/functions.cpp                   | 247 +++++++++++++++++
 src/{ => soap}/functions.hpp             |  19 +-
 src/{ => soap}/globals.cpp               |   2 +-
 src/{ => soap}/globals.hpp               |   2 +-
 src/{ => soap}/linalg/CMakeLists.txt     |   3 +
 src/{ => soap}/linalg/__init__.py        |   0
 src/{ => soap}/linalg/bindings.cpp       |   2 +-
 src/{ => soap}/linalg/bindings.hpp       |   4 +-
 src/{ => soap}/linalg/gsl/operations.cpp |   2 +-
 src/{ => soap}/linalg/matrix.cpp         |   2 +-
 src/{ => soap}/linalg/matrix.hpp         |   2 +-
 src/{ => soap}/linalg/mkl/operations.cpp |   2 +-
 src/{ => soap}/linalg/numpy.hpp          |   0
 src/{ => soap}/linalg/operations.hpp     |   0
 src/{ => soap}/linalg/types.hpp          |   2 +-
 src/{ => soap}/linalg/vec.hpp            |   0
 src/{ => soap}/options.cpp               |   2 +-
 src/{ => soap}/options.hpp               |   2 +-
 src/{ => soap}/power.cpp                 |   4 +-
 src/{ => soap}/power.hpp                 |   0
 src/{ => soap}/radialbasis.cpp           |   8 +-
 src/{ => soap}/radialbasis.hpp           |   8 +-
 src/{ => soap}/spectrum.cpp              |   2 +-
 src/{ => soap}/spectrum.hpp              |  16 +-
 src/{ => soap}/structure.cpp             |   2 +-
 src/{ => soap}/structure.hpp             |   6 +-
 src/{ => soap}/tools/CMakeLists.txt      |   0
 src/{ => soap}/tools/__init__.py         |   0
 src/{ => soap}/tools/extract.py          |   0
 src/{ => soap}/tools/inverse.py          |   0
 src/{ => soap}/tools/loadwrite.py        |   0
 src/{ => soap}/tools/loadwrite.py~       |   0
 src/{ => soap}/types.hpp                 |   4 +-
 67 files changed, 1274 insertions(+), 166 deletions(-)
 create mode 100644 cxx_tests/CMakeLists.txt
 create mode 100755 cxx_tests/build.sh
 create mode 100644 cxx_tests/cmake/FindFFTW3.cmake
 create mode 100644 cxx_tests/cmake/FindGSL.cmake
 create mode 100644 cxx_tests/cmake/FindGTEST.cmake
 create mode 100644 cxx_tests/cmake/FindMKL.cmake
 create mode 100644 cxx_tests/cmake/FindSOAP.cmake
 create mode 100755 cxx_tests/gtest.sh
 create mode 100644 cxx_tests/src/CMakeLists.txt
 create mode 100644 cxx_tests/src/gtest_dylm.cpp
 create mode 100644 cxx_tests/src/test_dylm.cpp
 create mode 100644 cxx_tests/src/test_dylm.py
 create mode 100644 cxx_tests/src/test_functions.cpp
 create mode 100644 cxx_tests/src/test_functions.py
 delete mode 100644 src/functions.cpp
 rename src/{ => soap}/CMakeLists.txt (82%)
 rename src/{ => soap}/SOAPRC.in (56%)
 rename src/{ => soap}/__init__.py (100%)
 rename src/{ => soap}/angularbasis.cpp (96%)
 rename src/{ => soap}/angularbasis.hpp (93%)
 rename src/{ => soap}/atomicspectrum.cpp (99%)
 rename src/{ => soap}/atomicspectrum.hpp (92%)
 rename src/{ => soap}/base/exceptions.hpp (100%)
 rename src/{ => soap}/base/logger.hpp (100%)
 rename src/{ => soap}/base/objectfactory.hpp (100%)
 rename src/{ => soap}/basis.cpp (99%)
 rename src/{ => soap}/basis.hpp (91%)
 rename src/{ => soap}/bindings.cpp (100%)
 rename src/{ => soap}/bindings.hpp (58%)
 rename src/{ => soap}/boundary.cpp (92%)
 rename src/{ => soap}/boundary.hpp (99%)
 rename src/{ => soap}/cutoff.cpp (95%)
 rename src/{ => soap}/cutoff.hpp (93%)
 create mode 100644 src/soap/functions.cpp
 rename src/{ => soap}/functions.hpp (66%)
 rename src/{ => soap}/globals.cpp (83%)
 rename src/{ => soap}/globals.hpp (86%)
 rename src/{ => soap}/linalg/CMakeLists.txt (91%)
 rename src/{ => soap}/linalg/__init__.py (100%)
 rename src/{ => soap}/linalg/bindings.cpp (88%)
 rename src/{ => soap}/linalg/bindings.hpp (63%)
 rename src/{ => soap}/linalg/gsl/operations.cpp (98%)
 rename src/{ => soap}/linalg/matrix.cpp (99%)
 rename src/{ => soap}/linalg/matrix.hpp (99%)
 rename src/{ => soap}/linalg/mkl/operations.cpp (98%)
 rename src/{ => soap}/linalg/numpy.hpp (100%)
 rename src/{ => soap}/linalg/operations.hpp (100%)
 rename src/{ => soap}/linalg/types.hpp (77%)
 rename src/{ => soap}/linalg/vec.hpp (100%)
 rename src/{ => soap}/options.cpp (98%)
 rename src/{ => soap}/options.hpp (98%)
 rename src/{ => soap}/power.cpp (98%)
 rename src/{ => soap}/power.hpp (100%)
 rename src/{ => soap}/radialbasis.cpp (98%)
 rename src/{ => soap}/radialbasis.hpp (95%)
 rename src/{ => soap}/spectrum.cpp (99%)
 rename src/{ => soap}/spectrum.hpp (95%)
 rename src/{ => soap}/structure.cpp (99%)
 rename src/{ => soap}/structure.hpp (98%)
 rename src/{ => soap}/tools/CMakeLists.txt (100%)
 rename src/{ => soap}/tools/__init__.py (100%)
 rename src/{ => soap}/tools/extract.py (100%)
 rename src/{ => soap}/tools/inverse.py (100%)
 rename src/{ => soap}/tools/loadwrite.py (100%)
 rename src/{ => soap}/tools/loadwrite.py~ (100%)
 rename src/{ => soap}/types.hpp (93%)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1cccf34..efb0587 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,6 +23,7 @@ endif(${CMAKE_VERSION} VERSION_GREATER 3.1)
 
 # LOCAL PACKAGES
 # ...
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src)
 
 # FIND PACKAGES
 find_package(PythonLibs)
@@ -57,5 +58,5 @@ foreach(lib ${local_dirs})
 endforeach()
 
 # SUBDIRECTORIES
-add_subdirectory(src)
+add_subdirectory(src/soap)
 
diff --git a/cxx_tests/CMakeLists.txt b/cxx_tests/CMakeLists.txt
new file mode 100644
index 0000000..25b3439
--- /dev/null
+++ b/cxx_tests/CMakeLists.txt
@@ -0,0 +1,61 @@
+# CONFIGURE CMAKE
+message("CMake version: ${CMAKE_VERSION}")
+cmake_minimum_required(VERSION 2.8.3)
+
+# PROJECT OPTIONS
+project(soapxx_tests)
+set(CMAKE_INSTALL_PREFIX ${CMAKE_SOURCE_DIR})
+set(LOCAL_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/soap_tests)
+set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
+
+# BUILD OPTIONS
+enable_language(CXX)
+#set(CMAKE_CXX_COMPILER "/usr/local/shared/intel/compilers_and_libraries_2016.0.109/linux/bin/intel64/icc")
+message("C++ compiler: " ${CMAKE_CXX_COMPILER} " " ${CMAKE_CXX_COMPILER_ID})
+option(BUILD_SHARED_LIBS "Build shared libs" ON)
+if(${CMAKE_VERSION} VERSION_GREATER 3.1)
+    message("Setting C++ standard 11 (CMake version > 3.1)")
+    set(CMAKE_CXX_STANDARD 11)
+else(${CMAKE_VERSION} VERSION_GREATER 3.1)
+    message("Setting C++ standard 11 (CMake version <= 3.1)")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+endif(${CMAKE_VERSION} VERSION_GREATER 3.1)
+
+# LOCAL PACKAGES
+find_package(SOAP REQUIRED)
+find_package(GTEST REQUIRED)
+
+# FIND PACKAGES
+find_package(PythonLibs)
+include_directories(${PYTHON_INCLUDE_DIRS})
+
+if(DEFINED ENV{BOOST_ROOT})
+    set(BOOST_ROOT "$ENV{BOOST_ROOT}")
+    message("-- BOOST_ROOT is set: ${BOOST_ROOT}")
+else(DEFINED ENV{BOOST_ROOT})
+    message("-- Note: BOOST_ROOT not set.")
+endif(DEFINED ENV{BOOST_ROOT})
+message("-- BOOST_ROOT is set: ${BOOST_ROOT}")
+find_package(Boost 1.60.0 COMPONENTS python mpi filesystem serialization)
+include_directories(${Boost_INCLUDE_DIRS})
+
+find_package(MPI REQUIRED)
+include_directories(${MPI_INCLUDE_PATH})
+
+
+# SUMMARIZE INCLUDES & LIBS
+get_property(local_dirs DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY INCLUDE_DIRECTORIES)
+get_property(local_libs DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY SOAPXX_LINK_LIBRARIES)
+
+message(STATUS "Include directories: ")
+foreach(dir ${local_dirs})
+  message(STATUS " o ${dir}")
+endforeach()
+message(STATUS "Linked libraries:    ")
+foreach(lib ${local_dirs})
+  message(STATUS " o ${lib}")
+endforeach()
+
+# SUBDIRECTORIES
+add_subdirectory(src)
+
diff --git a/cxx_tests/build.sh b/cxx_tests/build.sh
new file mode 100755
index 0000000..4e9c5c9
--- /dev/null
+++ b/cxx_tests/build.sh
@@ -0,0 +1,5 @@
+#! /bin/bash
+mkdir -p build
+cd build
+cmake .. && make && make install
+cd ..
diff --git a/cxx_tests/cmake/FindFFTW3.cmake b/cxx_tests/cmake/FindFFTW3.cmake
new file mode 100644
index 0000000..f471f40
--- /dev/null
+++ b/cxx_tests/cmake/FindFFTW3.cmake
@@ -0,0 +1,22 @@
+# - Find FFTW3
+# Find the native FFTW3 includes and library
+#
+# FFTW3_INCLUDES - where to find fftw3.h
+# FFTW3_LIBRARIES - List of libraries when using FFTW3.
+# FFTW3_FOUND - True if FFTW3 found.
+
+if (FFTW3_INCLUDES)
+  # Already in cache, be silent
+  set (FFTW3_FIND_QUIETLY TRUE)
+endif (FFTW3_INCLUDES)
+
+find_path (FFTW3_INCLUDES fftw3.h)
+
+find_library (FFTW3_LIBRARIES NAMES fftw3)
+
+# handle the QUIETLY and REQUIRED arguments and set FFTW3_FOUND to TRUE if
+# all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (FFTW3 DEFAULT_MSG FFTW3_LIBRARIES FFTW3_INCLUDES)
+
+mark_as_advanced (FFTW3_LIBRARIES FFTW3_INCLUDES)
diff --git a/cxx_tests/cmake/FindGSL.cmake b/cxx_tests/cmake/FindGSL.cmake
new file mode 100644
index 0000000..b1c3a1c
--- /dev/null
+++ b/cxx_tests/cmake/FindGSL.cmake
@@ -0,0 +1,59 @@
+# - Find gsl
+# Find the native GSL headers and libraries.
+#
+#  GSL_INCLUDE_DIRS - where to find gsl/gsl_linalg.h, etc.
+#  GSL_LIBRARIES    - List of libraries when using gsl.
+#  GSL_FOUND        - True if gsl found.
+#
+# Copyright 2009-2011 The VOTCA Development Team (http://www.votca.org)
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+find_package(PkgConfig)
+
+pkg_check_modules(PC_GSL gsl)
+find_path(GSL_INCLUDE_DIR gsl/gsl_linalg.h HINTS ${PC_GSL_INCLUDE_DIRS})
+find_library(GSL_LIBRARY NAMES gsl HINTS ${PC_GSL_LIBRARY_DIRS} )
+
+include(FindPackageHandleStandardArgs)
+# handle the QUIETLY and REQUIRED arguments and set GSL_FOUND to TRUE
+# if all listed variables are TRUE
+find_package_handle_standard_args(GSL DEFAULT_MSG GSL_LIBRARY GSL_INCLUDE_DIR)
+
+set(GSL_LIBRARIES "${GSL_LIBRARY}")
+set(GSL_INCLUDE_DIRS ${GSL_INCLUDE_DIR} )
+
+if (GSL_FOUND)
+  include(CheckLibraryExists)
+  #adding MATH_LIBRARIES here to allow static libs, this does not harm us as we are anyway using it
+  check_library_exists("${GSL_LIBRARIES};${MATH_LIBRARIES}" gsl_linalg_QR_decomp "" FOUND_QR_DECOMP)
+  if(NOT FOUND_QR_DECOMP)
+    #Some distributions like OpenSUSE need cblas in the linker flags, too
+    find_library(GSLCBLAS_LIBRARY NAMES gslcblas cblas HINTS ${PC_GSL_LIBRARY_DIRS} )
+    find_package_handle_standard_args(GSLCBLAS DEFAULT_MSG GSLCBLAS_LIBRARY)
+    if (GSLCBLAS_FOUND)
+      check_library_exists("${GSLCBLAS_LIBRARY};${MATH_LIBRARIES}" cblas_dsyrk "" FOUND_DSYRK)
+      if(NOT FOUND_DSYRK)
+        message(FATAL_ERROR "Could not find cblas_dsyrk in ${GSLCBLAS_LIBRARY}, take a look at the error message in ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log to find out what was going wrong. If you are using a static lib (.a) make sure you have specified all dependencies of libcblas in GSLCBLAS_LIBRARY by hand (i.e. -DGSLCBLAS_LIBRARY='/path/to/libcblas.so;/path/to/libm.so')! If your gsl was build against an different version of cblas, specify it in GSLCBLAS_LIBRARY")
+      endif(NOT FOUND_DSYRK)
+      set(GSL_LIBRARIES "${GSL_LIBRARY};${GSLCBLAS_LIBRARY}")
+      check_library_exists("${GSL_LIBRARIES};${MATH_LIBRARIES}" gsl_linalg_QR_lssolve "" FOUND_QR_DECOMP_AGAIN)
+    endif(GSLCBLAS_FOUND)
+  endif(NOT FOUND_QR_DECOMP)
+  if(NOT FOUND_QR_DECOMP AND NOT FOUND_QR_DECOMP_AGAIN)
+      message(FATAL_ERROR "Could not find gsl_linalg_QR_decompx in ${GSL_LIBRARY}, take a look at the error message in ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log to find out what was going wrong.  If you are using a static lib (.a) make sure you have specified all dependencies of libgsl in GSL_LIBRARY by hand (i.e. -DGSL_LIBRARY='/path/to/libgsl.so;/path/to/libm.so') !")
+  endif(NOT FOUND_QR_DECOMP AND NOT FOUND_QR_DECOMP_AGAIN)
+endif (GSL_FOUND)
+
+mark_as_advanced(GSL_INCLUDE_DIR GSL_LIBRARY)
diff --git a/cxx_tests/cmake/FindGTEST.cmake b/cxx_tests/cmake/FindGTEST.cmake
new file mode 100644
index 0000000..0821a5a
--- /dev/null
+++ b/cxx_tests/cmake/FindGTEST.cmake
@@ -0,0 +1,11 @@
+if(DEFINED ENV{GTEST_ROOT})
+    set(GTEST_ROOT "$ENV{GTEST_ROOT}")
+    message("-- GTEST_ROOT is set: ${GTEST_ROOT}")
+else()
+    message(FATAL_ERROR "-- Note: GTEST_ROOT not set.")
+endif()
+message("-- GTEST_ROOT is set: ${GTEST_ROOT}")
+set(GTEST_LIBRARIES "${GTEST_ROOT}/lib/libgtest_main.a" "${GTEST_ROOT}/lib/libgtest.a")
+set(GTEST_INCLUDE_DIRS "${GTEST_ROOT}/include")
+include_directories(${GTEST_INCLUDE_DIRS})
+
diff --git a/cxx_tests/cmake/FindMKL.cmake b/cxx_tests/cmake/FindMKL.cmake
new file mode 100644
index 0000000..b35d9d6
--- /dev/null
+++ b/cxx_tests/cmake/FindMKL.cmake
@@ -0,0 +1,339 @@
+# CMake script to detect Intel(R) Math Kernel Library (MKL)
+#
+# This will try to find Intel MKL libraries, and include path by automatic
+# search through typical install locations and if failed it will
+# examine MKLROOT environment variable.
+# Note, MKLROOT is not set by IPP installer, it should be set manually.
+#
+# Usage example:
+#   set(MKL_USE_STATIC_LIBS ON)
+#   find_package(MKL)
+#   if (MKL_FOUND)
+#      include_directories(${MKL_INCLUDE_DIRS})
+#      add_executable(foo foo.cc)
+#      target_link_libraries(foo ${MKL_LIBRARIES})
+#   endif()
+#
+# Variables used by this module, they can change the default behaviour and
+# need to be set before calling find_package:
+#
+#   MKL_ADDITIONAL_VERSIONS      A list of version numbers to use for searching
+#                                the MKL include directory.
+#
+#   MKL_USE_STATIC_LIBS          Can be set to ON to force the use of the static
+#                                boost libraries. Defaults to OFF.
+#
+#   MKL_FIND_DEBUG               Set this to TRUE to enable debugging output
+#                                of FindMKL.cmake if you are having problems.
+#
+# On return this will define:
+#   MKL_FOUND                   Indicates whether MKL was found (True/False)
+#   MKL_INCLUDE_DIRS            MKL include folder
+#   MKL_LIBRARY_DIRS            MKL libraries folder
+#   MKL_LIBRARIES               MKL libraries names
+#
+# NOTE: this script has only been tested with Intel(R) Parallel Studio XE 2011
+# and may need changes for compatibility with older versions.
+#
+# Adapted from OpenCV IPP detection script
+#   https://code.ros.org/trac/opencv/browser/trunk/opencv/OpenCVFindIPP.cmake
+# Many portions taken from FindBoost.cmake
+
+# TODO:
+# - caller needs to link with libiomp5md.lib or /Qopenmp...
+# - runtime DLLs:
+#   <Composer XE directory> -> C:\Program Files\Intel\ComposerXE-2011
+#     redist\ia32\mkl
+#     redist\intel64\mkl
+
+set(_MKL_IA32 FALSE)
+set(_MKL_INTEL64 FALSE)
+if (CMAKE_SIZEOF_VOID_P EQUAL 4)
+    set(_MKL_IA32 TRUE)
+elseif (CMAKE_SIZEOF_VOID_P EQUAL 8)
+    set(_MKL_INTEL64 TRUE)
+else()
+    message(FATAL_ERROR "Unsupported 'void *' size (${SIZEOF_VOID_P})")
+endif()
+
+# Versions should be listed is decreasing order of preference
+set(_MKL_TEST_VERSIONS ${MKL_ADDITIONAL_VERSIONS}
+    "2011"
+    # alternative form: "2011.xxx.y"
+    # (y is the release-update number and xxx is the package number)
+)
+
+set(MKL_FIND_QUIETLY true)
+if (MKL_FIND_VERSION AND NOT MKL_FIND_QUIETLY)
+    message(WARNING "Requesting a specific version of Intel(R) MKL is not supported")
+endif()
+
+# Use environment variables from Intel build scripts, if available
+if (NOT MKL_ROOT AND NOT $ENV{MKLROOT} STREQUAL "")
+  set(MKL_ROOT $ENV{MKLROOT})
+endif()
+
+if (MKL_ROOT)
+  file(TO_CMAKE_PATH ${MKL_ROOT} MKL_ROOT)
+endif()
+
+if (NOT INTEL_ROOT AND NOT $ENV{INTELROOT} STREQUAL "")
+  set(INTEL_ROOT $ENV{INTELROOT})
+endif()
+
+if (INTEL_ROOT)
+  file(TO_CMAKE_PATH ${INTEL_ROOT} INTEL_ROOT)
+endif()
+
+if (MKL_FIND_DEBUG)
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "_MKL_TEST_VERSIONS = ${_MKL_TEST_VERSIONS}")
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "MKL_ADDITIONAL_VERSIONS = ${MKL_ADDITIONAL_VERSIONS}")
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "MKL_USE_STATIC_LIBS = ${MKL_USE_STATIC_LIBS}")
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "MKL_ROOT = ${MKL_ROOT}")
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "INTEL_ROOT = ${INTEL_ROOT}")
+endif()
+
+# Find MKL include directory
+
+set(_MKL_ROOT_SEARCH_DIRS
+  ${MKL_ROOT}
+)
+
+foreach (_MKL_VER ${_MKL_TEST_VERSIONS})
+    if (WIN32)
+        list(APPEND _MKL_ROOT_SEARCH_DIRS "$ENV{ProgramFiles}/Intel/Composer XE/mkl")
+    else()
+        list(APPEND _MKL_ROOT_SEARCH_DIRS "/opt/intel/composerxe-${_MKL_VER}/mkl")
+    endif()
+endforeach()
+
+if (MKL_FIND_DEBUG)
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "_MKL_ROOT_SEARCH_DIRS = ${_MKL_ROOT_SEARCH_DIRS}")
+endif()
+
+find_path(MKL_INCLUDE_DIR
+    NAMES mkl.h
+    PATHS ${_MKL_ROOT_SEARCH_DIRS}
+    PATH_SUFFIXES include
+    DOC "The path to Intel(R) MKL header files"
+)
+
+if (MKL_INCLUDE_DIR)
+    if (MKL_FIND_DEBUG)
+        message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                       "location of mkl.h: ${MKL_INCLUDE_DIR}/mkl.h")
+    endif()
+else()
+    if (MKL_FIND_DEBUG)
+        message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                       "unable to find Intel(R) MKL header files. Please set MKLROOT"
+                       " to the root directory containing MKL.")
+    endif()
+endif()
+
+# Find MKL library directory
+
+set(_INTEL_LIBRARY_DIR_SUFFIXES "lib")
+if (_MKL_IA32)
+    list(APPEND _INTEL_LIBRARY_DIR_SUFFIXES "lib/ia32")
+elseif (_MKL_INTEL64)
+    list(APPEND _INTEL_LIBRARY_DIR_SUFFIXES "lib/intel64")
+else()
+    message(FATAL_ERROR "unreachable")
+endif()
+
+set(_MKL_LIBRARY_SEARCH_DIRS ${_MKL_ROOT_SEARCH_DIRS})
+if (MKL_INCLUDE_DIR)
+    list(APPEND _MKL_LIBRARY_SEARCH_DIRS "${MKL_INCLUDE_DIR}/..")
+endif()
+
+if (MKL_FIND_DEBUG)
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "_INTEL_LIBRARY_DIR_SUFFIXES = ${_INTEL_LIBRARY_DIR_SUFFIXES}")
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "_MKL_LIBRARY_SEARCH_DIRS = ${_MKL_LIBRARY_SEARCH_DIRS}")
+endif()
+
+set(MKL_LIB_PREFIX "mkl_")
+if (MKL_USE_STATIC_LIBS)
+    if (_MKL_IA32)
+        if (WIN32)
+            set(_MKL_LIBRARIES intel_c)
+        else()
+            set(_MKL_LIBRARIES intel)
+        endif()
+    elseif (_MKL_INTEL64)
+        set(_MKL_LIBRARIES intel_lp64)
+    else()
+        message(FATAL_ERROR "unreachable")
+    endif()
+
+    list(APPEND _MKL_LIBRARIES intel_thread)
+    list(APPEND _MKL_LIBRARIES core)
+else()
+    set(_MKL_LIBRARIES rt)
+endif()
+
+set(_MKL_MISSING_LIBRARIES "")
+set(MKL_LIBRARIES "")
+set(MKL_LIBRARY_DIRS "")
+# Find MKL libraries
+foreach (_MKL_LIB_RAW ${_MKL_LIBRARIES})
+    set(_MKL_LIB ${MKL_LIB_PREFIX}${_MKL_LIB_RAW})
+    string(TOUPPER ${_MKL_LIB} _MKL_LIB_UPPER)
+
+    find_library(${_MKL_LIB_UPPER}_LIBRARY
+        NAMES ${_MKL_LIB}
+        PATHS ${_MKL_LIBRARY_SEARCH_DIRS}
+        PATH_SUFFIXES ${_INTEL_LIBRARY_DIR_SUFFIXES}
+        DOC "The path to Intel(R) MKL ${_MKL_LIB_RAW} library"
+    )
+    mark_as_advanced(${_MKL_LIB_UPPER}_LIBRARY)
+
+    if (NOT ${_MKL_LIB_UPPER}_LIBRARY)
+        list(APPEND _MKL_MISSING_LIBRARIES ${_MKL_LIB})
+    else()
+        list(APPEND MKL_LIBRARIES ${${_MKL_LIB_UPPER}_LIBRARY})
+        if (MKL_FIND_DEBUG)
+            message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                           "Found ${_MKL_LIB}: ${${_MKL_LIB_UPPER}_LIBRARY}")
+        endif()
+
+        get_filename_component(_MKL_LIB_PATH "${${_MKL_LIB_UPPER}_LIBRARY}" PATH)
+        list(APPEND MKL_LIBRARY_DIRS ${_MKL_LIB_PATH})
+    endif()
+endforeach()
+
+## Find OpenMP, pthread and math libraries
+
+set(_INTEL_LIBRARY_SEARCH_DIRS
+  ${INTEL_ROOT}
+  ${INTEL_ROOT}/compiler
+)
+
+foreach(_MKL_DIR ${_MKL_ROOT_SEARCH_DIRS})
+    list(APPEND _INTEL_LIBRARY_SEARCH_DIRS "${_MKL_DIR}/..")
+    list(APPEND _INTEL_LIBRARY_SEARCH_DIRS "${_MKL_DIR}/../compiler")
+endforeach()
+
+if (MKL_FIND_DEBUG)
+    message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                   "_INTEL_LIBRARY_SEARCH_DIRS = ${_INTEL_LIBRARY_SEARCH_DIRS}")
+endif()
+
+if (NOT WIN32)
+    find_library(PTHREAD_LIBRARY pthread DOC "Path to POSIX threads library")
+endif()
+
+set(_IOMP5_LIB iomp5)
+if (WIN32)
+  if (MKL_USE_STATIC_LIBS)
+      list(APPEND _IOMP5_LIB libiomp5mt.lib)
+  else()
+      list(APPEND _IOMP5_LIB libiomp5md.lib)
+  endif()
+endif()
+
+find_library(IOMP5_LIBRARY
+    NAMES ${_IOMP5_LIB}
+    PATHS ${_INTEL_LIBRARY_SEARCH_DIRS}
+    PATH_SUFFIXES ${_INTEL_LIBRARY_DIR_SUFFIXES}
+    DOC "Path to OpenMP runtime library"
+)
+
+if (NOT IOMP5_LIBRARY)
+    # we could instead fallback to default library (via FindOpenMP.cmake)
+    list(APPEND _MKL_MISSING_LIBRARIES IOMP5)
+else()
+    list(APPEND MKL_LIBRARIES ${IOMP5_LIBRARY})
+    if (MKL_FIND_DEBUG)
+        message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                       "Found IOMP5_LIBRARY: ${IOMP5_LIBRARY}")
+    endif()
+
+    get_filename_component(_MKL_LIB_PATH "${IOMP5_LIBRARY}" PATH)
+    list(APPEND MKL_LIBRARY_DIRS ${_MKL_LIB_PATH})
+endif()
+
+# Optimized math library (optional)
+set(_MATH_LIB imf)  # linked by default with Intel compiler
+if (WIN32)
+  if (MKL_USE_STATIC_LIBS)
+    list(APPEND _MATH_LIB libmmds.lib)  # assumes (/MD) otherwise libmmt.lib (for /MT)
+  else()
+      list(APPEND _MATH_LIB libmmd.lib)
+  endif()
+endif()
+
+find_library(MATH_LIBRARY
+    NAMES ${_MATH_LIB}
+    PATHS ${_INTEL_LIBRARY_SEARCH_DIRS}
+    PATH_SUFFIXES ${_INTEL_LIBRARY_DIR_SUFFIXES}
+    DOC "Path to optimized math library"
+)
+
+if (NOT MATH_LIBRARY)
+    # we could instead fallback to default library (via FindOpenMP.cmake)
+    list(APPEND _MKL_MISSING_LIBRARIES MATH)
+else()
+    list(APPEND MKL_LIBRARIES ${MATH_LIBRARY})
+    if (MKL_FIND_DEBUG)
+        message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                       "Found MATH_LIBRARY: ${MATH_LIBRARY}")
+    endif()
+
+    get_filename_component(_MKL_LIB_PATH "${MATH_LIBRARY}" PATH)
+    list(APPEND MKL_LIBRARY_DIRS ${_MKL_LIB_PATH})
+endif()
+
+# Check all required libraries are available
+list(REMOVE_DUPLICATES MKL_LIBRARY_DIRS)
+
+set(MKL_INCLUDE_DIRS
+    ${MKL_INCLUDE_DIR}
+)
+
+set(MKL_FOUND TRUE)
+if (NOT MKL_INCLUDE_DIR)
+    set(MKL_FOUND FALSE)
+    if (MKL_FIND_DEBUG)
+        message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                       "MKL not found - MKL_INCLUDE_DIR was empty")
+    endif()
+elseif (_MKL_MISSING_LIBRARIES)
+    set(MKL_FOUND FALSE)
+    if (MKL_FIND_DEBUG)
+        message(STATUS "[ ${CMAKE_CURRENT_LIST_FILE}:${CMAKE_CURRENT_LIST_LINE} ] "
+                       "MKL not found - the following libraries are missing: "
+                       "${_MKL_MISSING_LIBRARIES}")
+    endif()
+endif()
+
+if (MKL_FOUND)
+    if (NOT MKL_FIND_QUIETLY OR MKL_FIND_DEBUG)
+        message(STATUS
+            "Intel(R) MKL was found:\n"
+            "  MKL_INCLUDE_DIRS: ${MKL_INCLUDE_DIRS}\n"
+            "  MKL_LIBRARY_DIRS: ${MKL_LIBRARY_DIRS}\n"
+            "  MKL_LIBRARIES: ${MKL_LIBRARIES}"
+        )
+    endif()
+else()
+    if (MKL_FIND_REQUIRED)
+        message(SEND_ERROR "Intel(R) MKL could not be found.")
+    else()
+        message(STATUS "Intel(R) MKL could not be found.")
+    endif()
+endif()
+
+mark_as_advanced(
+    MKL_INCLUDE_DIR
+    MKL_INCLUDE_DIRS
+    MKL_LIBRARY_DIRS
+)
diff --git a/cxx_tests/cmake/FindSOAP.cmake b/cxx_tests/cmake/FindSOAP.cmake
new file mode 100644
index 0000000..4331a0d
--- /dev/null
+++ b/cxx_tests/cmake/FindSOAP.cmake
@@ -0,0 +1,11 @@
+if(DEFINED ENV{SOAP_ROOT})
+    set(SOAP_ROOT "$ENV{SOAP_ROOT}")
+    message("-- SOAP_ROOT is set: ${SOAP_ROOT}")
+else()
+    message(FATAL_ERROR "-- Note: SOAP_ROOT not set.")
+endif()
+message("-- SOAP_ROOT is set: ${SOAP_ROOT}")
+set(SOAP_LIBRARIES "${SOAP_ROOT}/soap/_soapxx.so")
+set(SOAP_INCLUDE_DIRS "${SOAP_ROOT}/soap/include")
+include_directories(${SOAP_INCLUDE_DIRS})
+
diff --git a/cxx_tests/gtest.sh b/cxx_tests/gtest.sh
new file mode 100755
index 0000000..9939dc5
--- /dev/null
+++ b/cxx_tests/gtest.sh
@@ -0,0 +1,4 @@
+#! /bin/bash
+
+soap_tests/test.exe
+
diff --git a/cxx_tests/src/CMakeLists.txt b/cxx_tests/src/CMakeLists.txt
new file mode 100644
index 0000000..9be1e25
--- /dev/null
+++ b/cxx_tests/src/CMakeLists.txt
@@ -0,0 +1,17 @@
+# COMPILE LIBRARIES
+set(LD_LIBRARIES ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${MPI_LIBRARIES} ${SOAP_LIBRARIES} ${GTEST_LIBRARIES})
+
+add_executable(test_functions.exe test_functions.cpp)
+target_link_libraries(test_functions.exe ${LD_LIBRARIES})
+install(TARGETS test_functions.exe DESTINATION ${LOCAL_INSTALL_DIR})
+
+add_executable(test_dylm.exe test_dylm.cpp)
+target_link_libraries(test_dylm.exe ${LD_LIBRARIES})
+install(TARGETS test_dylm.exe DESTINATION ${LOCAL_INSTALL_DIR})
+
+file(GLOB local_sources gtest_*.cpp)
+
+add_executable(test.exe ${local_sources})
+target_link_libraries(test.exe ${LD_LIBRARIES})
+install(TARGETS test.exe DESTINATION ${LOCAL_INSTALL_DIR})
+
diff --git a/cxx_tests/src/gtest_dylm.cpp b/cxx_tests/src/gtest_dylm.cpp
new file mode 100644
index 0000000..fbf82f7
--- /dev/null
+++ b/cxx_tests/src/gtest_dylm.cpp
@@ -0,0 +1,190 @@
+#include <iostream>
+#include <vector>
+#include <boost/format.hpp>
+#include <gtest/gtest.h>
+#include <soap/functions.hpp>
+
+TEST(TestFunctions, PowerExponentZero) {
+    // CAREFUL, THIS CAN LEAD TO NAN'S DEPENDING ON IMPLEMENTATION!
+    double p = pow(0.,0);
+    std::complex<double> q = soap::pow_nnan(std::complex<double>(0.,0.),0);
+    
+    EXPECT_EQ(p, 1);
+    EXPECT_EQ(q, std::complex<double>(1.,0.));
+}
+
+TEST(TestFunctions, GradientYlm) {
+
+    // SETUP
+    std::vector< soap::vec > r_list;    
+    r_list.push_back(soap::vec(2.,0.,0.));
+    r_list.push_back(soap::vec(0.,1.,0.));
+    r_list.push_back(soap::vec(0.,0.,1.));
+    r_list.push_back(soap::vec(std::sqrt(0.5),0.,std::sqrt(0.5)));
+    r_list.push_back(soap::vec(-0.2,0.3,0.7));
+    r_list.push_back(soap::vec(-0.2,-0.4,-0.1));
+
+    std::vector< std::pair<int,int> > lm_list;
+    lm_list.push_back(std::pair<int,int>(1,0));
+    lm_list.push_back(std::pair<int,int>(1,1));
+    lm_list.push_back(std::pair<int,int>(2,0));
+    lm_list.push_back(std::pair<int,int>(2,1));
+    lm_list.push_back(std::pair<int,int>(2,-1));
+    lm_list.push_back(std::pair<int,int>(2,2));
+    lm_list.push_back(std::pair<int,int>(2,-2));
+
+    std::vector<std::complex<double> > results;
+    results.push_back(std::complex<double>(-1.4959138e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+2.4430126e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-1.8319660e-33, +0.0000000e+00));
+    results.push_back(std::complex<double>(-2.9918275e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(+4.8860251e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(-2.4430126e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+2.4430126e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+1.4011873e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-2.1017810e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+1.3011025e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-1.0154458e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-2.0308916e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+1.0154458e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(-6.4769779e-34, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, -1.7274707e-01));
+    results.push_back(std::complex<double>(+1.0577708e-17, -0.0000000e+00));
+    results.push_back(std::complex<double>(-3.4549415e-01, +2.1154717e-17));
+    results.push_back(std::complex<double>(+2.1155415e-17, -2.5907484e-33));
+    results.push_back(std::complex<double>(+1.2953528e-33, +2.1155415e-17));
+    results.push_back(std::complex<double>(-3.4549415e-01, +6.9868116e-22));
+    results.push_back(std::complex<double>(-6.9868116e-22, -3.4549415e-01));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(-1.7274707e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, -3.4549415e-01));
+    results.push_back(std::complex<double>(+1.7274707e-01, -0.0000000e+00));
+    results.push_back(std::complex<double>(-4.1046975e-01, -4.2462388e-02));
+    results.push_back(std::complex<double>(-4.2462388e-02, -3.7508443e-01));
+    results.push_back(std::complex<double>(-9.9078905e-02, +1.4861836e-01));
+    results.push_back(std::complex<double>(-6.1032432e-01, +2.8721145e-01));
+    results.push_back(std::complex<double>(+2.8721145e-01, -1.7950715e-01));
+    results.push_back(std::complex<double>(+7.1802861e-02, +1.4360572e-01));
+    results.push_back(std::complex<double>(-3.5475869e-33, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+5.7936491e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(-4.3445409e-49, +0.0000000e+00));
+    results.push_back(std::complex<double>(-7.0951738e-33, +0.0000000e+00));
+    results.push_back(std::complex<double>(+1.1587298e-16, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(-6.6904654e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+6.6904654e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+4.8244079e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-7.2366119e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+4.4798074e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+8.5820834e-02, +0.0000000e+00));
+    results.push_back(std::complex<double>(+1.7164167e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-8.5820834e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+2.3652473e-17, -0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, -2.3652473e-17));
+    results.push_back(std::complex<double>(-3.8627420e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-4.7304947e-17, +5.7930895e-33));
+    results.push_back(std::complex<double>(+5.7930895e-33, +4.7304947e-17));
+    results.push_back(std::complex<double>(-4.7303384e-17, -7.7254840e-01));
+    results.push_back(std::complex<double>(-7.7254840e-01, +1.5622986e-21));
+    results.push_back(std::complex<double>(-1.5622986e-21, -7.7254840e-01));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(-3.3449648e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(-0.0000000e+00, -5.4627422e-01));
+    results.push_back(std::complex<double>(+3.3449648e-17, -0.0000000e+00));
+    results.push_back(std::complex<double>(-7.5968600e-01, -1.6881911e-01));
+    results.push_back(std::complex<double>(-1.6881911e-01, -6.1900340e-01));
+    results.push_back(std::complex<double>(-1.4470209e-01, +2.1705314e-01));
+    results.push_back(std::complex<double>(+2.2773536e-01, -2.8028967e-01));
+    results.push_back(std::complex<double>(-2.8028967e-01, -1.9269915e-01));
+    results.push_back(std::complex<double>(+6.6568797e-01, +1.3313759e+00));
+    results.push_back(std::complex<double>(-2.3652473e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, -2.3652473e-17));
+    results.push_back(std::complex<double>(+3.8627420e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+4.7304947e-17, +5.7930895e-33));
+    results.push_back(std::complex<double>(-5.7930895e-33, +4.7304947e-17));
+    results.push_back(std::complex<double>(+4.7303384e-17, -7.7254840e-01));
+    results.push_back(std::complex<double>(+7.7254840e-01, +1.5622986e-21));
+    results.push_back(std::complex<double>(+1.5622986e-21, -7.7254840e-01));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+3.3449648e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, -5.4627422e-01));
+    results.push_back(std::complex<double>(-3.3449648e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(+7.5968600e-01, -1.6881911e-01));
+    results.push_back(std::complex<double>(+1.6881911e-01, -6.1900340e-01));
+    results.push_back(std::complex<double>(+1.4470209e-01, +2.1705314e-01));
+    results.push_back(std::complex<double>(-2.2773536e-01, -2.8028967e-01));
+    results.push_back(std::complex<double>(+2.8028967e-01, -1.9269915e-01));
+    results.push_back(std::complex<double>(-6.6568797e-01, +1.3313759e+00));
+    results.push_back(std::complex<double>(+1.4482963e-33, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +3.8627420e-01));
+    results.push_back(std::complex<double>(-2.3652473e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(+9.4606768e-17, +7.7254840e-01));
+    results.push_back(std::complex<double>(-8.6895864e-33, -4.7304947e-17));
+    results.push_back(std::complex<double>(+4.7304947e-17, -5.7929938e-33));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +0.0000000e+00));
+    results.push_back(std::complex<double>(+2.7313711e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, +5.4627422e-01));
+    results.push_back(std::complex<double>(-2.7313711e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-2.6930668e-01, +3.2557971e-01));
+    results.push_back(std::complex<double>(-3.4366747e-01, -1.7685812e-01));
+    results.push_back(std::complex<double>(+7.0341296e-02, +1.6881911e-01));
+    results.push_back(std::complex<double>(-1.1561949e+00, -9.1094143e-01));
+    results.push_back(std::complex<double>(+6.3065176e-01, +3.8539830e-01));
+    results.push_back(std::complex<double>(-2.1021725e-01, +2.8028967e-01));
+    results.push_back(std::complex<double>(+1.4482963e-33, -0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, -3.8627420e-01));
+    results.push_back(std::complex<double>(-2.3652473e-17, +0.0000000e+00));
+    results.push_back(std::complex<double>(+9.4606768e-17, -7.7254840e-01));
+    results.push_back(std::complex<double>(-8.6895864e-33, +4.7304947e-17));
+    results.push_back(std::complex<double>(+4.7304947e-17, +5.7929938e-33));
+    results.push_back(std::complex<double>(+0.0000000e+00, -0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, -0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, -0.0000000e+00));
+    results.push_back(std::complex<double>(+2.7313711e-01, -0.0000000e+00));
+    results.push_back(std::complex<double>(+0.0000000e+00, -5.4627422e-01));
+    results.push_back(std::complex<double>(-2.7313711e-01, +0.0000000e+00));
+    results.push_back(std::complex<double>(-2.6930668e-01, -3.2557971e-01));
+    results.push_back(std::complex<double>(-3.4366747e-01, +1.7685812e-01));
+    results.push_back(std::complex<double>(+7.0341296e-02, -1.6881911e-01));
+    results.push_back(std::complex<double>(-1.1561949e+00, +9.1094143e-01));
+    results.push_back(std::complex<double>(+6.3065176e-01, -3.8539830e-01));
+    results.push_back(std::complex<double>(-2.1021725e-01, -2.8028967e-01));   
+    
+    // COMPUTE
+    int res_idx = -1;
+    for (int lm = 0; lm < lm_list.size(); ++lm) {
+        int l = lm_list[lm].first;
+        int m = lm_list[lm].second;
+
+        for (int n = 0; n < r_list.size(); ++n) {
+            soap::vec r = r_list[n];    
+            std::vector<std::complex<double> > dylm = soap::GradSphericalYlm::eval(l, m, r);
+
+            // VERIFY
+            res_idx += 1;
+            EXPECT_NEAR(dylm[0].real(), results[res_idx].real(), 1e-7);
+            EXPECT_NEAR(dylm[0].imag(), results[res_idx].imag(), 1e-7);
+            res_idx += 1;
+            EXPECT_NEAR(dylm[1].real(), results[res_idx].real(), 1e-7);
+            EXPECT_NEAR(dylm[1].imag(), results[res_idx].imag(), 1e-7);            
+            res_idx += 1;
+            EXPECT_NEAR(dylm[2].real(), results[res_idx].real(), 1e-7);
+            EXPECT_NEAR(dylm[2].imag(), results[res_idx].imag(), 1e-7);
+        }        
+    }
+}
+
+
+
+
diff --git a/cxx_tests/src/test_dylm.cpp b/cxx_tests/src/test_dylm.cpp
new file mode 100644
index 0000000..37ff96b
--- /dev/null
+++ b/cxx_tests/src/test_dylm.cpp
@@ -0,0 +1,52 @@
+#include <iostream>
+#include <vector>
+#include <boost/format.hpp>
+#include <soap/functions.hpp>
+
+int main() {
+
+    std::cout << "soapxx/test/functions" << std::endl;
+
+    std::vector< soap::vec > r_list;    
+    r_list.push_back(soap::vec(2.,0.,0.));
+    r_list.push_back(soap::vec(0.,1.,0.));
+    r_list.push_back(soap::vec(0.,0.,1.));
+    r_list.push_back(soap::vec(std::sqrt(0.5),0.,std::sqrt(0.5)));
+    r_list.push_back(soap::vec(-0.2,0.3,0.7));
+
+    std::vector< std::pair<int,int> > lm_list;
+    lm_list.push_back(std::pair<int,int>(1,0));
+    lm_list.push_back(std::pair<int,int>(1,-1));
+    lm_list.push_back(std::pair<int,int>(1,1));
+    lm_list.push_back(std::pair<int,int>(2,0));
+    lm_list.push_back(std::pair<int,int>(2,1));
+    lm_list.push_back(std::pair<int,int>(2,2));
+    
+    // CAREFUL, THIS CAN LEAD TO NAN'S DEPENDING ON IMPLEMENTATION!
+    double p = pow(0.,0);
+    std::cout << p << std::endl;
+    
+    std::complex<double> q = soap::pow_nnan(std::complex<double>(0.,0.),0);
+    std::cout << q << std::endl;
+
+    for (int lm = 0; lm < lm_list.size(); ++lm) {
+        int l = lm_list[lm].first;
+        int m = lm_list[lm].second;
+
+        std::cout << "====" << l << m << "====" << std::endl;
+
+        for (int n = 0; n < r_list.size(); ++n) {
+            soap::vec r = r_list[n];    
+            
+            std::vector<std::complex<double> > dylm = soap::GradSphericalYlm::eval(l, m, r);
+            
+            std::cout << "r = " << r << std::flush;
+            std::cout << "=>" << std::flush;
+            std::cout << dylm[0] << dylm[1] << dylm[2] << std::endl;        
+        }        
+    }
+}
+
+
+
+
diff --git a/cxx_tests/src/test_dylm.py b/cxx_tests/src/test_dylm.py
new file mode 100644
index 0000000..a8d8075
--- /dev/null
+++ b/cxx_tests/src/test_dylm.py
@@ -0,0 +1,116 @@
+import scipy.special
+import numpy as np
+
+
+def theta_hat(theta, phi):
+    x = np.cos(theta)*np.cos(phi)
+    y = np.cos(theta)*np.sin(phi)
+    z = -np.sin(theta)
+    return np.array([x,y,z])
+    
+def phi_hat(theta, phi):
+    x = -np.sin(phi)
+    y = +np.cos(phi)
+    z = 0.
+    return np.array([x,y,z])
+
+def R_theta_phi(r):
+    r = np.array(r)
+    R = np.dot(r,r)**0.5
+    theta = np.arccos(r[2]/R)
+    
+    if r[0] == 0.:
+        if r[1] <= 0.:
+            phi = -0.5*np.pi
+        else:
+            phi = +0.5*np.pi
+    else:
+        phi = np.arctan2(r[1],r[0])
+    return R, theta, phi
+    
+
+def psi_10(r):
+    R, theta, phi = R_theta_phi(r)    
+    pre = - (3./(4.*np.pi))**0.5 * np.sin(theta)    
+    direction = theta_hat(theta, phi)
+    return pre*direction
+
+def psi_11(r):
+    R, theta, phi = R_theta_phi(r)    
+    pre = - (3./(8.*np.pi))**0.5 * np.exp(phi*1.j)
+    direction = np.cos(theta)*theta_hat(theta, phi) + 1.j*phi_hat(theta, phi)
+    return pre*direction
+    
+def psi_20(r):
+    R, theta, phi = R_theta_phi(r)
+    return -1.5*(5./np.pi)**0.5 * np.sin(theta)*np.cos(theta)*theta_hat(theta, phi)
+
+def psi_21(r):
+    R, theta, phi = R_theta_phi(r)
+    return -1.*(15./8./np.pi)**0.5*np.exp(phi*1.j)*( np.cos(2*theta)*theta_hat(theta,phi) + 1.j*np.cos(theta)*phi_hat(theta,phi) )
+    
+def psi_2_1(r):
+    dylm_21 = psi_21(r)
+    return -1*np.conjugate(dylm_21)
+    
+def psi_22(r):
+    R, theta, phi = R_theta_phi(r)
+    return (15./8./np.pi)**0.5*np.sin(theta)*np.exp(2.j*phi)*( np.cos(theta)*theta_hat(theta,phi) + 1.j*phi_hat(theta, phi) )
+    
+def psi_2_2(r):
+    dylm_22 = psi_22(r)
+    return np.conjugate(dylm_22)
+    
+r_list = []
+r_list.append([2.,0.,0.])
+r_list.append([0.,1.,0.])
+r_list.append([0.,0.,1.])
+r_list.append([0.5**0.5,0.,0.5**0.5])
+r_list.append([-0.2,0.3,0.7])
+r_list.append([-0.2,-0.4,-0.1])
+
+
+results = []
+
+print "==== THETA, PHI ===="
+for r in r_list:
+    R, theta, phi = R_theta_phi(r)
+    print r, "=>", theta, phi
+
+print "==== 10 ===="
+for r in r_list:
+    print r, "=>", " => unnorm", psi_10(r)/np.dot(r,r)**0.5
+    results.append(psi_10(r)/np.dot(r,r)**0.5) 
+    
+print "==== 11 ===="  
+for r in r_list:
+    print r, "=>", " => unnorm", psi_11(r)/np.dot(r,r)**0.5
+    results.append(psi_11(r)/np.dot(r,r)**0.5)
+    
+print "==== 20 ===="  
+for r in r_list:
+    print r, "=>", " => unnorm", psi_20(r)/np.dot(r,r)**0.5
+    results.append(psi_20(r)/np.dot(r,r)**0.5)
+    
+print "==== 21 ===="  
+for r in r_list:
+    print r, "=>", " => unnorm", psi_21(r)/np.dot(r,r)**0.5
+    results.append(psi_21(r)/np.dot(r,r)**0.5)
+    
+for r in r_list:
+    results.append(psi_2_1(r)/np.dot(r,r)**0.5)
+
+print "==== 22 ===="  
+for r in r_list:
+    print r, "=>", " => unnorm", psi_22(r)/np.dot(r,r)**0.5
+    results.append(psi_22(r)/np.dot(r,r)**0.5)
+    
+for r in r_list:
+    results.append(psi_2_2(r)/np.dot(r,r)**0.5)
+    
+for r in results:
+    print "results.push_back(std::complex<double>(%+1.7e, %+1.7e));" % (r[0].real, r[0].imag)
+    print "results.push_back(std::complex<double>(%+1.7e, %+1.7e));" % (r[1].real, r[1].imag)
+    print "results.push_back(std::complex<double>(%+1.7e, %+1.7e));" % (r[2].real, r[2].imag)
+    
+
diff --git a/cxx_tests/src/test_functions.cpp b/cxx_tests/src/test_functions.cpp
new file mode 100644
index 0000000..e672e3a
--- /dev/null
+++ b/cxx_tests/src/test_functions.cpp
@@ -0,0 +1,26 @@
+#include <iostream>
+#include <vector>
+#include <soap/functions.hpp>
+#include <boost/format.hpp>
+
+int main() {
+
+    std::cout << "soapxx/test/functions" << std::endl;
+
+    int N = 10;
+    std::vector<double> r_list;
+    r_list.push_back(0.);
+    r_list.push_back(0.1);
+    r_list.push_back(1.);
+
+    soap::ModifiedSphericalBessel1stKind sph_in(N);
+
+    for (int i = 0; i < r_list.size(); ++i) {
+        sph_in.evaluate(r_list[i], true); 
+        std::cout << "r = " << r_list[i] << std::endl;    
+        for (int n = 0; n <= 10; ++n) {
+            std::cout << boost::format("%1$2d %2$+1.7e %3$+1.7e") % n % sph_in._in[n] % sph_in._din[n] << std::endl;
+        }
+    }
+
+}
diff --git a/cxx_tests/src/test_functions.py b/cxx_tests/src/test_functions.py
new file mode 100644
index 0000000..0b95163
--- /dev/null
+++ b/cxx_tests/src/test_functions.py
@@ -0,0 +1,10 @@
+import scipy.special
+
+n = 10
+
+for r in [0.,0.1,1.]:
+    sph_in, sph_din = scipy.special.sph_in(n, r)
+    print "r =", r
+
+    for i in range(n):
+        print "%2d %+1.7e %+1.7e" % (i, sph_in[i], sph_din[i])
diff --git a/src/functions.cpp b/src/functions.cpp
deleted file mode 100644
index f0586bd..0000000
--- a/src/functions.cpp
+++ /dev/null
@@ -1,95 +0,0 @@
-#include "functions.hpp"
-
-#include <boost/math/special_functions/erf.hpp>
-
-namespace soap {
-
-// ======================
-// RadialGaussian
-// ======================
-
-RadialGaussian::RadialGaussian(double r0, double sigma)
-: _r0(r0),
-  _sigma(sigma),
-  _alpha(1./(2.*sigma*sigma)) {
-
-	// COMPUTE NORMALIZATION S g^2 r^2 dr
-	// This normalization is to be used for radial basis functions
-    double w = 2*_alpha;
-    double W0 = 2*_alpha*_r0;
-    _integral_r2_g2_dr =
-        1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
-            2*sqrt(w)*W0 +
-			sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
-			    1 - boost::math::erf<double>(-W0/sqrt(w))
-			)
-        );
-    _norm_r2_g2_dr = 1./sqrt(_integral_r2_g2_dr);
-
-    // COMPUTE NORMALIZATION S g r^2 dr
-    // This normalization is to be used for "standard" radial Gaussians
-    w = _alpha;
-	W0 = _alpha*_r0;
-	_integral_r2_g_dr =
-		1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
-			2*sqrt(w)*W0 +
-			sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
-				1 - boost::math::erf<double>(-W0/sqrt(w))
-			)
-		);
-	_norm_r2_g_dr = 1./_integral_r2_g_dr;
-}
-
-RadialGaussian::RadialGaussian() :
-	_r0(0.),
-	_sigma(0.),
-	_alpha(1./0.),
-	_integral_r2_g2_dr(0.),
-	_norm_r2_g2_dr(0.),
-	_integral_r2_g_dr(0.),
-	_norm_r2_g_dr(0.) {
-	;
-}
-
-double RadialGaussian::at(double r) {
-	double p = _alpha*(r-_r0)*(r-_r0);
-	if (p < 40) return _norm_r2_g2_dr * exp(-p);
-	else return 0.0;
-}
-
-// ======================
-// SphericalGaussian
-// ======================
-
-SphericalGaussian::SphericalGaussian(vec r0, double sigma) :
-	_r0(r0), _sigma(sigma), _alpha(1./(2*sigma*sigma)) {
-	_norm_g_dV = pow(_alpha/M_PI, 1.5);
-}
-
-// ======================
-// ModSphBessel1stKind
-// ======================
-
-std::vector<double> ModifiedSphericalBessel1stKind::eval(int degree, double r) {
-	std::vector<double> il;
-	if (r < RADZERO) {
-		il.push_back(1.);
-		il.push_back(0.);
-	}
-	else {
-		il.push_back(sinh(r)/r);
-		il.push_back(cosh(r)/r - sinh(r)/(r*r));
-	}
-	for (int l = 2; l <= degree; ++l) {
-		if (r < RADZERO) {
-			il.push_back(0.);
-		}
-		else {
-			if (il[l-1] < SPHZERO) il.push_back(0.);
-			else il.push_back( il[l-2] - (2*(l-1)+1)/r*il[l-1] );
-		}
-	}
-	return il;
-}
-
-} /* CLOSE NAMESPACE */
diff --git a/src/CMakeLists.txt b/src/soap/CMakeLists.txt
similarity index 82%
rename from src/CMakeLists.txt
rename to src/soap/CMakeLists.txt
index 7438f18..48b5e90 100644
--- a/src/CMakeLists.txt
+++ b/src/soap/CMakeLists.txt
@@ -1,6 +1,7 @@
 file(GLOB LOCAL_SOURCES *.cpp)
 
 # SUBDIRECTORIES
+add_subdirectory(base)
 add_subdirectory(linalg)
 add_subdirectory(tools)
 
@@ -11,6 +12,10 @@ message("Sources: soap/")
 foreach(item ${LOCAL_SOURCES})
     message(STATUS " o " ${item})
 endforeach()
+message("Headers: soap/")
+foreach(item ${HEADERS})
+    message(STATUS " o " ${item})
+endforeach()
 
 # COMPILE LIBRARIES
 set(LD_LIBRARIES ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${MPI_LIBRARIES})
@@ -25,5 +30,8 @@ configure_file(SOAPRC.in SOAPRC @ONLY)
 
 install(TARGETS _soapxx LIBRARY DESTINATION ${LOCAL_INSTALL_DIR})
 install(FILES __init__.py ${CMAKE_CURRENT_BINARY_DIR}/SOAPRC DESTINATION ${LOCAL_INSTALL_DIR})
-#install(DIRECTORY linalg DESTINATION ${LOCAL_INSTALL_DIR})
 
+# HEADERS
+file(GLOB HEADERS *.hpp)
+install(FILES ${HEADERS} DESTINATION ${LOCAL_INSTALL_DIR}/include/soap)
+#install(DIRECTORY linalg DESTINATION ${LOCAL_INSTALL_DIR})
diff --git a/src/SOAPRC.in b/src/soap/SOAPRC.in
similarity index 56%
rename from src/SOAPRC.in
rename to src/soap/SOAPRC.in
index 5bda694..8352beb 100644
--- a/src/SOAPRC.in
+++ b/src/soap/SOAPRC.in
@@ -4,7 +4,11 @@ test "$is_csh" = 101 && goto CSH
 
 # BASH
 export PYTHONPATH="${PYTHONPATH}:@CMAKE_INSTALL_PREFIX@"
+export SOAP_ROOT="@CMAKE_INSTALL_PREFIX@"
+export LD_LIBRARY_PATH="{LD_LIBRARY_PATH}:@CMAKE_INSTALL_PREFIX@/soap"
 
 # CSH/TCSH
 CSH:
 setenv PYTHONPATH="${PYTHONPATH}:@CMAKE_INSTALL_PREFIX@"
+setenv SOAP_ROOT="@CMAKE_INSTALL_PREFIX@"
+
diff --git a/src/__init__.py b/src/soap/__init__.py
similarity index 100%
rename from src/__init__.py
rename to src/soap/__init__.py
diff --git a/src/angularbasis.cpp b/src/soap/angularbasis.cpp
similarity index 96%
rename from src/angularbasis.cpp
rename to src/soap/angularbasis.cpp
index 89da443..f0262a9 100644
--- a/src/angularbasis.cpp
+++ b/src/soap/angularbasis.cpp
@@ -1,4 +1,4 @@
-#include "angularbasis.hpp"
+#include "soap/angularbasis.hpp"
 
 #include <math.h>
 #include <boost/math/special_functions/spherical_harmonic.hpp>
diff --git a/src/angularbasis.hpp b/src/soap/angularbasis.hpp
similarity index 93%
rename from src/angularbasis.hpp
rename to src/soap/angularbasis.hpp
index d5affe2..3d8b336 100644
--- a/src/angularbasis.hpp
+++ b/src/soap/angularbasis.hpp
@@ -5,9 +5,9 @@
 #include <math.h>
 #include <vector>
 
-#include "base/exceptions.hpp"
-#include "base/objectfactory.hpp"
-#include "options.hpp"
+#include "soap/base/exceptions.hpp"
+#include "soap/base/objectfactory.hpp"
+#include "soap/options.hpp"
 
 namespace soap {
 
diff --git a/src/atomicspectrum.cpp b/src/soap/atomicspectrum.cpp
similarity index 99%
rename from src/atomicspectrum.cpp
rename to src/soap/atomicspectrum.cpp
index 9f855f7..6d645c9 100644
--- a/src/atomicspectrum.cpp
+++ b/src/soap/atomicspectrum.cpp
@@ -3,7 +3,7 @@
 #include <boost/archive/binary_oarchive.hpp>
 #include <boost/archive/binary_iarchive.hpp>
 
-#include "atomicspectrum.hpp"
+#include "soap/atomicspectrum.hpp"
 
 namespace soap {
 
diff --git a/src/atomicspectrum.hpp b/src/soap/atomicspectrum.hpp
similarity index 92%
rename from src/atomicspectrum.hpp
rename to src/soap/atomicspectrum.hpp
index 9b50c1b..93b5ca0 100644
--- a/src/atomicspectrum.hpp
+++ b/src/soap/atomicspectrum.hpp
@@ -7,13 +7,13 @@
 #include <boost/serialization/complex.hpp>
 #include <boost/serialization/base_object.hpp>
 
-#include "base/logger.hpp"
-#include "types.hpp"
-#include "globals.hpp"
-#include "options.hpp"
-#include "structure.hpp"
-#include "basis.hpp"
-#include "power.hpp"
+#include "soap/base/logger.hpp"
+#include "soap/types.hpp"
+#include "soap/globals.hpp"
+#include "soap/options.hpp"
+#include "soap/structure.hpp"
+#include "soap/basis.hpp"
+#include "soap/power.hpp"
 
 
 namespace soap {
diff --git a/src/base/exceptions.hpp b/src/soap/base/exceptions.hpp
similarity index 100%
rename from src/base/exceptions.hpp
rename to src/soap/base/exceptions.hpp
diff --git a/src/base/logger.hpp b/src/soap/base/logger.hpp
similarity index 100%
rename from src/base/logger.hpp
rename to src/soap/base/logger.hpp
diff --git a/src/base/objectfactory.hpp b/src/soap/base/objectfactory.hpp
similarity index 100%
rename from src/base/objectfactory.hpp
rename to src/soap/base/objectfactory.hpp
diff --git a/src/basis.cpp b/src/soap/basis.cpp
similarity index 99%
rename from src/basis.cpp
rename to src/soap/basis.cpp
index 8c7a36d..d557043 100644
--- a/src/basis.cpp
+++ b/src/soap/basis.cpp
@@ -1,7 +1,7 @@
 #include <math.h>
 
-#include "linalg/numpy.hpp"
-#include "basis.hpp"
+#include "soap/linalg/numpy.hpp"
+#include "soap/basis.hpp"
 
 
 namespace soap {
diff --git a/src/basis.hpp b/src/soap/basis.hpp
similarity index 91%
rename from src/basis.hpp
rename to src/soap/basis.hpp
index 7eeafd7..65169a7 100644
--- a/src/basis.hpp
+++ b/src/soap/basis.hpp
@@ -6,13 +6,13 @@
 #include <vector>
 #include <fstream>
 
-#include "base/exceptions.hpp"
-#include "options.hpp"
-#include "globals.hpp"
-#include "structure.hpp"
-#include "angularbasis.hpp"
-#include "radialbasis.hpp"
-#include "cutoff.hpp"
+#include "soap/base/exceptions.hpp"
+#include "soap/options.hpp"
+#include "soap/globals.hpp"
+#include "soap/structure.hpp"
+#include "soap/angularbasis.hpp"
+#include "soap/radialbasis.hpp"
+#include "soap/cutoff.hpp"
 
 namespace soap {
 
diff --git a/src/bindings.cpp b/src/soap/bindings.cpp
similarity index 100%
rename from src/bindings.cpp
rename to src/soap/bindings.cpp
diff --git a/src/bindings.hpp b/src/soap/bindings.hpp
similarity index 58%
rename from src/bindings.hpp
rename to src/soap/bindings.hpp
index a920ca6..db40f4c 100644
--- a/src/bindings.hpp
+++ b/src/soap/bindings.hpp
@@ -1,9 +1,9 @@
 #ifndef _SOAP_LINALG_BINDINGS_H
 #define	_SOAP_LINALG_BINDINGS_H
 #include <boost/python.hpp>
-#include "structure.hpp"
-#include "options.hpp"
-#include "spectrum.hpp"
+#include "soap/structure.hpp"
+#include "soap/options.hpp"
+#include "soap/spectrum.hpp"
 
 namespace soap {
 
diff --git a/src/boundary.cpp b/src/soap/boundary.cpp
similarity index 92%
rename from src/boundary.cpp
rename to src/soap/boundary.cpp
index e817b56..4354396 100644
--- a/src/boundary.cpp
+++ b/src/soap/boundary.cpp
@@ -1,4 +1,4 @@
-#include "boundary.hpp"
+#include "soap/boundary.hpp"
 #include <boost/archive/binary_iarchive.hpp>
 #include <boost/archive/binary_oarchive.hpp>
 #include <boost/archive/text_iarchive.hpp>
diff --git a/src/boundary.hpp b/src/soap/boundary.hpp
similarity index 99%
rename from src/boundary.hpp
rename to src/soap/boundary.hpp
index 04687cd..acf3f32 100644
--- a/src/boundary.hpp
+++ b/src/soap/boundary.hpp
@@ -4,7 +4,7 @@
 #include <boost/serialization/base_object.hpp>
 #include <boost/serialization/export.hpp>
 
-#include "types.hpp"
+#include "soap/types.hpp"
 
 namespace soap {
 
diff --git a/src/cutoff.cpp b/src/soap/cutoff.cpp
similarity index 95%
rename from src/cutoff.cpp
rename to src/soap/cutoff.cpp
index 8d46398..a954824 100644
--- a/src/cutoff.cpp
+++ b/src/soap/cutoff.cpp
@@ -1,4 +1,4 @@
-#include "cutoff.hpp"
+#include "soap/cutoff.hpp"
 
 namespace soap {
 
diff --git a/src/cutoff.hpp b/src/soap/cutoff.hpp
similarity index 93%
rename from src/cutoff.hpp
rename to src/soap/cutoff.hpp
index 439b8eb..b06e90f 100644
--- a/src/cutoff.hpp
+++ b/src/soap/cutoff.hpp
@@ -5,10 +5,10 @@
 #include <math.h>
 #include <vector>
 
-#include "base/exceptions.hpp"
-#include "base/objectfactory.hpp"
-#include "globals.hpp"
-#include "options.hpp"
+#include "soap/base/exceptions.hpp"
+#include "soap/base/objectfactory.hpp"
+#include "soap/globals.hpp"
+#include "soap/options.hpp"
 
 namespace soap {
 
diff --git a/src/soap/functions.cpp b/src/soap/functions.cpp
new file mode 100644
index 0000000..71a25df
--- /dev/null
+++ b/src/soap/functions.cpp
@@ -0,0 +1,247 @@
+#include "soap/functions.hpp"
+#include "soap/base/exceptions.hpp"
+
+#include <math.h>
+#include <boost/math/special_functions/erf.hpp>
+#include <boost/math/special_functions/spherical_harmonic.hpp>
+
+namespace soap {
+
+// ======================
+// RadialGaussian
+// ======================
+
+RadialGaussian::RadialGaussian(double r0, double sigma)
+: _r0(r0),
+  _sigma(sigma),
+  _alpha(1./(2.*sigma*sigma)) {
+
+	// COMPUTE NORMALIZATION S g^2 r^2 dr
+	// This normalization is to be used for radial basis functions
+    double w = 2*_alpha;
+    double W0 = 2*_alpha*_r0;
+    _integral_r2_g2_dr =
+        1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
+            2*sqrt(w)*W0 +
+			sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
+			    1 - boost::math::erf<double>(-W0/sqrt(w))
+			)
+        );
+    _norm_r2_g2_dr = 1./sqrt(_integral_r2_g2_dr);
+
+    // COMPUTE NORMALIZATION S g r^2 dr
+    // This normalization is to be used for "standard" radial Gaussians
+    w = _alpha;
+	W0 = _alpha*_r0;
+	_integral_r2_g_dr =
+		1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
+			2*sqrt(w)*W0 +
+			sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
+				1 - boost::math::erf<double>(-W0/sqrt(w))
+			)
+		);
+	_norm_r2_g_dr = 1./_integral_r2_g_dr;
+}
+
+RadialGaussian::RadialGaussian() :
+	_r0(0.),
+	_sigma(0.),
+	_alpha(1./0.),
+	_integral_r2_g2_dr(0.),
+	_norm_r2_g2_dr(0.),
+	_integral_r2_g_dr(0.),
+	_norm_r2_g_dr(0.) {
+	;
+}
+
+double RadialGaussian::at(double r) {
+	double p = _alpha*(r-_r0)*(r-_r0);
+	if (p < 40) return _norm_r2_g2_dr * exp(-p);
+	else return 0.0;
+}
+
+// ======================
+// SphericalGaussian
+// ======================
+
+SphericalGaussian::SphericalGaussian(vec r0, double sigma) :
+	_r0(r0), _sigma(sigma), _alpha(1./(2*sigma*sigma)) {
+	_norm_g_dV = pow(_alpha/M_PI, 1.5);
+}
+
+// ======================
+// ModSphBessel1stKind
+// ======================
+
+ModifiedSphericalBessel1stKind::ModifiedSphericalBessel1stKind(int degree) :
+    _degree(degree) {
+    _in.reserve(degree);
+    _din.reserve(degree);
+}
+
+void ModifiedSphericalBessel1stKind::evaluate(double r, bool differentiate) {
+    _in.clear();
+    _din.clear();
+
+    _in = ModifiedSphericalBessel1stKind::eval(_degree, r);
+
+    if (differentiate) {
+        if (r < RADZERO) {
+            _din.push_back(0.0);
+            _din.push_back(1./3.);
+            for (int n = 1; n <= _degree; ++n) {
+                _din.push_back(0.);
+            }
+        }
+        else {
+            _din.push_back( _in[1] );
+            for (int n = 1; n <= _degree; ++n) {
+                _din.push_back( _in[n-1] - (n+1.)/r*_in[n] );
+            }
+        }
+    }
+
+    return;
+}
+
+std::vector<double> ModifiedSphericalBessel1stKind::eval(int degree, double r) {
+	std::vector<double> il;
+	if (r < RADZERO) {
+		il.push_back(1.);
+		il.push_back(0.);
+	}
+	else {
+		il.push_back(sinh(r)/r);
+		il.push_back(cosh(r)/r - sinh(r)/(r*r));
+	}
+	for (int l = 2; l <= degree; ++l) {
+		if (r < RADZERO) {
+			il.push_back(0.);
+		}
+		else {
+			if (il[l-1] < SPHZERO) il.push_back(0.);
+			il.push_back( il[l-2] - (2*(l-1)+1)/r*il[l-1] );
+		}
+	}
+	return il;
+}
+
+// ==============================
+// GradSphericalYlm \Nabla Y_{lm}
+// ==============================
+
+std::vector<std::complex<double> > GradSphericalYlm::eval(int l, int m, vec &r) {
+
+    std::vector<std::complex<double> > dylm;
+    dylm.push_back(std::complex<double>(0.,0.));
+    dylm.push_back(std::complex<double>(0.,0.));
+    dylm.push_back(std::complex<double>(0.,0.));
+
+    // TODO WHAT IF RADIUS IS ZERO?
+    double R = soap::linalg::abs(r);
+    if (R < RADZERO) {
+        throw soap::base::NotImplemented("<GradSphericalYlm::eval> R < RADZERO");
+    }
+    vec d = r / R;
+    double x = r.getX();
+    double y = r.getY();
+    double z = r.getZ();
+
+    // COMPUTE YLM
+    double theta = acos(d.getZ());
+    double phi = atan2(d.getY(), d.getX());
+    // Shift [-pi, -0] to [pi, 2*pi]
+    if (phi < 0.) phi += 2*M_PI;
+    cmplx ylm = boost::math::spherical_harmonic<double,double>(l, m, theta, phi);
+
+    // COMPUTE DYLM
+    int p, q, s;
+    for (p = 0; p <= l; ++p) {
+         q = p - m;
+         s = l - p - q;
+
+         //std::cout << p << "_" << q << "_" << s << std::endl;
+         //std::cout << "DYLM" << dylm[0] << dylm[1] << dylm[2] << std::endl;
+
+         double fact_p;
+         double fact_q;
+         double fact_s;
+         if ((p >= 0)) {
+             fact_p = factorial(p);
+         }
+         if (q >= 0) {
+              fact_q = factorial(q);
+         }
+         if ((s >= 0)) {
+              fact_s = factorial(s);
+         }
+
+         if ((p >= 1) && (q >= 0) && (s >= 0)) {
+             double fact_p_1 = factorial(p-1);
+            dylm[0] = dylm[0]
+                - (( pow_nnan(cmplx(-0.5 * x, -0.5 * y),(p - 1)) )
+                * ( pow_nnan(cmplx(0.5 * x, -0.5 * y),q) )
+                * ( pow(z,s) )
+                * 0.5
+                / (  fact_p_1*fact_q*fact_s  ));
+            dylm[1] = dylm[1]
+                - (( pow_nnan(cmplx(-0.5 * x, -0.5 * y),(p - 1)) )
+                * ( pow_nnan(cmplx(0.5 * x, -0.5 * y),q) )
+                * ( pow(z,s) )
+                * 0.5 * cmplx(0.0, 1.0)
+                / (  fact_p_1*fact_q*fact_s  ));
+         }
+
+         if ((p >= 0) && (q >= 1) && (s >= 0)) {
+             double fact_q_1 = factorial(q-1);
+             dylm[0] = dylm[0]
+                + (( pow_nnan(cmplx(-0.5 * x, -0.5 * y),p) )
+                * ( pow_nnan(cmplx(0.5 * x, -0.5 * y),(q - 1)) )
+                * ( pow(z,s) )
+                * 0.5
+                / ( fact_p*fact_q_1*fact_s ));
+             dylm[1] = dylm[1]
+                - (( pow_nnan(cmplx(-0.5 * x, -0.5 * y),p) )
+                * ( pow_nnan(cmplx(0.5 * x, -0.5 * y),(q - 1)) )
+                * ( pow(z,s) )
+                * 0.5 * cmplx(0.0, 1.0)
+                / ( fact_p*fact_q_1*fact_s ));
+         }
+
+         if ((p >= 0) && (q >= 0) && (s >= 1)) {
+             double fact_s_1 = factorial(s-1);
+             dylm[2] = dylm[2]
+                + (( pow_nnan(cmplx(-0.5 * x, -0.5 * y),p) )
+                * ( pow_nnan(cmplx(0.5 * x, -0.5 * y),q) )
+                * ( pow(z,(s - 1)) )
+                / ( fact_p*fact_q*fact_s_1 ));
+         }
+         //std::cout << "DYLM-OUT" << dylm[0] << dylm[1] << dylm[2] << std::endl;
+    }
+
+    double prefac = sqrt(factorial(l + m) * factorial(l - m) * ((2.0 * l) + 1) / (4.0 * M_PI)) * ( pow(R*R,(-0.5 * l)) );
+    dylm[0] = dylm[0]*prefac;
+    dylm[1] = dylm[1]*prefac;
+    dylm[2] = dylm[2]*prefac;
+
+    dylm[0] = dylm[0] - (l*x*ylm/(R*R));
+    dylm[1] = dylm[1] - (l*y*ylm/(R*R));
+    dylm[2] = dylm[2] - (l*z*ylm/(R*R));
+
+    return dylm;
+}
+
+std::complex<double> pow_nnan(std::complex<double> z, double a) {
+    if (a == 0) {
+        return std::complex<double>(1.,0.);
+    }
+    else {
+        return pow(z,a);
+    }
+}
+
+int factorial(int n) {
+    return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
+}
+
+} /* CLOSE NAMESPACE */
diff --git a/src/functions.hpp b/src/soap/functions.hpp
similarity index 66%
rename from src/functions.hpp
rename to src/soap/functions.hpp
index c368075..4571f96 100644
--- a/src/functions.hpp
+++ b/src/soap/functions.hpp
@@ -4,7 +4,7 @@
 #include <math.h>
 #include <vector>
 
-#include "types.hpp"
+#include "soap/types.hpp"
 
 namespace soap {
 
@@ -46,11 +46,28 @@ struct SphericalGaussian
 
 struct ModifiedSphericalBessel1stKind
 {
+    ModifiedSphericalBessel1stKind(int degree);
+    void evaluate(double r, bool differentiate);
+
 	static std::vector<double> eval(int degree, double r);
 	static constexpr double RADZERO = 1e-10;
 	static constexpr double SPHZERO = 1e-4;
+
+    int _degree;
+    std::vector<double> _in;
+    std::vector<double> _din;
+};
+
+struct GradSphericalYlm
+{
+    typedef std::complex<double> cmplx;
+    static std::vector<cmplx > eval(int l, int m, vec &r);
+
+    static constexpr double RADZERO = 1e-10;
 };
 
+std::complex<double> pow_nnan(std::complex<double> z, double a);
+int factorial(int n);
 
 }
 
diff --git a/src/globals.cpp b/src/soap/globals.cpp
similarity index 83%
rename from src/globals.cpp
rename to src/soap/globals.cpp
index 4da87c1..38a9332 100644
--- a/src/globals.cpp
+++ b/src/soap/globals.cpp
@@ -1,4 +1,4 @@
-#include "globals.hpp"
+#include "soap/globals.hpp"
 
 namespace soap {
 
diff --git a/src/globals.hpp b/src/soap/globals.hpp
similarity index 86%
rename from src/globals.hpp
rename to src/soap/globals.hpp
index 1fc0970..8d67aa1 100644
--- a/src/globals.hpp
+++ b/src/soap/globals.hpp
@@ -1,7 +1,7 @@
 #ifndef _SOAP_GLOBALS_H
 #define	_SOAP_GLOBALS_H
 
-#include "base/logger.hpp"
+#include "soap/base/logger.hpp"
 
 namespace soap {
 
diff --git a/src/linalg/CMakeLists.txt b/src/soap/linalg/CMakeLists.txt
similarity index 91%
rename from src/linalg/CMakeLists.txt
rename to src/soap/linalg/CMakeLists.txt
index 6e24b07..3ae84aa 100644
--- a/src/linalg/CMakeLists.txt
+++ b/src/soap/linalg/CMakeLists.txt
@@ -35,3 +35,6 @@ install(TARGETS _linalg LIBRARY DESTINATION ${LOCAL_INSTALL_DIR}/linalg)
 install(FILES __init__.py DESTINATION ${LOCAL_INSTALL_DIR}/linalg)
 #install(DIRECTORY linalg DESTINATION ${LOCAL_INSTALL_DIR})
 
+# HEADERS
+file(GLOB HEADERS *.hpp)
+install(FILES ${HEADERS} DESTINATION ${LOCAL_INSTALL_DIR}/include/soap/linalg)
diff --git a/src/linalg/__init__.py b/src/soap/linalg/__init__.py
similarity index 100%
rename from src/linalg/__init__.py
rename to src/soap/linalg/__init__.py
diff --git a/src/linalg/bindings.cpp b/src/soap/linalg/bindings.cpp
similarity index 88%
rename from src/linalg/bindings.cpp
rename to src/soap/linalg/bindings.cpp
index ed353e4..6a9af97 100644
--- a/src/linalg/bindings.cpp
+++ b/src/soap/linalg/bindings.cpp
@@ -1,4 +1,4 @@
-#include "bindings.hpp"
+#include "soap/linalg/bindings.hpp"
 
 namespace soap { namespace linalg {
 
diff --git a/src/linalg/bindings.hpp b/src/soap/linalg/bindings.hpp
similarity index 63%
rename from src/linalg/bindings.hpp
rename to src/soap/linalg/bindings.hpp
index 66e70e6..3510a3b 100644
--- a/src/linalg/bindings.hpp
+++ b/src/soap/linalg/bindings.hpp
@@ -1,7 +1,7 @@
 #ifndef _SOAP_LINALG_BINDINGS_H
 #define	_SOAP_LINALG_BINDINGS_H
-#include "vec.hpp"
-#include "matrix.hpp"
+#include "soap/linalg/vec.hpp"
+#include "soap/linalg/matrix.hpp"
 
 namespace soap { namespace linalg {
 
diff --git a/src/linalg/gsl/operations.cpp b/src/soap/linalg/gsl/operations.cpp
similarity index 98%
rename from src/linalg/gsl/operations.cpp
rename to src/soap/linalg/gsl/operations.cpp
index 8f23365..42221c6 100644
--- a/src/linalg/gsl/operations.cpp
+++ b/src/soap/linalg/gsl/operations.cpp
@@ -1,4 +1,4 @@
-#include "../operations.hpp"
+#include "soap/linalg/operations.hpp"
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_errno.h>
diff --git a/src/linalg/matrix.cpp b/src/soap/linalg/matrix.cpp
similarity index 99%
rename from src/linalg/matrix.cpp
rename to src/soap/linalg/matrix.cpp
index d00d87b..711ae34 100644
--- a/src/linalg/matrix.cpp
+++ b/src/soap/linalg/matrix.cpp
@@ -1,6 +1,6 @@
 #include <math.h>
 #include <stdlib.h>
-#include "matrix.hpp"
+#include "soap/linalg/matrix.hpp"
 
 #define PITIMES2 2*3.141592654
 
diff --git a/src/linalg/matrix.hpp b/src/soap/linalg/matrix.hpp
similarity index 99%
rename from src/linalg/matrix.hpp
rename to src/soap/linalg/matrix.hpp
index 55b979d..dd8b5f2 100644
--- a/src/linalg/matrix.hpp
+++ b/src/soap/linalg/matrix.hpp
@@ -2,7 +2,7 @@
 #define	_SOAP_MATRIX_HPP
 
 #include <ostream>
-#include "vec.hpp"
+#include "soap/linalg/vec.hpp"
 
 namespace soap { namespace linalg {
 
diff --git a/src/linalg/mkl/operations.cpp b/src/soap/linalg/mkl/operations.cpp
similarity index 98%
rename from src/linalg/mkl/operations.cpp
rename to src/soap/linalg/mkl/operations.cpp
index faa7d48..283ad81 100644
--- a/src/linalg/mkl/operations.cpp
+++ b/src/soap/linalg/mkl/operations.cpp
@@ -1,4 +1,4 @@
-#include "../operations.hpp"
+#include "soap/linalg/operations.hpp"
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <mkl.h>
 #include <mkl_lapacke.h>
diff --git a/src/linalg/numpy.hpp b/src/soap/linalg/numpy.hpp
similarity index 100%
rename from src/linalg/numpy.hpp
rename to src/soap/linalg/numpy.hpp
diff --git a/src/linalg/operations.hpp b/src/soap/linalg/operations.hpp
similarity index 100%
rename from src/linalg/operations.hpp
rename to src/soap/linalg/operations.hpp
diff --git a/src/linalg/types.hpp b/src/soap/linalg/types.hpp
similarity index 77%
rename from src/linalg/types.hpp
rename to src/soap/linalg/types.hpp
index ceeef01..27ac48e 100644
--- a/src/linalg/types.hpp
+++ b/src/soap/linalg/types.hpp
@@ -1,7 +1,7 @@
 #ifndef _SOAP_LINALG_TYPES_H
 #define	_SOAP_LINALG_TYPES_H
 
-#include "vec.hpp"
+#include "soap/linalg/vec.hpp"
 
 namespace soap { namespace linalg {
 
diff --git a/src/linalg/vec.hpp b/src/soap/linalg/vec.hpp
similarity index 100%
rename from src/linalg/vec.hpp
rename to src/soap/linalg/vec.hpp
diff --git a/src/options.cpp b/src/soap/options.cpp
similarity index 98%
rename from src/options.cpp
rename to src/soap/options.cpp
index 9e1ebe1..8388c3b 100644
--- a/src/options.cpp
+++ b/src/soap/options.cpp
@@ -1,4 +1,4 @@
-#include "options.hpp"
+#include "soap/options.hpp"
 
 
 namespace soap {
diff --git a/src/options.hpp b/src/soap/options.hpp
similarity index 98%
rename from src/options.hpp
rename to src/soap/options.hpp
index 3697bc4..4972805 100644
--- a/src/options.hpp
+++ b/src/soap/options.hpp
@@ -4,7 +4,7 @@
 #include <map>
 #include <boost/format.hpp>
 
-#include "types.hpp"
+#include "soap/types.hpp"
 
 namespace soap {
 
diff --git a/src/power.cpp b/src/soap/power.cpp
similarity index 98%
rename from src/power.cpp
rename to src/soap/power.cpp
index 46409ad..a04bae0 100644
--- a/src/power.cpp
+++ b/src/soap/power.cpp
@@ -1,5 +1,5 @@
-#include "power.hpp"
-#include "linalg/numpy.hpp"
+#include "soap/power.hpp"
+#include "soap/linalg/numpy.hpp"
 
 namespace soap {
 
diff --git a/src/power.hpp b/src/soap/power.hpp
similarity index 100%
rename from src/power.hpp
rename to src/soap/power.hpp
diff --git a/src/radialbasis.cpp b/src/soap/radialbasis.cpp
similarity index 98%
rename from src/radialbasis.cpp
rename to src/soap/radialbasis.cpp
index ae1b537..954f1bc 100644
--- a/src/radialbasis.cpp
+++ b/src/soap/radialbasis.cpp
@@ -8,10 +8,10 @@
 #include <boost/archive/text_oarchive.hpp>
 #include <boost/serialization/vector.hpp>
 
-#include "base/exceptions.hpp"
-#include "linalg/operations.hpp"
-#include "globals.hpp"
-#include "radialbasis.hpp"
+#include "soap/base/exceptions.hpp"
+#include "soap/linalg/operations.hpp"
+#include "soap/globals.hpp"
+#include "soap/radialbasis.hpp"
 
 namespace soap {
 
diff --git a/src/radialbasis.hpp b/src/soap/radialbasis.hpp
similarity index 95%
rename from src/radialbasis.hpp
rename to src/soap/radialbasis.hpp
index d1ab523..504f48d 100644
--- a/src/radialbasis.hpp
+++ b/src/soap/radialbasis.hpp
@@ -6,10 +6,10 @@
 #include <boost/serialization/base_object.hpp>
 #include <boost/serialization/export.hpp>
 
-#include "base/exceptions.hpp"
-#include "base/objectfactory.hpp"
-#include "options.hpp"
-#include "functions.hpp"
+#include "soap/base/exceptions.hpp"
+#include "soap/base/objectfactory.hpp"
+#include "soap/options.hpp"
+#include "soap/functions.hpp"
 
 namespace soap {
 
diff --git a/src/spectrum.cpp b/src/soap/spectrum.cpp
similarity index 99%
rename from src/spectrum.cpp
rename to src/soap/spectrum.cpp
index e11b414..db23c9a 100644
--- a/src/spectrum.cpp
+++ b/src/soap/spectrum.cpp
@@ -3,7 +3,7 @@
 #include <boost/archive/binary_oarchive.hpp>
 #include <boost/archive/binary_iarchive.hpp>
 
-#include "spectrum.hpp"
+#include "soap/spectrum.hpp"
 
 namespace soap {
 
diff --git a/src/spectrum.hpp b/src/soap/spectrum.hpp
similarity index 95%
rename from src/spectrum.hpp
rename to src/soap/spectrum.hpp
index ab15986..84ac31a 100644
--- a/src/spectrum.hpp
+++ b/src/soap/spectrum.hpp
@@ -7,14 +7,14 @@
 #include <boost/serialization/complex.hpp>
 #include <boost/serialization/base_object.hpp>
 
-#include "base/logger.hpp"
-#include "types.hpp"
-#include "globals.hpp"
-#include "options.hpp"
-#include "structure.hpp"
-#include "basis.hpp"
-#include "power.hpp"
-#include "atomicspectrum.hpp"
+#include "soap/base/logger.hpp"
+#include "soap/types.hpp"
+#include "soap/globals.hpp"
+#include "soap/options.hpp"
+#include "soap/structure.hpp"
+#include "soap/basis.hpp"
+#include "soap/power.hpp"
+#include "soap/atomicspectrum.hpp"
 
 
 namespace soap {
diff --git a/src/structure.cpp b/src/soap/structure.cpp
similarity index 99%
rename from src/structure.cpp
rename to src/soap/structure.cpp
index b0c9e6a..f3eeb7d 100644
--- a/src/structure.cpp
+++ b/src/soap/structure.cpp
@@ -1,4 +1,4 @@
-#include "structure.hpp"
+#include "soap/structure.hpp"
 
 
 namespace soap {
diff --git a/src/structure.hpp b/src/soap/structure.hpp
similarity index 98%
rename from src/structure.hpp
rename to src/soap/structure.hpp
index a95b5be..1624b12 100644
--- a/src/structure.hpp
+++ b/src/soap/structure.hpp
@@ -7,9 +7,9 @@
 #include <boost/python/suite/indexing/vector_indexing_suite.hpp>
 #include <iostream>
 
-#include "base/exceptions.hpp"
-#include "types.hpp"
-#include "boundary.hpp"
+#include "soap/base/exceptions.hpp"
+#include "soap/types.hpp"
+#include "soap/boundary.hpp"
 
 namespace soap {
 
diff --git a/src/tools/CMakeLists.txt b/src/soap/tools/CMakeLists.txt
similarity index 100%
rename from src/tools/CMakeLists.txt
rename to src/soap/tools/CMakeLists.txt
diff --git a/src/tools/__init__.py b/src/soap/tools/__init__.py
similarity index 100%
rename from src/tools/__init__.py
rename to src/soap/tools/__init__.py
diff --git a/src/tools/extract.py b/src/soap/tools/extract.py
similarity index 100%
rename from src/tools/extract.py
rename to src/soap/tools/extract.py
diff --git a/src/tools/inverse.py b/src/soap/tools/inverse.py
similarity index 100%
rename from src/tools/inverse.py
rename to src/soap/tools/inverse.py
diff --git a/src/tools/loadwrite.py b/src/soap/tools/loadwrite.py
similarity index 100%
rename from src/tools/loadwrite.py
rename to src/soap/tools/loadwrite.py
diff --git a/src/tools/loadwrite.py~ b/src/soap/tools/loadwrite.py~
similarity index 100%
rename from src/tools/loadwrite.py~
rename to src/soap/tools/loadwrite.py~
diff --git a/src/types.hpp b/src/soap/types.hpp
similarity index 93%
rename from src/types.hpp
rename to src/soap/types.hpp
index 4ade029..079d98d 100644
--- a/src/types.hpp
+++ b/src/soap/types.hpp
@@ -6,8 +6,8 @@
 #include <stdexcept>
 #include <boost/python.hpp>
 #include <boost/lexical_cast.hpp>
-#include "linalg/types.hpp"
-#include "linalg/matrix.hpp"
+#include "soap/linalg/types.hpp"
+#include "soap/linalg/matrix.hpp"
 
 namespace soap {
 
-- 
GitLab