@@ -624,6 +624,54 @@ void initAtoms(Eigen::VectorXd & charges_, Eigen::Matrix3Xd & sphereCenter_)
     collect_atoms(chg, centers);
 }
 
+void initMolecule(Molecule & molecule_)
+{
+    // Gather information necessary to build molecule_
+    // 1. number of atomic centers
+    int nuclei;
+    collect_nctot(&nuclei);
+    // 2. position and charges of atomic centers
+    Eigen::Matrix3Xd centers = Eigen::Matrix3Xd::Zero(Eigen::NoChange, nuclei);
+    Eigen::VectorXd charges  = Eigen::VectorXd::Zero(nuclei);
+    double * chg = charges.data();
+    double * pos = centers.data();
+    collect_atoms(chg, pos);
+    // 3. list of atoms and list of spheres
+    bool scaling = parsedInput->scaling();
+    std::string set = parsedInput->radiiSet();
+    double factor = angstromToBohr(parsedInput->CODATAyear());
+    std::vector<Atom> radiiSet, atoms;
+    if ( set == "UFF" ) {
+        radiiSet = Atom::initUFF();
+    } else {
+        radiiSet = Atom::initBondi();
+    }
+    std::vector<Sphere> spheres;
+    for (int i = 0; i < charges.size(); ++i) {
+        int index = int(charges(i)) - 1;
+	atoms.push_back(radiiSet[index]);
+        double radius = radiiSet[index].atomRadius() * factor;
+        if (scaling) {
+            radius *= radiiSet[index].atomRadiusScaling();
+        }
+        spheres.push_back(Sphere(centers.col(i), radius));
+    }
+    // 4. masses
+    Eigen::VectorXd masses = Eigen::VectorXd::Zero(nuclei);
+    for (int i = 0; i < masses.size(); ++i) {
+	 masses(i) = atoms[i].atomMass();
+    }
+    // Based on the creation mode (Implicit or Atoms)
+    // the spheres list might need postprocessing
+    std::string _mode = parsedInput->mode();
+    if ( _mode == "ATOMS" ) {
+       initSpheresAtoms(centers, spheres);	 
+    }
+
+    // OK, now get molecule_
+    molecule_ = Molecule(nuclei, charges, masses, centers, atoms, spheres);
+}
+
 void initSpheresAtoms(const Eigen::Matrix3Xd & sphereCenter_,
                       std::vector<Sphere> & spheres_)
 {
 
@@ -40,6 +40,7 @@
 #include <Eigen/Dense>
 
 #include "InputManager.hpp"
+#include "Molecule.hpp"
 #include "Sphere.hpp"
 
 /*
@@ -311,6 +312,13 @@ void initWaveletCavity();
  */
 void initAtoms(Eigen::VectorXd & charges_, Eigen::Matrix3Xd & sphereCenter_);
 
+/*! \fn void initMolecule(Molecule & molecule_, Eigen::Matrix3Xd & sphereCenter_)
+ *  \param[out] molecule_  a Molecule object 
+ *
+ *  Initializes Molecule object with all the information on charges, masses and positions of atomic centers
+ */
+void initMolecule(Molecule & molecule_);
+
 /*! \fn void initSpheresImplicit(const Eigen::VectorXd & charges_, const Eigen::Matrix3Xd & sphereCenter_, std::vector<Sphere> & spheres_)
  *  \param[in] charges_      contains charges of atomic centers
  *  \param[in] sphereCenter_ contains coordinates of atomic centers
 
@@ -142,14 +142,15 @@ void Input::reader(const char * pythonParsed)
     if (mode_ == "EXPLICIT") {
         std::vector<double> spheresInput = cavity.getDblVec("SPHERES");
         int j = 0;
-        int upperBound = int(spheresInput.size() / 4);
-        for (int i = 0; i < upperBound; ++i) {
+        int nAtoms = int(spheresInput.size() / 4);
+        for (int i = 0; i < nAtoms; ++i) {
             Eigen::Vector3d center;
             center << spheresInput[j], spheresInput[j+1], spheresInput[j+2];
             Sphere sph(center, spheresInput[j+3]);
             spheres_.push_back(sph);
             j += 4;
         }
+	molecule_ = Molecule(nAtoms, spheres_);
     } else if (mode_ == "ATOMS") {
         atoms_ = cavity.getIntVec("ATOMS");
         radii_ = cavity.getDblVec("RADII");
 
@@ -35,6 +35,7 @@
 
 #include "Compression.hpp"
 #include "InputManager.hpp"
+#include "Molecule.hpp"
 #include "Solvent.hpp"
 #include "Sphere.hpp"
 
@@ -86,6 +87,8 @@ class Input
     std::vector<Sphere> spheres() { return spheres_; }
     Sphere spheres(int i) { return spheres_[i]; }
     void spheres(const std::vector<Sphere> & sph) { spheres_ = sph; }
+    Molecule molecule() { return molecule_; }
+    void molecule(const Molecule & m) { molecule_ = m; }
     // Medium section input
     Solvent solvent() { return solvent_; }
     bool fromSolvent() { return hasSolvent_; }
@@ -168,6 +171,8 @@ class Input
     std::vector<double> radii_;
     /// List of spheres for fully custom cavity generation
     std::vector<Sphere> spheres_;
+    /// Molecule or atomic aggregate
+    Molecule molecule_;
     /// The solvent for a vacuum/uniform dielectric run
     Solvent solvent_;
     /// Whether the medium was initialized from a solvent object
 
@@ -39,12 +39,29 @@
 #include "MathUtils.hpp"
 
 Molecule::Molecule(int nat, const Eigen::VectorXd & chg, const Eigen::VectorXd & m, 
-             const Eigen::MatrixX3d & geo, const std::vector<Atom> & at, const std::vector<Sphere> & sph)
+             const Eigen::Matrix3Xd & geo, const std::vector<Atom> & at, const std::vector<Sphere> & sph)
 	: nAtoms_(nat), charges_(chg), masses_(m), geometry_(geo), atoms_(at), spheres_(sph) 
 {
     rotor_ = findRotorType();
 }
 
+Molecule::Molecule(int nat, const std::vector<Sphere> & sph)
+	: nAtoms_(nat), spheres_(sph) 
+{
+    // This constructor is used when initializing the Molecule in EXPLICIT mode
+    // charges are set to 1.0; masses are set based on the radii; geometry is set from the list of spheres
+    charges_ = Eigen::VectorXd::Ones(nAtoms_);
+    for (int i = 0; i < nAtoms_; ++i) {
+        masses_(i) = spheres_[i].radius();
+	geometry_.col(i) = spheres_[i].center();
+	double charge = charges_(i);
+	double mass = masses_(i);
+	// All the atoms are dummies
+	atoms_[i] = Atom("Dummy", "Du", charge, mass, mass, geometry_.col(i));
+    }
+    rotor_ = findRotorType();
+}
+
 
 Molecule::Molecule(const Molecule &other){
     *this = other;
@@ -65,14 +82,14 @@ Eigen::Matrix3d Molecule::inertiaTensor(){
 
     for (int i = 0; i < nAtoms_; ++i){
         // Diagonal
-        inertia(0,0) += masses_(i) * (geometry_(i,1) * geometry_(i,1) + geometry_(i,2) * geometry_(i,2));
-        inertia(1,1) += masses_(i) * (geometry_(i,0) * geometry_(i,0) + geometry_(i,2) * geometry_(i,2));
-        inertia(2,2) += masses_(i) * (geometry_(i,0) * geometry_(i,0) + geometry_(i,1) * geometry_(i,1));
+        inertia(0,0) += masses_(i) * (geometry_(1,i) * geometry_(1,i) + geometry_(2,i) * geometry_(2,i));
+        inertia(1,1) += masses_(i) * (geometry_(0,i) * geometry_(0,i) + geometry_(2,i) * geometry_(2,i));
+        inertia(2,2) += masses_(i) * (geometry_(0,i) * geometry_(0,i) + geometry_(1,i) * geometry_(1,i));
 
         // Off-diagonal
-        inertia(0,1) -= masses_(i) * (geometry_(i,0) * geometry_(i,1));
-        inertia(0,2) -= masses_(i) * (geometry_(i,0) * geometry_(i,2));
-        inertia(1,2) -= masses_(i) * (geometry_(i,1) * geometry_(i,2));
+        inertia(0,1) -= masses_(i) * (geometry_(0,i) * geometry_(1,i));
+        inertia(0,2) -= masses_(i) * (geometry_(0,i) * geometry_(2,i));
+        inertia(1,2) -= masses_(i) * (geometry_(1,i) * geometry_(2,i));
     }
     // Now symmetrize
     hermitivitize(inertia);
@@ -134,8 +151,8 @@ rotorType Molecule::findRotorType(){
 void Molecule::translate(const Eigen::Vector3d &translationVector){
     // Translate the geometry_ matrix and update the geometric data in atoms_.
     for (int i = 0; i < nAtoms_; ++i){
-        geometry_.row(i) -= translationVector;
-	Eigen::Vector3d tmp = geometry_.row(i).transpose();
+        geometry_.col(i) -= translationVector;
+	Eigen::Vector3d tmp = geometry_.col(i);
         atoms_[i].atomCoord(tmp);
     }
 }
@@ -149,7 +166,7 @@ void Molecule::rotate(const Eigen::Matrix3d &rotationMatrix){
     // Rotate the geometry_ matrix and update the geometric data in atoms_.
     geometry_ *= rotationMatrix; // The power of Eigen: geometry_ = geometry_ * rotationMatrix;
     for (int i = 0; i < nAtoms_; ++i){
-	Eigen::Vector3d tmp = geometry_.row(i).transpose();
+	Eigen::Vector3d tmp = geometry_.col(i);
         atoms_[i].atomCoord(tmp);
     }
 }
@@ -190,7 +207,7 @@ std::ostream & operator<<(std::ostream &os, const Molecule &m){
         os << "       Center              X                  Y                   Z       " << std::endl;
         os << "    ------------   -----------------  -----------------  -----------------" << std::endl;
         for (int i = 0; i < m.nAtoms_; ++i){
-            os << std::setw(10) << m.atoms_[i].atomSymbol() << std::setw(15) <<m.geometry_.row(i).format(CleanFmt) << std::endl;
+            os << std::setw(10) << m.atoms_[i].atomSymbol() << std::setw(15) <<m.geometry_.col(i).transpose().format(CleanFmt) << std::endl;
         }
     } else {
         os << "  No atoms in this molecule!" << std::endl;
 
@@ -60,7 +60,7 @@ class Molecule {
     Eigen::VectorXd masses_;
     /// Molecular geometry, in cartesian coordinates. The dimensions are (# atoms * 3)
     /// Units are Bohr.
-    Eigen::MatrixX3d geometry_;
+    Eigen::Matrix3Xd geometry_;
     /// A container for all the atoms composing the molecule
     std::vector<Atom> atoms_;
     /// A container for the spheres composing the molecule
@@ -69,17 +69,17 @@ class Molecule {
     rotorType rotor_;
 public:
     Molecule() {}
-    /// Constructor
     Molecule(int nat, const Eigen::VectorXd & chg, const Eigen::VectorXd & masses, 
-             const Eigen::MatrixX3d & geo, const std::vector<Atom> & at, const std::vector<Sphere> & sph); 
+             const Eigen::Matrix3Xd & geo, const std::vector<Atom> & at, const std::vector<Sphere> & sph); 
+    Molecule(int nat, const std::vector<Sphere> & sph); 
     /// Copy constructor.
     Molecule(const Molecule &other);
     ~Molecule(){}
 
     int nAtoms() { return nAtoms_; }
     Eigen::VectorXd charges() { return charges_; }
     Eigen::VectorXd masses() { return masses_; }
-    Eigen::MatrixX3d geometry() { return geometry_; }
+    Eigen::Matrix3Xd geometry() { return geometry_; }
     std::vector<Atom> atoms() { return atoms_; }
     std::vector<Sphere> spheres() { return spheres_; }