PEDRA fuzz solved. Same number of tesserae as DALTON are now obtained.

The problem was in the formation of the inertia tensor, used in rotating the cavity to standard orientation

Author: Roberto Di Remigio

src/cavity/GePolCavity.cpp

@@ -97,7 +97,10 @@ void GePolCavity::build(int maxts, int maxsph, int maxvert)
     double * ze = zv.data();
 
     double * rin = radii_scratch.data();
-    double * masses = molecule_.masses().data();
+    double * mass = new double[molecule_.nAtoms()]; 
+    for (int i = 0; i < molecule_.nAtoms(); ++i) {
+	    mass[i] = molecule_.masses(i);
+    }
 
     addedSpheres = 0;
     // Number of generators and generators of the point group
@@ -110,7 +113,7 @@ void GePolCavity::build(int maxts, int maxsph, int maxvert)
     generatecavity_cpp(&maxts, &maxsph, &maxvert,
                        xtscor, ytscor, ztscor, ar, xsphcor, ysphcor, zsphcor, rsph,
                        &nts, &ntsirr, &nSpheres_, &addedSpheres,
-                       xe, ye, ze, rin, masses,
+                       xe, ye, ze, rin, mass,
 		       &averageArea, &probeRadius, &minimalRadius,
                        &nr_gen, &gen1, &gen2, &gen3);
 
@@ -198,6 +201,7 @@ void GePolCavity::build(int maxts, int maxsph, int maxvert)
     delete[] ysphcor;
     delete[] zsphcor;
     delete[] rsph;
+    delete[] mass;
 
     built = true;
 }

src/pedra/pedra_cavity.F90

@@ -370,7 +370,7 @@ subroutine polyhedra(intsph,vert,centr,newsph,icav1,icav2,xval,yval, &
     geom = 0.0d0
 
     do i = 1, nesfp
-        mass(i)   = masses(i) !rin(i)
+        mass(i)   = masses(i)
         geom(i,1) = xe(i)
         geom(i,2) = ye(i)
         geom(i,3) = ze(i)
@@ -2275,7 +2275,8 @@ subroutine pcmtns(vmat, geom, amass, katom)
 #include "pcm_pcm.h"
 
     integer(kind=regint_k), intent(in) :: katom
-    real(kind=dp),    intent(in) :: geom(katom, 3), amass(katom)
+    real(kind=dp), intent(inout) :: geom(katom, 3)
+    real(kind=dp),    intent(in) :: amass(katom)
     real(kind=dp), intent(inout) :: vmat(3, 3)
 
     real(kind=dp) :: eigval(3), eigvec(3, 3), tinert(3, 3)
@@ -2285,28 +2286,27 @@ subroutine pcmtns(vmat, geom, amass, katom)
     integer(kind=regint_k) :: i, j, k, jax, nmax
     integer(kind=regint_k) :: nopax, nshift
     real(kind=dp) :: com(3), total_mass
-    real(kind=dp) :: local_geom(katom, 3)
-
     real(kind=dp) :: dij
 
     iax = [1, 2, 3, 3, 2, 1]
 
     angmom = [1.0d0, 1.0d0, 1.0d0]
 
     ! Calculate center-of-mass (com) and translate
-    total_mass = sum(amass)
+    total_mass = sum(amass) 
     do i = 1, 3
+      com(i) = 0.0_dp
       do j = 1, katom
         com(i) = com(i) + geom(j, i) * amass(j)
       end do
       com(i) = com(i) / total_mass
     end do
-    ! Translation
+    ! Translation to COM
     do i = 1, katom
-      local_geom(i, :) = geom(i, :) - com(:)
+      geom(i, :) = geom(i, :) - com(:)
     end do
 
-    call wlkdin(local_geom, amass, nesfp, angmom, tinert, omegad, eigval, eigvec, .true., planar, linear)
+    call wlkdin(geom, amass, nesfp, angmom, tinert, omegad, eigval, eigvec, .true., planar, linear)
 
     do i = 1, 3
         do j = 1, 3

src/utils/Molecule.hpp

@@ -78,10 +78,14 @@ class Molecule {
 
     int nAtoms() { return nAtoms_; }
     Eigen::VectorXd charges() { return charges_; }
+    double charges(int i) { return charges_(i); }
     Eigen::VectorXd masses() { return masses_; }
+    double masses(int i) { return masses_(i); }
     Eigen::Matrix3Xd geometry() { return geometry_; }
     std::vector<Atom> atoms() { return atoms_; }
+    Atom atoms(int i) const { return atoms_[i]; }
     std::vector<Sphere> spheres() const { return spheres_; }
+    Sphere spheres(int i) const { return spheres_[i]; }
 
     rotorType rotor();
     rotorType findRotorType();