Modify way lenght conversions are handled.

This is an ugly development (since we have to rely on a global object)
We should rather rely on Boost.Units, which has the advantage of performing compile-time
dimensional analysis

Author: Roberto Di Remigio

cmake/custom/autogenerated.cmake

@@ -1,7 +1,6 @@
 # Configure the header with library-wide preprocessor definitions
 configure_file(${PROJECT_SOURCE_DIR}/include/Config.hpp.in
                ${PROJECT_BINARY_DIR}/include/tmp-config-hpp @ONLY)
-             #file(STRINGS ${PROJECT_BINARY_DIR}/include/tmp-config-hpp _contents NEWLINE_CONSUME)
 add_custom_command(
   OUTPUT
     ${PROJECT_BINARY_DIR}/include/Config.hpp
@@ -10,7 +9,7 @@ add_custom_command(
   VERBATIM
   )
 set_source_files_properties(${PROJECT_BINARY_DIR}/include/Config.hpp PROPERTIES GENERATED TRUE)
-add_custom_target(generate-config-hpp DEPENDS ${PROJECT_BINARY_DIR}/include/Config.hpp)
+add_custom_target(generate-config-hpp ALL DEPENDS ${PROJECT_BINARY_DIR}/include/Config.hpp)
 if(BUILD_CUSTOM_BOOST)
   add_dependencies(generate-config-hpp custom_boost)
 endif()

include/Config.hpp.in

@@ -46,6 +46,7 @@
 #include "Cxx11Workarounds.hpp"
 #include "ErrorHandling.hpp"
 #include "LoggerInterface.hpp"
+#include "PhysicalConstants.hpp"
 #include "TimerInterface.hpp"
 
 #endif // CONFIG_HPP

include/PhysicalConstants.hpp

@@ -0,0 +1,57 @@
+/* pcmsolver_copyright_start */
+/*
+ *     PCMSolver, an API for the Polarizable Continuum Model
+ *     Copyright (C) 2013-2015 Roberto Di Remigio, Luca Frediani and contributors
+ *
+ *     This file is part of PCMSolver.
+ *
+ *     PCMSolver is free software: you can redistribute it and/or modify
+ *     it under the terms of the GNU Lesser General Public License as published by
+ *     the Free Software Foundation, either version 3 of the License, or
+ *     (at your option) any later version.
+ *
+ *     PCMSolver is distributed in the hope that it will be useful,
+ *     but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *     GNU Lesser General Public License for more details.
+ *
+ *     You should have received a copy of the GNU Lesser General Public License
+ *     along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *     For information on the complete list of contributors to the
+ *     PCMSolver API, see: <http://pcmsolver.readthedocs.org/>
+ */
+/* pcmsolver_copyright_end */
+
+#ifndef PHYSICALCONSTANTS_HPP
+#define PHYSICALCONSTANTS_HPP
+
+#include <cmath>
+
+/*! \file PhysicalConstants.hpp
+ *  \brief This header file contains physical constants to be used throughout the module.
+ *  \author Roberto Di Remigio
+ *  \date 2013
+ */
+
+struct LengthConversion
+{
+  double operator()() { return BOHR_TO_ANGSTROM; }
+  double BOHR_TO_ANGSTROM;
+};
+
+extern LengthConversion bohrToAngstrom;
+
+void initBohrToAngstrom(LengthConversion & conversion, int year = 2010);
+
+double angstromToBohr();
+
+double bohr2ToAngstrom2();
+
+double angstrom2ToBohr2();
+
+double bohr3ToAngstrom3();
+
+double angstrom3ToBohr3();
+
+#endif // PHYSICALCONSTANTS_HPP

src/bin/plot_green_spherical.cpp

@@ -35,7 +35,6 @@
 #include "DerivativeTypes.hpp"
 #include "GeneralPurpose.hpp"
 #include "OneLayerTanh.hpp"
-#include "PhysicalConstants.hpp"
 #include "SphericalDiffuse.hpp"
 #include "UniformDielectric.hpp"
 

src/bin/run_pcm.cpp

@@ -51,7 +51,6 @@
 #include "utils/Atom.hpp"
 #include "Citation.hpp"
 #include "utils/cnpy.hpp"
-#include "utils/PhysicalConstants.hpp"
 #include "utils/Solvent.hpp"
 #include "utils/Sphere.hpp"
 

src/interface/CMakeLists.txt

@@ -2,7 +2,7 @@
 list(APPEND headers_list Meddle.hpp Input.hpp)
 
 # List of sources
-list(APPEND sources_list Meddle.cpp Input.cpp)
+list(APPEND sources_list Meddle.cpp Input.cpp PhysicalConstants.cpp)
 
 add_library(interface OBJECT ${sources_list} ${headers_list})
 set_target_properties(interface PROPERTIES POSITION_INDEPENDENT_CODE 1 INTERPROCEDURAL_OPTIMIZATION 1)

src/interface/Input.cpp

@@ -41,7 +41,6 @@
 #include "green/GreenData.hpp"
 #include "solver/SolverData.hpp"
 #include "PCMInput.h"
-#include "utils/PhysicalConstants.hpp"
 #include "utils/Solvent.hpp"
 #include "utils/Sphere.hpp"
 
@@ -65,6 +64,7 @@ void Input::reader(const std::string & filename)
     Getkw input_ = Getkw(filename, false, true);
     units_      = input_.getStr("UNITS");
     CODATAyear_ = input_.getInt("CODATA");
+    initBohrToAngstrom(bohrToAngstrom, CODATAyear_);
 
     const Section & mol = input_.getSect("MOLECULE");
     if (mol.isDefined()) geometry_ = mol.getDblVec("GEOMETRY");
@@ -149,7 +149,7 @@ void Input::reader(const std::string & filename)
         } else {
             solvent_ = solvents()[name];
         }
-        probeRadius_ = solvent_.probeRadius * angstromToBohr(CODATAyear_);
+        probeRadius_ = solvent_.probeRadius * angstromToBohr();
         // Specification of the solvent by name means isotropic PCM
         // We have to initialize the Green's functions data here, Solvent class
         // is an helper class and should not be used in the core classes.
@@ -187,27 +187,29 @@ std::string trim_and_upper(const char * src)
 void Input::reader(const PCMInput & host_input)
 {
     CODATAyear_ = 2010;
+    std::cout << "Calling initBohrToAngstrom" << std::endl;
+    initBohrToAngstrom(bohrToAngstrom, CODATAyear_);
 
     type_ = trim_and_upper(host_input.cavity_type);
-    area_ = host_input.area * angstrom2ToBohr2(CODATAyear_);
+    area_ = host_input.area * angstrom2ToBohr2();
     patchLevel_ = host_input.patch_level;
-    coarsity_ = host_input.coarsity * angstromToBohr(CODATAyear_);
-    minDistance_ = host_input.min_distance * angstromToBohr(CODATAyear_);
+    coarsity_ = host_input.coarsity * angstromToBohr();
+    minDistance_ = host_input.min_distance * angstromToBohr();
     derOrder_ = host_input.der_order;
     if (type_ == "RESTART") {
         cavFilename_ = trim(host_input.restart_name); // No case conversion here!
     }
 
     scaling_ = host_input.scaling;
     radiiSet_ = trim_and_upper(host_input.radii_set);
-    minimalRadius_ = host_input.min_radius * angstromToBohr(CODATAyear_);
+    minimalRadius_ = host_input.min_radius * angstromToBohr();
     mode_ = std::string("IMPLICIT");
 
     std::string name = trim_and_upper(host_input.solvent);
     if (name.empty()) {
         hasSolvent_ = false;
         // Get the probe radius
-        probeRadius_ = host_input.probe_radius * angstromToBohr(CODATAyear_);
+        probeRadius_ = host_input.probe_radius * angstromToBohr();
         // Get the contents of the Green<inside> section...
         // ...and initialize the data members
         greenInsideType_ = trim_and_upper(host_input.inside_type);
@@ -223,7 +225,7 @@ void Input::reader(const PCMInput & host_input)
         // Just initialize the solvent object in this class
         hasSolvent_ = true;
         solvent_ = solvents()[name];
-        probeRadius_ = solvent_.probeRadius* angstromToBohr(CODATAyear_);
+        probeRadius_ = solvent_.probeRadius* angstromToBohr();
         // Specification of the solvent by name means isotropic PCM
         // We have to initialize the Green's functions data here, Solvent class
         // is an helper class and should not be used in the core classes.
@@ -266,7 +268,7 @@ void Input::initMolecule()
         j += 4;
     }
     // 3. list of atoms and list of spheres
-    double factor = angstromToBohr(CODATAyear_);
+    double factor = angstromToBohr();
     std::vector<Atom> radiiSet, atoms;
     if ( radiiSet_ == "UFF" ) {
         radiiSet = initUFF();

src/interface/Meddle.cpp

@@ -45,7 +45,6 @@
 #include "utils/Atom.hpp"
 #include "Citation.hpp"
 #include "utils/cnpy.hpp"
-#include "utils/PhysicalConstants.hpp"
 #include "utils/Solvent.hpp"
 #include "utils/Sphere.hpp"
 
@@ -426,7 +425,7 @@ namespace pcm {
   {
     bool scaling = inp.scaling();
     std::string set = inp.radiiSet();
-    double factor = angstromToBohr(inp.CODATAyear());
+    double factor = angstromToBohr();
     std::vector<Atom> radiiSet, atoms;
     if ( set == "UFF" ) {
       radiiSet = initUFF();

src/interface/PhysicalConstants.cpp

@@ -0,0 +1,88 @@
+/* pcmsolver_copyright_start */
+/*
+ *     PCMSolver, an API for the Polarizable Continuum Model
+ *     Copyright (C) 2013-2015 Roberto Di Remigio, Luca Frediani and contributors
+ *
+ *     This file is part of PCMSolver.
+ *
+ *     PCMSolver is free software: you can redistribute it and/or modify
+ *     it under the terms of the GNU Lesser General Public License as published by
+ *     the Free Software Foundation, either version 3 of the License, or
+ *     (at your option) any later version.
+ *
+ *     PCMSolver is distributed in the hope that it will be useful,
+ *     but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *     GNU Lesser General Public License for more details.
+ *
+ *     You should have received a copy of the GNU Lesser General Public License
+ *     along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *     For information on the complete list of contributors to the
+ *     PCMSolver API, see: <http://pcmsolver.readthedocs.org/>
+ */
+/* pcmsolver_copyright_end */
+
+#include "PhysicalConstants.hpp"
+
+LengthConversion bohrToAngstrom;
+
+void initBohrToAngstrom(LengthConversion & conversion, int year)
+{
+  switch(year) {
+    case 2010:
+      /*
+       * P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys. 84(4), 1527-1605 (2012)
+       */
+      conversion.BOHR_TO_ANGSTROM = 0.52917721092;
+      break;
+    case 2006:
+      /*
+       * P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys. 80(2), 633-730 (2008)
+       */
+      conversion.BOHR_TO_ANGSTROM = 0.52917720859;
+      break;
+    case 2002:
+      /* Choose this to match DALTON values for physical constants
+       * P. J. Mohr and B. N. Taylor, Rev. Mod. Phys. 77(1), 1-107 (2005)
+       */
+      conversion.BOHR_TO_ANGSTROM = 0.5291772108;
+      break;
+    case 1998:
+      /*
+       * P. J. Mohr and B. N. Taylor, Rev. Mod. Phys. 72(2), 351-495 (2000)
+       */
+      conversion.BOHR_TO_ANGSTROM = 0.5291772083;
+      break;
+    default:
+      /*
+       * Use latest values by default */
+      conversion.BOHR_TO_ANGSTROM = 0.52917721092;
+      break;
+  }
+}
+
+double angstromToBohr()
+{
+  return (1.0 / bohrToAngstrom());
+}
+
+double bohr2ToAngstrom2()
+{
+  return std::pow(bohrToAngstrom(), 2);
+}
+
+double angstrom2ToBohr2()
+{
+  return (1.0 / bohr2ToAngstrom2());
+}
+
+double bohr3ToAngstrom3()
+{
+  return std::pow(bohrToAngstrom(), 3);
+}
+
+double angstrom3ToBohr3()
+{
+  return (1.0 / bohr3ToAngstrom3());
+}

src/utils/CMakeLists.txt

@@ -1,7 +1,7 @@
 # List of headers
 list(APPEND headers_list Atom.hpp Factory.hpp ForId.hpp FortranCUtils.hpp Logger.hpp LoggerImpl.hpp MathUtils.hpp Molecule.hpp
-                         PhysicalConstants.hpp QuadratureRules.hpp Solvent.hpp Sphere.hpp SplineFunction.hpp
-                         Stencils.hpp Symmetry.hpp TestingMolecules.hpp Timer.hpp cnpy.hpp)
+                         QuadratureRules.hpp Solvent.hpp Sphere.hpp SplineFunction.hpp
+                         Stencils.hpp Symmetry.hpp Timer.hpp cnpy.hpp)
 set_source_files_properties(cnpy.hpp PROPERTIES COMPILE_FLAGS -w)
 
 # List of sources