class MomentOfInertiaDescriptor extends AbstractMolecularDescriptor implements IMolecularDescriptor {

    @Override
    public DescriptorValue calculate(IAtomContainer container) {
        if (!GeometryUtil.has3DCoordinates(container))
            return getDummyDescriptorValue(new CDKException("Molecule must have 3D coordinates"));

        IAtomContainer clone;
        IsotopeFactory factory;
        try {
            clone = (IAtomContainer) container.clone();
            factory = Isotopes.getInstance();
            factory.configureAtoms(clone);
        } catch (CloneNotSupportedException | IOException e) {
            logger.debug(e);
            return getDummyDescriptorValue(e);
        }

        DoubleArrayResult retval = new DoubleArrayResult(7);

        double ccf = 1.000138;
        double eps = 1e-5;

        double[][] imat = new double[3][3];
        Point3d centerOfMass = GeometryUtil.get3DCentreOfMass(clone);

        double xdif;
        double ydif;
        double zdif;
        double xsq;
        double ysq;
        double zsq;
        for (int i = 0; i < clone.getAtomCount(); i++) {
            IAtom currentAtom = clone.getAtom(i);

            Double mass = factory.getMajorIsotope(currentAtom.getSymbol()).getExactMass();
            if (mass == null)
                mass = factory.getNaturalMass(currentAtom);

            xdif = currentAtom.getPoint3d().x - centerOfMass.x;
            ydif = currentAtom.getPoint3d().y - centerOfMass.y;
            zdif = currentAtom.getPoint3d().z - centerOfMass.z;
            xsq = xdif * xdif;
            ysq = ydif * ydif;
            zsq = zdif * zdif;

            imat[0][0] += mass * (ysq + zsq);
            imat[1][1] += mass * (xsq + zsq);
            imat[2][2] += mass * (xsq + ysq);

            imat[1][0] += -1 * mass * ydif * xdif;
            imat[0][1] = imat[1][0];

            imat[2][0] += -1 * mass * xdif * zdif;
            imat[0][2] = imat[2][0];

            imat[2][1] += -1 * mass * ydif * zdif;
            imat[1][2] = imat[2][1];
        }

        // diagonalize the MI tensor
        Matrix tmp = new Matrix(imat);
        EigenvalueDecomposition eigenDecomp = tmp.eig();
        double[] eval = eigenDecomp.getRealEigenvalues();

        retval.add(eval[2]);
        retval.add(eval[1]);
        retval.add(eval[0]);

        double etmp = eval[0];
        eval[0] = eval[2];
        eval[2] = etmp;

        if (Math.abs(eval[1]) > 1e-3)
            retval.add(eval[0] / eval[1]);
        else
            retval.add(1000);

        if (Math.abs(eval[2]) > 1e-3) {
            retval.add(eval[0] / eval[2]);
            retval.add(eval[1] / eval[2]);
        } else {
            retval.add(1000);
            retval.add(1000);
        }

        // finally get the radius of gyration
        double pri;
        IMolecularFormula formula = MolecularFormulaManipulator.getMolecularFormula(clone);
        if (Math.abs(eval[2]) > eps)
            pri = Math.pow(eval[0] * eval[1] * eval[2], 1.0 / 3.0);
        else
            pri = Math.sqrt(eval[0] * ccf / MolecularFormulaManipulator.getTotalExactMass(formula));
        retval.add(Math.sqrt(Math.PI * 2 * pri * ccf / MolecularFormulaManipulator.getTotalExactMass(formula)));

        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                getDescriptorNames());
    }

}