diff --git a/CMakeLists.txt b/CMakeLists.txt
index 44fa0a2..6b83dea 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,8 +1,19 @@
 # Cmake config for NanoShaper
-project(NanoShaper LANGUAGES CXX C)
 cmake_minimum_required(VERSION 3.6)
+project(NanoShaper LANGUAGES CXX C)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
+cmake_policy(SET CMP0042 NEW) # MacOS RPATH on by default
+cmake_policy(SET CMP0048 NEW) # project() command manages the VERSION variables
+cmake_policy(SET CMP0054 NEW) # Only interpret if() arguments as variables or keywords when unquoted
+cmake_policy(SET CMP0077 NEW) # option() honors normal variables (i.e. does nothing if a normal variable with the same name exists)
+cmake_policy(SET CMP0083 NEW) # Pass flags needed for position-independent executables
+cmake_policy(SET CMP0091 NEW) # MSVC runtime library flags are selected by an abstraction (i.e. CMAKE_MSVC_RUNTIME_LIBRARY)
+set(CMAKE_MACOSX_RPATH 1)
+
+include(FetchContent)
+include(ExternalProject)
+
 ################################################################################
 # Build type
 ################################################################################
@@ -17,29 +28,21 @@ if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
     "Debug" "Release" "RelWithDebInfo")
 endif ()
 
-################################################################################
-# Git submodules
-################################################################################
-option(ENABLE_GIT_SUBMODULES "init and update git submodules" OFF)
-
-if(ENABLE_GIT_SUBMODULES)
-  include(${CMAKE_CURRENT_SOURCE_DIR}/CMakeLists_git.txt)
-endif()
-
-
 ################################################################################
 # SPDLog
 ################################################################################
 option(ENABLE_SPDLOG "SPDLog" OFF)
 
 if(ENABLE_SPDLOG)
+  FetchContent_Declare( spdlog
+      GIT_REPOSITORY https://github.com/gabime/spdlog.git
+      GIT_TAG 27cb4c76708608465c413f6d0e6b8d99a4d84302 # v1.14.1
+  )
+  FetchContent_MakeAvailable( spdlog )
   add_definitions(-DSPDLOG_ENABLED)
   link_libraries(
     spdlog::spdlog_header_only
   )
-  if(ENABLE_GIT_SUBMODULES)
-    add_subdirectory(lib/spdlog)
-  endif()
 endif()
 
 ################################################################################
@@ -48,19 +51,45 @@ endif()
 option(ENABLE_JSON "JSON" OFF)
 
 if(ENABLE_JSON)
+  FetchContent_Declare( json
+      GIT_REPOSITORY https://github.com/nlohmann/json.git
+      GIT_TAG bc889afb4c5bf1c0d8ee29ef35eaaf4c8bef8a5d # v3.11.2
+  )
+  FetchContent_MakeAvailable( json )
   set(JSON_BuildTests OFF CACHE INTERNAL "")
   add_definitions(-DJSON_ENABLED)
   link_libraries(
     nlohmann_json::nlohmann_json
   )
-  if(ENABLE_GIT_SUBMODULES)
-    add_subdirectory(lib/json)
-  endif()
 endif()
 
-# if (${UNIX}) 
-#   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -mtune=native") 
-# endif (${UNIX}) 
+################################################################################
+# TinyPLY
+################################################################################
+option(ENABLE_PLY "PLY" OFF)
+
+if(ENABLE_PLY)
+
+  FetchContent_Declare( tinyply
+      GIT_REPOSITORY https://github.com/ddiakopoulos/tinyply.git
+      GIT_TAG 40aa4a0ae9e9c203e11893f78b8bcaf8a50e65f0 # 2.3.4
+  )
+
+  FetchContent_GetProperties(tinyply)
+  if(NOT tinyply_POPULATED)
+    FetchContent_Populate(tinyply)
+    add_library(tinyply STATIC
+      ${tinyply_SOURCE_DIR}/source/tinyply.cpp
+      ${tinyply_SOURCE_DIR}/source/tinyply.h
+    )
+  include_directories(${tinyply_SOURCE_DIR}/source)
+  endif()
+  
+  add_definitions(-DPLY_ENABLED)
+  link_libraries(
+    tinyply
+  )
+endif()
 
 ################################################################################
 # BOOST
diff --git a/src/Configuration.h b/src/Configuration.h
index 0ffb32e..c539a1c 100644
--- a/src/Configuration.h
+++ b/src/Configuration.h
@@ -92,6 +92,7 @@ struct Configuration {
   bool vaFlag = false;
   bool computeNormals = false;
   bool saveMSMS = false;
+  bool savePLY = false;
   double sternLayer = -1.;
   int Max_Atoms_Multi_Grid = 100;
   std::string surfName;
@@ -110,7 +111,7 @@ NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(
     maxMeshDim2D, maxMeshPatches2D, NumMSMSfiles, skin_s, maxSkinDim,
     maxSkinPatches, maxSkinDim2D, maxSkinPatches2D, useFastProjection,
     savePovRay, checkDuplicatedVertices, wellShaped, probeRadius, lb, vaFlag,
-    computeNormals, saveMSMS, sternLayer, Max_Atoms_Multi_Grid, surfName)
+    computeNormals, saveMSMS, savePLY, sternLayer, Max_Atoms_Multi_Grid, surfName)
 #endif
 
 using ConfigurationOP = std::shared_ptr<Configuration>;
diff --git a/src/Surface.cpp b/src/Surface.cpp
index e1e4237..bc65b8b 100644
--- a/src/Surface.cpp
+++ b/src/Surface.cpp
@@ -8,10 +8,17 @@
 #include <Surface.h>
 #include <logging.h>
 #include <array>
+#include <cstdint>
+#include <ostream>
 #include <stdexcept>
+#include <streambuf>
 #include <tuple>
 #include <vector>
 
+#ifdef PLY_ENABLED
+#include <tinyply.h>
+#endif // PLY_ENABLED
+
 namespace nanoshaper {
 
 void Surface::init() {
@@ -53,6 +60,7 @@ void Surface::init() {
   vertexAtomsMapFlag = false;
   computeNormals = false;
   saveMSMS = false;
+  savePLY = false;
   providesAnalyticalNormals = false;
   activeCubes = NULL;
   MAX_ATOMS_MULTI_GRID = 100;
@@ -68,6 +76,7 @@ void Surface::init(ConfigurationOP cf) {
   bool vaFlag = cf->vaFlag;
   bool computeNormals = cf->computeNormals;
   bool saveMSMS = cf->saveMSMS;
+  bool savePLY = cf->savePLY;
   double sternLayer = cf->sternLayer;
   MAX_ATOMS_MULTI_GRID = cf->Max_Atoms_Multi_Grid;
 
@@ -80,6 +89,7 @@ void Surface::init(ConfigurationOP cf) {
   setVertexAtomsMap(vaFlag);
   setComputeNormals(computeNormals);
   setSaveMSMS(saveMSMS);
+  setSavePLY(savePLY);
 
   // if >0 enable stern layer, else disabled by default
   if (sternLayer > 0)
@@ -4405,6 +4415,60 @@ void Surface::approximateNormals(vector<int>& appNormals, bool doOnlyList) {
 
   deleteMatrix2D<double>(nt, planes);
 }
+
+bool Surface::savePLYMesh(int format, bool revert, const char* fileName,
+                       vector<double*>& vertList, vector<int*>& triList,
+                       vector<double*>& normalsList) {
+  int numVertexes = (int)vertList.size();
+  int numTriangles = (int)triList.size();
+  char fullName[100];
+
+  snprintf(fullName, sizeof(fullName), "%s.ply", fileName);
+
+  std::filebuf fb;
+  fb.open(fullName, std::ios::out | std::ios::binary);
+  std::ostream outstream(&fb);
+  if (outstream.fail()) {
+    logging::log<logging::level::warn>("Cannot write file {}", fileName);
+    return false;
+  }
+
+  struct double3 { double x, y, z; };
+  struct int3 { int32_t x, y, z; };
+
+  std::vector<double3> vertStructVec;
+  for (double* ptr : vertList) {
+    vertStructVec.push_back(*reinterpret_cast<double3*>(ptr));
+  }
+  std::vector<int3> triStructVec;
+  for (int* ptr : triList) {
+    triStructVec.push_back(*reinterpret_cast<int3*>(ptr));
+  }
+
+  tinyply::PlyFile mesh_ply;
+  mesh_ply.add_properties_to_element("vertex", {"x", "y", "z"} ,
+                                      tinyply::Type::FLOAT64, numVertexes,
+                                      reinterpret_cast<uint8_t *>(vertStructVec.data()),
+                                      tinyply::Type::INVALID, 0);
+  mesh_ply.add_properties_to_element("face", { "vertex_indices" },
+                                      tinyply::Type::INT32, numTriangles,
+                                      reinterpret_cast<uint8_t *>(triStructVec.data()),
+                                      tinyply::Type::UINT8, 3);
+  if (normalsList.size() != 0) {
+    std::vector<double3> normalsStructVec;
+    for (double* ptr : vertList) {
+      normalsStructVec.push_back(*reinterpret_cast<double3*>(ptr));
+    }
+    mesh_ply.add_properties_to_element("vertex", { "nx", "ny", "nz" },
+                                         tinyply::Type::FLOAT64, numVertexes,
+                                         reinterpret_cast<uint8_t*>(normalsStructVec.data()),
+                                         tinyply::Type::INVALID, 0);
+  }
+
+  mesh_ply.write(outstream, true);
+  return true;
+}
+
 bool Surface::saveMesh(int format, bool revert, const char* fileName,
                        vector<double*>& vertList, vector<int*>& triList,
                        vector<double*>& normalsList) {
@@ -4412,6 +4476,10 @@ bool Surface::saveMesh(int format, bool revert, const char* fileName,
   int numTriangles = (int)triList.size();
   char fullName[100];
 
+  if (format == PLY) {
+    return savePLYMesh(format, revert, fileName, vertList, triList, normalsList);
+  }
+
   if (format == OFF || format == OFF_A || format == OFF_N ||
       format == OFF_N_A) {
     // save all in OFF format
@@ -5204,22 +5272,26 @@ void Surface::smoothSurface(const char* fn, bool revert) {
 
 int Surface::deduceFormat() {
   int format = -1;
-  if (!saveMSMS) {
-    if (vertexAtomsMapFlag && computeNormals)
-      format = OFF_N_A;
-    else {
-      if (vertexAtomsMapFlag)
-        format = OFF_A;
-      else if (computeNormals)
-        format = OFF_N;
-      else
-        format = OFF;
-    }
-  } else {
+  if (savePLY) {
+    format = PLY;
+    return format;
+  }
+  if (saveMSMS) {
     if (vertexAtomsMapFlag)
       format = MSMS;
     else
       format = MSMS_NO_A;
+    return format;
+  }
+  if (vertexAtomsMapFlag && computeNormals)
+    format = OFF_N_A;
+  else {
+    if (vertexAtomsMapFlag)
+      format = OFF_A;
+    else if (computeNormals)
+      format = OFF_N;
+    else
+      format = OFF;
   }
   return format;
 }
diff --git a/src/Surface.h b/src/Surface.h
index 0898fb1..31b4da3 100644
--- a/src/Surface.h
+++ b/src/Surface.h
@@ -135,6 +135,7 @@ enum FileFormat : int {
   OFF_N_A = 3,
   MSMS_NO_A = 4,
   MSMS = 5,
+  PLY = 6,
   N_FORMATS
 };
 // #define DEDUCE -1
@@ -254,6 +255,7 @@ class Surface {
   bool useLoadBalancing;
   bool vertexAtomsMapFlag;
   bool saveMSMS;
+  bool savePLY;
   //////////////////////////////////////////////////////////////////////////////////////
 
   // current panel under analysis
@@ -562,6 +564,10 @@ class Surface {
   virtual bool saveMesh(const char* filename, bool revert = false,
                         int format = FileFormat::DEDUCE);
 
+  virtual bool savePLYMesh(int format, bool revert, const char* fileName,
+                        vector<double*>& vertList, vector<int*>& triList,
+                        vector<double*>& normalsList);
+
   /** smooth a given mesh and overwrites the given file name.
   The input/output mesh is in .off format*/
   virtual void smoothSurface(const char* fn = "triangulatedSurf",
@@ -645,6 +651,10 @@ class Surface {
 
   bool getSaveMSMS() { return saveMSMS; }
 
+  void setSavePLY(bool m) { savePLY = m; }
+
+  bool getSavePLY() { return savePLY; }
+
   void setTriangulationFlag(bool flag) { accurateTriangulation = flag; }
   bool getTriangulationFlag() { return accurateTriangulation; }
 
diff --git a/src/main_functions.cpp b/src/main_functions.cpp
index de8185c..bfe7ab5 100644
--- a/src/main_functions.cpp
+++ b/src/main_functions.cpp
@@ -193,6 +193,7 @@ ConfigurationOP parse(ConfigFileOP cf) {
   conf->vaFlag = cf->read<bool>("Vertex_Atom_Info", false);
   conf->computeNormals = cf->read<bool>("Compute_Vertex_Normals", false);
   conf->saveMSMS = cf->read<bool>("Save_Mesh_MSMS_Format", false);
+  conf->savePLY = cf->read<bool>("Save_Mesh_PLY_Format", false);
   conf->sternLayer = cf->read<double>("Stern_layer", -1.);
   conf->Max_Atoms_Multi_Grid = cf->read<int>("Max_Atoms_Multi_Grid", 100);
   conf->surfName = cf->read<string>("Surface");