Fix some inconsistencies and a potential memory leak

Author: Roberto Di Remigio

CHANGELOG.md

@@ -17,6 +17,9 @@
 
 ### Fixed
 
+- Some inconsistencies in input reading from host and a related memory leak in the radii
+  initialization.
+
 ### Security
 
 ## [Version 1.1.5] - 2016-07-19

src/interface/Input.cpp

@@ -217,7 +217,7 @@ void Input::reader(const PCMInput & host_input)
     mode_ = std::string("IMPLICIT");
 
     std::string name = trim_and_upper(host_input.solvent);
-    if (name.empty()) {
+    if (name.empty() || name == "EXPLICIT") {
         hasSolvent_ = false;
         // Get the probe radius
         probeRadius_ = host_input.probe_radius * angstromToBohr();
@@ -259,6 +259,13 @@ void Input::reader(const PCMInput & host_input)
     isDynamic_ = false;
 
     providedBy_ = std::string("host-side");
+
+    // Fill the input wrapping structs
+    insideGreenData_ = greenData(derivativeInsideType_, integratorType_, profileType_, epsilonInside_, integratorScaling_);
+    outsideStaticGreenData_ = greenData(derivativeOutsideType_,
+        integratorType_, profileType_,  epsilonStaticOutside_, integratorScaling_);
+    outsideDynamicGreenData_ = greenData(derivativeOutsideType_,
+        integratorType_, profileType_, epsilonDynamicOutside_, integratorScaling_);
 }
 
 void Input::semanticCheck()
@@ -280,9 +287,9 @@ void Input::initMolecule()
     j += 4;
   }
   // 3. list of atoms and list of spheres
-  double factor = angstromToBohr();
   std::vector<Atom> radiiSet;
   std::vector<Atom> atoms;
+  atoms.reserve(nuclei);
   if ( radiiSet_ == "UFF" ) {
     radiiSet = initUFF();
     radiiSetName_ = "UFF";
@@ -296,15 +303,15 @@ void Input::initMolecule()
   for (int i = 0; i < charges.size(); ++i) {
     int index = int(charges(i)) - 1;
     atoms.push_back(radiiSet[index]);
-    if (scaling_) atoms[index].radiusScaling = 1.2;
+    if (scaling_) atoms[i].radiusScaling = 1.2;
   }
   // Based on the creation mode (Implicit or Atoms)
   // the spheres list might need postprocessing
   if ( mode_ == "IMPLICIT" || mode_ == "ATOMS") {
     for (int i = 0; i < charges.size(); ++i) {
-      int index = int(charges(i)) - 1;
-      double radius = radiiSet[index].radius * factor;
-      if (scaling_) radius *= 1.2;
+      // Convert to Bohr and multiply by scaling factor (alpha)
+      double radius = atoms[i].radius * angstromToBohr()
+        * atoms[i].radiusScaling;
       spheres_.push_back(Sphere(centers.col(i), radius));
     }
     if (mode_ == "ATOMS") {
@@ -348,48 +355,45 @@ cavityData Input::cavityParams()
 greenData Input::insideGreenParams()
 {
     if (insideGreenData_.empty) {
-        int profile = profilePolicy("UNIFORM");
-        insideGreenData_ = greenData(derivativeInsideType_, integratorType_, profile, epsilonInside_, integratorScaling_);
+      insideGreenData_ = greenData(derivativeInsideType_, integratorType_, profileType_, epsilonInside_, integratorScaling_);
     }
     return insideGreenData_;
 }
 
 greenData Input::outsideStaticGreenParams()
 {
-    if (outsideStaticGreenData_.empty) {
-        int profile = profilePolicy("UNIFORM");
-        outsideStaticGreenData_ = greenData(derivativeOutsideType_,
-                                            integratorType_, profile,  epsilonStaticOutside_, integratorScaling_);
-        if (not hasSolvent_) {
-           outsideStaticGreenData_.howProfile = profileType_;
-           outsideStaticGreenData_.epsilon1 = epsilonStatic1_;
-           outsideStaticGreenData_.epsilon2 = epsilonStatic2_;
-           outsideStaticGreenData_.center   = center_;
-           outsideStaticGreenData_.width    = width_;
-           outsideStaticGreenData_.origin   << origin_[0], origin_[1], origin_[2];
-           outsideStaticGreenData_.maxL = maxL_;
-        }
+  if (outsideStaticGreenData_.empty) {
+    outsideStaticGreenData_ = greenData(derivativeOutsideType_,
+        integratorType_, profileType_,  epsilonStaticOutside_, integratorScaling_);
+    if (not hasSolvent_) {
+      outsideStaticGreenData_.howProfile = profileType_;
+      outsideStaticGreenData_.epsilon1 = epsilonStatic1_;
+      outsideStaticGreenData_.epsilon2 = epsilonStatic2_;
+      outsideStaticGreenData_.center   = center_;
+      outsideStaticGreenData_.width    = width_;
+      outsideStaticGreenData_.origin   << origin_[0], origin_[1], origin_[2];
+      outsideStaticGreenData_.maxL = maxL_;
     }
-    return outsideStaticGreenData_;
+  }
+  return outsideStaticGreenData_;
 }
 
 greenData Input::outsideDynamicGreenParams()
 {
-    if (outsideDynamicGreenData_.empty) {
-        int profile = profilePolicy("UNIFORM");
-        outsideDynamicGreenData_ = greenData(derivativeOutsideType_,
-                                             integratorType_, profile, epsilonDynamicOutside_, integratorScaling_);
-        if (not hasSolvent_) {
-           outsideDynamicGreenData_.howProfile  = profileType_;
-           outsideDynamicGreenData_.epsilon1 = epsilonDynamic1_;
-           outsideDynamicGreenData_.epsilon2 = epsilonDynamic2_;
-           outsideDynamicGreenData_.center   = center_;
-           outsideDynamicGreenData_.width    = width_;
-           outsideDynamicGreenData_.origin   << origin_[0], origin_[1], origin_[2];
-           outsideDynamicGreenData_.maxL = maxL_;
-        }
+  if (outsideDynamicGreenData_.empty) {
+    outsideDynamicGreenData_ = greenData(derivativeOutsideType_,
+        integratorType_, profileType_, epsilonDynamicOutside_, integratorScaling_);
+    if (not hasSolvent_) {
+      outsideDynamicGreenData_.howProfile  = profileType_;
+      outsideDynamicGreenData_.epsilon1 = epsilonDynamic1_;
+      outsideDynamicGreenData_.epsilon2 = epsilonDynamic2_;
+      outsideDynamicGreenData_.center   = center_;
+      outsideDynamicGreenData_.width    = width_;
+      outsideDynamicGreenData_.origin   << origin_[0], origin_[1], origin_[2];
+      outsideDynamicGreenData_.maxL = maxL_;
     }
-    return outsideDynamicGreenData_;
+  }
+  return outsideDynamicGreenData_;
 }
 
 solverData Input::solverParams()

src/interface/Meddle.cpp

@@ -459,7 +459,6 @@ namespace pcm {
   {
     bool scaling = inp.scaling();
     std::string set = inp.radiiSet();
-    double factor = angstromToBohr();
     std::vector<Atom> radiiSet;
     std::vector<Atom> atoms;
     if ( set == "UFF" ) {
@@ -470,12 +469,15 @@ namespace pcm {
       radiiSet = initAllinger();
     }
     std::vector<Sphere> spheres;
+    atoms.reserve(nuclei);
+    spheres.reserve(nuclei);
     for (int i = 0; i < charges.size(); ++i) {
       int index = int(charges(i)) - 1;
       atoms.push_back(radiiSet[index]);
-      if (scaling) atoms[index].radiusScaling = 1.2;
-      double radius = radiiSet[index].radius * factor;
-      if (scaling) radius *= 1.2;
+      if (scaling) atoms[i].radiusScaling = 1.2;
+      // Convert to Bohr and multiply by scaling factor (alpha)
+      double radius = atoms[i].radius * angstromToBohr()
+        * atoms[i].radiusScaling;
       spheres.push_back(Sphere(centers.col(i), radius));
     }
     Eigen::VectorXd masses = Eigen::VectorXd::Zero(nuclei);

src/utils/Atom.hpp

@@ -2,22 +2,22 @@
 /*
  *     PCMSolver, an API for the Polarizable Continuum Model
  *     Copyright (C) 2013-2016 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.io/>
  */
@@ -57,10 +57,11 @@ struct Atom
   /*! Atomic symbol */
   std::string symbol;
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW /* See http://eigen.tuxfamily.org/dox/group__TopicStructHavingEigenMembers.html */
-  Atom() {}
+  Atom() : charge(0.0), mass(0.0), radius(0.0), radiusScaling(0.0),
+      position(Eigen::Vector3d::Zero()), element("Dummy"), symbol("Du") {}
   Atom(const std::string & elem, const std::string & sym,
       double c, double m, double r,
-      const Eigen::Vector3d & coord, double scal = 1.2)
+      const Eigen::Vector3d & coord, double scal = 1.0)
     : charge(c), mass(m), radius(r), radiusScaling(scal),
       position(coord), element(elem), symbol(sym) {}
 };