ProSHADE.cpp

Go to the documentation of this file.

//============================================ Clipper
#include <clipper/clipper.h>
#include <clipper/clipper-contrib.h>
#include <clipper/clipper-ccp4.h>
#include <clipper/clipper-mmdb.h>
#include <clipper/clipper-minimol.h>

//============================================ FFTW3 + SOFT
#ifdef __cplusplus
extern "C" {
#endif
#include <fftw3.h>
#include <wrap_fftw.h>
#include <makeweights.h>
#include <s2_primitive.h>
#include <s2_cospmls.h>
#include <s2_legendreTransforms.h>
#include <s2_semi_fly.h>
#include <rotate_so3_utils.h>
#include <utils_so3.h>
#include <soft_fftw.h>
#include <rotate_so3_fftw.h>
#ifdef __cplusplus
}
#endif

//============================================ CMAPLIB
#include <cmaplib.h>

//============================================ RVAPI
#include <rvapi_interface.h>

//============================================ ProSHADE
#include "ProSHADE.h"
#include "ProSHADE_internal.h"
#include "ProSHADE_files.h"
#include "ProSHADE_misc.h"
#include "ProSHADE_legendre.h"
#include "ProSHADE_rvapi.h"

ProSHADE::ProSHADE_settings::ProSHADE_settings ( )
{
    //======================================== Settings regarding resolutions
    this->mapResolution                       = 8.0;
    this->wasResolutionGiven                  = false;
    this->bandwidth                           = 0;
    this->wasBandGiven                        = false;
    this->glIntegOrder                        = 0;
    this->theta                               = 0;
    this->phi                                 = 0;
    this->mapResDefault                       = true;
    
    //======================================== Settings regarding B factors
    this->bFactorValue                        = 80.0;
    this->bFactorChange                       = 0.0;
    this->wasBChangeGiven                     = false;
    
    //======================================== Settings regarding phase
    this->usePhase                            = true;
    this->saveWithAndWithout                  = false;
    
    //======================================== Settings regarding concentric shells
    this->shellSpacing                        = 0.0;
    this->wasShellSpacingGiven                = false;
    this->manualShells                        = 0;
    
    //======================================== Settings regarding map with phases
    this->useCOM                              = true;
    this->firstLineCOM                        = false;
    
    //======================================== Settings regarding space around the structure in lattice
    this->extraSpace                          = 8.0;
    this->wasExtraSpaceGiven                  = false;
    
    //======================================== Settings regarding weighting the distances
    this->alpha                               = 1.0;
    this->mPower                              = 1.0;
    
    //======================================== Settings regarding the map masking threshold
    this->noIQRsFromMap                       = 4.0;
    
    //======================================== Settings regarding peak finding
    this->peakHeightNoIQRs                    = 3.0;
    this->peakDistanceForReal                 = 0.2;
    this->peakSurroundingPoints               = 1;
    
    //======================================== Settings regarding which distances to compute
    this->energyLevelDist                     = true;
    this->traceSigmaDist                      = true;
    this->fullRotFnDist                       = true;
    
    //======================================== Settings regarding symmetry detection
    this->aaErrorTolerance                    = 0.1;
    this->symGapTolerance                     = 0.2;
    
    //======================================== Settings regarding hierarchical distances computation
    this->enLevelsThreshold                   = -999.9;
    this->trSigmaThreshold                    = -999.9;
    
    //======================================== Settings regarding bands to ignore
    this->ignoreLs.emplace_back ( 0 );
    
    //======================================== Settings regarding the task to perform
    this->taskToPerform                       = NA;
    
    //======================================== Settings regarding where and if to save the clear map
    this->clearMapFile                        = "";
    this->useCubicMaps                        = false;
    this->clearMapData                        = false;
    this->maskBlurFactor                      = 500.0;
    this->maskBlurFactorGiven                 = false;
    
    //======================================== Settings regarding map fragmentation
    this->mapFragBoxSize                      = 0.0;
    this->mapFragName                         = "./mapFrag";
    this->mapFragBoxFraction                  = 0.5;
    
    //======================================== Settings regarding the database
    this->databaseName                        = "";
    this->databaseMaxVolume                   = 0.0;
    this->databaseMinVolume                   = 0.0;
    this->volumeTolerance                     = 0.2;
    
    //======================================== Settings regarding the symmetry type required
    this->symmetryFold                        = 0;
    this->symmetryType                        = "";
    
    //======================================== Settings regarding the map rotation mode
    this->rotAngle                            = 0.0;
    this->rotXAxis                            = 0.0;
    this->rotYAxis                            = 0.0;
    this->rotZAxis                            = 0.0;
    this->rotChangeDefault                    = false;
    this->xTranslation                        = 0.0;
    this->yTranslation                        = 0.0;
    this->zTranslation                        = 0.0;
    
    //======================================== Settings regarding the map overlay mode
    this->overlayDefaults                     = false;
    this->dbDistOverlay                       = false;
    this->resizeMovingStructure               = false;
    this->maxRotError                         = 0;
    this->deleteModels                        = std::vector<std::string> ( );
    
    //======================================== Settings regarding the map saving mode
    this->axisOrder                           = "xyz";
    
    //======================================== Settings regarding the verbosity of the program
    this->htmlReport                          = true;
    this->htmlReportLine                      = 0;
    this->htmlReportLineProgress              = 0;
    this->verbose                             = 1;
    
}

void ProSHADE::ProSHADE_settings::printSettings ( )
{
    if ( this->verbose > 3 )
    {
        std::cout << "-----------------------------------------------------------" << std::endl;
        std::cout << "|                        SETTINGS:                        |" << std::endl;
        std::cout << "-----------------------------------------------------------" << std::endl;
        std::stringstream strstr2;
        strstr2 << this->mapResolution;
        printf (     "Resolution          : %37s\n", strstr2.str().c_str() );
        
        std::stringstream strstr3;
        strstr3 << this->bandwidth;
        printf (     "Bandwidth           : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->glIntegOrder;
        printf (     "Integration order   : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->theta;
        printf (     "Theta angle         : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->phi;
        printf (     "Phi angle           : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->bFactorValue;
        printf (     "B-factor value      : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->bFactorChange;
        printf (     "B-factor change     : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->usePhase ) { strstr3 << "TRUE"; }
        else                   { strstr3 << "FALSE"; }
        printf (     "Use phase           : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->useCOM   ) { strstr3 << "TRUE"; }
        else                   { strstr3 << "FALSE"; }
        printf (     "Use COM to center   : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->firstLineCOM ) { strstr3 << "TRUE"; }
        else                       { strstr3 << "FALSE"; }
        printf (     "Use first line COM  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->shellSpacing;
        printf (     "Shell distances in A: %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->manualShells;
        printf (     "Manual no shells    : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->extraSpace;
        printf (     "Extra cell space    : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->alpha;
        printf (     "Raise |F| to power  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->mPower;
        printf (     "Cross-cor weighting : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->noIQRsFromMap;
        printf (     "IQRs masking thres  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->peakHeightNoIQRs;
        printf (     "IQRs peak thres     : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->peakDistanceForReal;
        printf (     "Peak similarity thr : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->peakSurroundingPoints;
        printf (     "Peak width          : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->energyLevelDist ) { strstr3 << "TRUE"; }
        else                          { strstr3 << "FALSE"; }
        printf (     "Get cross-cor dists : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->traceSigmaDist  ) { strstr3 << "TRUE"; }
        else                          { strstr3 << "FALSE"; }
        printf (     "Get tr sigma dists  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->fullRotFnDist  ) { strstr3 << "TRUE"; }
        else                         { strstr3 << "FALSE"; }
        printf (     "Get rot func dists  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->aaErrorTolerance;
        printf (     "Symmetry axis toler : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->symGapTolerance;
        printf (     "Symmetry gap toler  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->enLevelsThreshold;
        printf (     "Cross-cor hierar thr: %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->trSigmaThreshold;
        printf (     "Tr sigma hierar thr : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << 0;
        for ( int i = 0; i < static_cast<int> ( this->ignoreLs.size() ); i++ )
        {
            if ( this->ignoreLs.at(i) != 0 )
            {
                strstr3 << ", " << this->ignoreLs.at(i);
            }
        }
        printf (     "Bands to ignore     : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << "\"" << this->clearMapFile << "\"";
        printf (     "Map saving location : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->useCubicMaps  ) { strstr3 << "TRUE"; }
        else                        { strstr3 << "FALSE"; }
        printf (     "Use cubic maps      : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->clearMapData  ) { strstr3 << "TRUE"; }
        else                        { strstr3 << "FALSE"; }
        printf (     "Map noice removal   : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->maskBlurFactor;
        printf (     "Map mask blurring   : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->maskBlurFactorGiven ) { strstr3 << "TRUE"; }
        else                             { strstr3 << "FALSE"; }
        printf (     "Blur. factor given  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->maxRotError;
        printf (     "Max Rotation Error  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->mapFragBoxSize;
        printf (     "Map fragm box size  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->mapFragName;
        printf (     "Map fragm save name : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->mapFragBoxFraction;
        printf (     "Map frag min density: %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << "\"" << this->databaseName << "\"";
        printf (     "Database name       : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->volumeTolerance;
        printf (     "Database volume tol : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->symmetryFold;
        printf (     "Required sym fold   : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << "\'" << this->symmetryType << "\'";
        printf (     "Required sym type   : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->rotAngle;
        printf (     "Map rotation angle  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->rotXAxis;
        printf (     "Map rotation X-axis : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->rotYAxis;
        printf (     "Map rotation Y-axis : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->rotZAxis;
        printf (     "Map rotation Z-axis : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->xTranslation;
        printf (     "Map transl.  X-axis : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->yTranslation;
        printf (     "Map transl.  Y-axis : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->zTranslation;
        printf (     "Map transl.  Z-axis : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->axisOrder;
        printf (     "Out map axis order  : %37s\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        strstr3 << this->verbose;
        printf (     "Verbose             : %37s\n\n", strstr3.str().c_str() );
        
        strstr3.str(std::string());
        if ( this->taskToPerform == Symmetry ) { strstr3 << "SYMMETRY DETECTION"; }
        if ( this->taskToPerform == Distances ) { strstr3 << "DISTANCES COMPUTATION"; }
        if ( this->taskToPerform == Features ) { strstr3 << "MAP FEATURES EXTRACTION"; }
        if ( this->taskToPerform == BuildDB ) { strstr3 << "BUILD STRUCTURE DATABASE"; }
        if ( this->taskToPerform == DistancesFrag ) { strstr3 << "FRAGMENT DISTANCES"; }
        if ( this->taskToPerform == HalfMaps ) { strstr3 << "HALF MAPS RE-BOXING"; }
        if ( this->taskToPerform == RotateMap ) { strstr3 << "MAP ROTATION"; }
        if ( this->taskToPerform == OverlayMap ) { { strstr3 << "MAP OVERLAY"; } }
        if ( this->taskToPerform == SimpleRebox ) { { strstr3 << "MAP RE-BOXING"; } }
        printf (     "Task to perform     : %37s\n", strstr3.str().c_str() );
        
        std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
    }
    
}

void ProSHADE::ProSHADE_settings::getCommandLineParams ( int argc,
                                                         char** argv )
{
    //======================================== If no command line arguments, print help
    if ( argc == 1 ) { ProSHADE_internal_misc::printHelp ( ); }
    
    //======================================== Long options struct
    const struct option longopts[] =
    {
        { "version",         no_argument,        nullptr, 'v' },
        { "help",            no_argument,        nullptr, 'h' },
        { "distances",       no_argument,        nullptr, 'D' },
        { "buildDB",         no_argument,        nullptr, 'B' },
        { "features",        no_argument,        nullptr, 'F' },
        { "symmetry",        no_argument,        nullptr, 'S' },
        { "distancesFrag",   no_argument,        nullptr, 'M' },
        { "halpMaps",        no_argument,        nullptr, 'H' },
        { "strOverlay",      no_argument,        nullptr, 'O' },
        { "reBox",           no_argument,        nullptr, 'E' },
        { "rotate",          no_argument,        nullptr, 'R' },
        { "no_phase",        no_argument,        nullptr, 'p' },
        { "no_COM",          no_argument,        nullptr, 'c' },
        { "fl_Cen",          no_argument,        nullptr, 'l' },
        { "no_crCorr",       no_argument,        nullptr, 'e' },
        { "no_trSig",        no_argument,        nullptr, 't' },
        { "no_rotFn",        no_argument,        nullptr, 'r' },
        { "clear",           no_argument,        nullptr, 'y' },
        { "cubic",           no_argument,        nullptr, 'U' },
        { "no_report",       no_argument,        nullptr, 'A' },
        { "resizeOverlay",   no_argument,        nullptr, 'C' },
        { "f",               required_argument,  nullptr, 'f' },
        { "fileList",        required_argument,  nullptr, 'i' },
        { "delModel",        required_argument,  nullptr, 'P' },
        { "bandList",        required_argument,  nullptr, 'b' },
        { "resolution",      required_argument,  nullptr, 's' },
        { "bandWidth",       required_argument,  nullptr, 'a' },
        { "integration",     required_argument,  nullptr, 'n' },
        { "maxRotErr",       required_argument,  nullptr, 'L' },
        { "sym",             required_argument,  nullptr, 'u' },
        { "clearMap",        required_argument,  nullptr, 'm' },
        { "axisOrder",       required_argument,  nullptr, 'g' },
        { "dbFile",          required_argument,  nullptr, 'x' },
        { "bValues",         required_argument,  nullptr, '{' },
        { "bChange",         required_argument,  nullptr, '}' },
        { "sphDistance",     required_argument,  nullptr, '[' },
        { "sphNumber",       required_argument,  nullptr, ']' },
        { "mapMaskThres",    required_argument,  nullptr, '(' },
        { "mapMaskBlur",     required_argument,  nullptr, 'K' },
        { "cellBorderSpace", required_argument,  nullptr, ')' },
        { "fPower",          required_argument,  nullptr, '@' },
        { "mPower",          required_argument,  nullptr, '%' },
        { "peakSize",        required_argument,  nullptr, '^' },
        { "peakThres",       required_argument,  nullptr, '&' },
        { "peakMinSim",      required_argument,  nullptr, '*' },
        { "axisTolerance",   required_argument,  nullptr, 'Q' },
        { "symGap",          required_argument,  nullptr, 'q' },
        { "CCThres",         required_argument,  nullptr, 'W' },
        { "TSThres",         required_argument,  nullptr, 'w' },
        { "mFrag",           required_argument,  nullptr, 'o' },
        { "mBoxFrac",        required_argument,  nullptr, 'G' },
        { "mBoxPath",        required_argument,  nullptr, 'T' },
        { "dbSizeLim",       required_argument,  nullptr, 'Y' },
        { "rotAng",          required_argument,  nullptr, '<' },
        { "rotX",            required_argument,  nullptr, '>' },
        { "rotY",            required_argument,  nullptr, '~' },
        { "rotZ",            required_argument,  nullptr, '|' },
        { "trsX",            required_argument,  nullptr, 'N' },
        { "trsY",            required_argument,  nullptr, 'J' },
        { "trsZ",            required_argument,  nullptr, 'I' },
        { "verbose",         optional_argument,  nullptr, 'V' },
        { nullptr,           0,                  nullptr,  0  },
    };
    
    //======================================== Change the defaults for Patterson data, unless user supplied their own values
    bool userRes                              = false;
    bool userBVal                             = false;
    
    //======================================== Short options string
    const char* const shortopts               = "vhDBFSMHOERpcletryUACf:i:P:b:s:a:n:L:u:m:g:x:{:}:[:]:(:K:):@:^:&:*:Q:q:W:w:o:G:T:Y:<:>:~:|:N:J:I:V::?";
    
    //======================================== Parsing the options
    while ( true )
    {
        //==================================== Read the next option
        const auto opt                        = getopt_long ( argc, argv, shortopts, longopts, nullptr );
        
        //==================================== Done parsing
        if ( opt == -1 )
        {
            break;
        }
        
        //==================================== For each option, set the internal values appropriately
        const char *tmp_optarg                = optarg;
        switch (opt)
        {
            //================================ Patterson should be used instead of phased maps
            case 'D':
            {
                this->taskToPerform           = Distances;
                continue;
            }
                
            //================================ Build the database from the input files
            case 'B':
            {
                this->taskToPerform           = BuildDB;
                continue;
            }
               
            //================================ Find the features of a MAP file
            case 'F':
            {
                this->taskToPerform           = Features;
                continue;
            }
                
            //================================Detect symmetry
            case 'S':
            {
                this->taskToPerform           = Symmetry;
                this->useCOM                  = false;
                this->extraSpace              = 0.0;
                this->mapResolution           = 6.0;
                continue;
            }
                
            //================================ Fragment distances to database required
            case 'M':
            {
                this->taskToPerform           = DistancesFrag;
                continue;
            }
                
            //================================ Half maps will need to be cleared
            case 'H':
            {
                this->taskToPerform           = HalfMaps;
                continue;
            }
                
            //================================ Overlay mode selected
            case 'O':
            {
                this->taskToPerform           = OverlayMap;
                this->clearMapData            = false;
                this->useCOM                  = false;
                this->mapResolution           = 2.0;
                this->rotChangeDefault        = true;
                this->extraSpace              = -777.7;
                this->overlayDefaults         = true;
                continue;
            }
                
            //================================ Simple reboxing selected
            case 'E':
            {
                this->taskToPerform           = SimpleRebox;
                this->clearMapData            = false;
                this->useCOM                  = false;
                this->mapResolution           = 2.0;
                this->rotChangeDefault        = true;
                this->extraSpace              = 3.0;
                this->overlayDefaults         = true;
                continue;
            }
                
            //================================ Rotate map
            case 'R':
            {
                this->taskToPerform           = RotateMap;
                this->clearMapData            = false;
                this->useCOM                  = false;
                this->mapResolution           = 1.0;
                this->ignoreLs                = std::vector<int> ();
                this->rotChangeDefault        = true;
                this->extraSpace              = 0.0;
                continue;
            } 
                
            //================================ Patterson should be used instead of phased maps
            case 'p':
            {
                this->usePhase                = false;
                if ( !userBVal )
                {
                    this->bFactorValue        = 120;
                }
                if ( !userRes )
                {
                    this->mapResolution       = 9.0;
                }
                continue;
            }

            //================================ Centreing has already been done
            case 'c':
            {
                this->useCOM                  = false;
                continue;
            }
                
            //================================ Centre on the first line coordinates
            case 'l':
            {
                this->firstLineCOM            = true;
                continue;
            }
            
            //================================ No interest in cross-correlation distance
            case 'e':
            {
                this->energyLevelDist         = false;
                continue;
            }
                
            //================================ No interest in trace sigma distance
            case 't':
            {
                this->traceSigmaDist          = false;
                continue;
            }
                
            //================================ No interest in rotation function distance
            case 'r':
            {
                this->fullRotFnDist           = false;
                continue;
            }
                
            //================================ No interest in rotation function distance
            case 'y':
            {
                this->clearMapData            = true;
                continue;
            }
                
            //================================ Do use the cubic clear maps instead of rectangular
            case 'U':
            {
                this->useCubicMaps            = true;
                continue;
            }
                
            //================================ Do not produce the HTML report
            case 'A':
            {
                this->htmlReport              = false;
                continue;
            }
                
            //================================ Do not resize the moving molecule to have same dimensions as the static molecule
            case 'C':
            {
                this->resizeMovingStructure   = true;
                continue;
            }
                
            //================================ File name of input file
            case 'f':
            {
                this->structFiles.emplace_back ( std::string ( optarg ) );
                continue;
            }
                
            //================================ File list of input files
            case 'i':
            {
                //============================ Open the file
                std::ifstream textFile ( std::string ( optarg ), std::ios::in );
                
                //============================ Check
                if ( textFile.fail() )
                {
                    std::cout << "!!! ProSHADE ERROR !!! The filename supplied with the \'-i\' parameter cannot be opened. Terminating..." << std::endl;
                    exit ( -1 );
                }
                
                //============================ Read
                std::string hlpStr;
                while ( std::getline ( textFile, hlpStr ) )
                {
                    this->structFiles.emplace_back ( hlpStr );
                }
                
                //============================ Close
                textFile.close ();
                continue;
            }
                
            //================================ File name of model to be deleted
            case 'P':
            {
                this->deleteModels.emplace_back ( std::string ( optarg ) );
                continue;
            }
                
            //================================ List of bands to ignore
            case 'b':
            {
                this->ignoreLs.emplace_back ( stoi ( std::string ( optarg ) ) );
                continue;
            }
                
            //================================ List of bands to ignore
            case 's':
            {
                this->mapResolution           = stof ( std::string ( optarg ) );
                userRes                       = true;
                this->mapResDefault           = false;
                this->wasResolutionGiven      = true;
                continue;
            }
                
            //================================ Maximum bandwidth
            case 'a':
            {
                this->bandwidth               = stof ( std::string ( optarg ) );
                this->wasBandGiven            = true;
                continue;
            }
                
            //================================ Maximum integration order
            case 'n':
            {
                this->glIntegOrder            = stoi ( std::string ( optarg ) );
                continue;
            }
             
            //================================ Maximum rotation error
            case 'L':
            {
                this->maxRotError             = static_cast<unsigned int> ( stoi ( std::string ( optarg ) ) );
                continue;
            }
                
            //================================ Symmetry required
            case 'u':
            {
                std::string hlpStr            = std::string ( optarg );
                if ( hlpStr.at(0) == 'C' )
                {
                    this->symmetryType        = "C";
                }
                else
                {
                    if ( hlpStr.at(0) == 'D' )
                    {
                        this->symmetryType    = "D";
                    }
                    else
                    {
                        if ( hlpStr.at(0) == 'T' )
                        {
                            this->symmetryType = "T";
                        }
                        else
                        {
                            if ( hlpStr.at(0) == 'O' )
                            {
                                this->symmetryType = "O";
                            }
                            else
                            {
                                if ( hlpStr.at(0) == 'I' )
                                {
                                    this->symmetryType = "I";
                                }
                                else
                                {
                                    ProSHADE_internal_misc::printHelp ( );
                                    exit ( -1 );
                                }
                            }
                        }
                    }
                }
                
                if ( ( this->symmetryType == "C" ) || ( this->symmetryType == "D" ) )
                {
                    std::string numHlp ( hlpStr.begin()+1, hlpStr.end() );
                    if ( numHlp.length() > 0 )
                    {
                        this->symmetryFold    = stoi ( numHlp );
                    }
                    else
                    {
                        std::cerr << "!!! ProSHADE ERROR !!! The input argument requests search for Cyclic/Dihedral symmetry, but does not specify the requested fold." << std::endl;
                        exit                  ( -1 );
                    }
                }
                else
                {
                    this->symmetryFold        = 0;
                    std::string numHlp ( hlpStr.begin()+1, hlpStr.end() );
                    if ( numHlp != "" )
                    {
                        ProSHADE_internal_misc::printHelp ( );
                        exit ( -1 );
                    }
                }
                
                continue;
            }
                
            //================================ Clear map location
            case 'm':
            {
                this->clearMapFile            = std::string ( optarg );
                continue;
            }
                
            //================================ Map axes order
            case 'g':
            {
                this->axisOrder               = std::string ( optarg );
                std::transform                ( this->axisOrder.begin(),
                                                this->axisOrder.end(),
                                                this->axisOrder.begin(),
                                                ::tolower);
                if ( this->axisOrder.length() != 3 )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string of three characters only." << std::endl;
                    exit ( -1 );
                }
                if ( ( this->axisOrder[0] != 'x' ) && ( this->axisOrder[0] != 'y' ) && ( this->axisOrder[0] != 'z' )  )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string consisting of only \'x\', \'y\' and \'z\' characters." << std::endl;
                    exit ( -1 );
                }
                if ( ( this->axisOrder[1] != 'x' ) && ( this->axisOrder[1] != 'y' ) && ( this->axisOrder[1] != 'z' )  )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string consisting of only \'x\', \'y\' and \'z\' characters." << std::endl;
                    exit ( -1 );
                }
                if ( ( this->axisOrder[2] != 'x' ) && ( this->axisOrder[2] != 'y' ) && ( this->axisOrder[2] != 'z' )  )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string consisting of only \'x\', \'y\' and \'z\' characters." << std::endl;
                    exit ( -1 );
                }
                if ( ( this->axisOrder[0] != 'x' ) && ( this->axisOrder[1] != 'x' ) && ( this->axisOrder[2] != 'x' )  )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string with at least one \'x\' character." << std::endl;
                    exit ( -1 );
                }
                if ( ( this->axisOrder[0] != 'y' ) && ( this->axisOrder[1] != 'y' ) && ( this->axisOrder[2] != 'y' )  )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string with at least one \'y\' character." << std::endl;
                    exit ( -1 );
                }
                if ( ( this->axisOrder[0] != 'z' ) && ( this->axisOrder[1] != 'z' ) && ( this->axisOrder[2] != 'z' )  )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Input parameter error !!! The axisOrder parameter requires a string with at least one \'z\' character." << std::endl;
                    exit ( -1 );
                }
                continue;
            }
                
            //================================ Database file location
            case 'x':
            {
                this->databaseName            = std::string ( optarg );
                continue;
            }
                
            //================================ B factor value to which all atoms are changed to
            case '{':
            {
                this->bFactorValue            = stof ( std::string ( optarg ) );
                userBVal                      = true;
                continue;
            }
                
            //================================ B factor value change for maps (sharpening/blurring)
            case '}':
            {
                this->bFactorChange           = stof ( std::string ( optarg ) );
                this->wasBChangeGiven         = true;
                continue;
            }
                
            //================================ Distance between concentric spheres
            case '[':
            {
                this->shellSpacing            = stof ( std::string ( optarg ) );
                this->wasShellSpacingGiven    = true;
                continue;
            }
                
            //================================ Maximum number of concentric spheres
            case ']':
            {
                this->manualShells            = stoi ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Number of inter-qusrtile ranges from the third quartile to use as threshold for masking
            case '(':
            {
                this->noIQRsFromMap           = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Amount of blurring to apply to create a mask
            case 'K':
            {
                this->maskBlurFactor          = stof ( std::string ( optarg ) );
                this->maskBlurFactorGiven     = true;
                continue;
            }
                
            //================================ Extra cellular space
            case ')':
            {
                this->extraSpace              = stof ( std::string ( optarg ) );
                this->wasExtraSpaceGiven      = true;
                continue;
            }
                
            //================================ Fourier coefficient power
            case '@':
            {
                this->alpha                   = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Cross-correlation descriptor weight
            case '%':
            {
                this->mPower                  = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Peak size
            case '^':
            {
                this->peakSurroundingPoints   = stoi ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Peak threshold IQRs
            case '&':
            {
                this->peakHeightNoIQRs        = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Minimum missing peak similarity to become real
            case '*':
            {
                this->peakDistanceForReal     = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Axis tolerance
            case 'Q':
            {
                this->aaErrorTolerance        = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Symmetry gap
            case 'q':
            {
                this->symGapTolerance         = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Cross-correlation hierarchical threshold
            case 'W':
            {
                this->enLevelsThreshold       = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Trace sigma hierarchical threshold
            case 'w':
            {
               this->trSigmaThreshold         = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Map fragmentation box size
            case 'o':
            {
                this->mapFragBoxSize          = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Map box fraction full
            case 'G':
            {
                this->mapFragBoxFraction      = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Map box save path
            case 'T':
            {
                this->mapFragName             = std::string ( optarg );
                continue;
            }
                
            //================================ Tolerance on volume for struct against db
            case 'Y':
            {
                this->volumeTolerance         = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ The angle by which to rotate the structure
            case '<':
            {
                this->rotAngle                = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ The X-axis of the rotation axis vector
            case '>':
            {
                this->rotXAxis                = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================The Y-axis of the rotation axis vector
            case '~':
            {
                this->rotYAxis                = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ The Z-axis of the rotation axis vector
            case '|':
            {
                this->rotZAxis                = stof ( std::string ( optarg ) );
                continue;
            }

            //================================ The X-axis of the translation vector
            case 'N':
            {
                this->xTranslation            = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ The Y-axis of the translation vector
            case 'J':
            {
                this->yTranslation            = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ The Z-axis of the translation vector
            case 'I':
            {
                this->zTranslation            = stof ( std::string ( optarg ) );
                continue;
            }
                
            //================================ Print version info
            case 'V':
            {
                if ( !optarg && argv[optind] != NULL &&  argv[optind][0] != '-' )
                {
                    tmp_optarg                = argv[optind++];
                }
                
                if ( tmp_optarg )
                {
                    this->verbose             = stoi ( std::string ( tmp_optarg ) );
                }
                else
                {
                    this->verbose             = 3;
                }
                
                if ( ( this->verbose < -1 ) || ( this->verbose > 4 ) )
                {
                    printf ( "!!! ProSHADE ERROR !!! The \'--verbose\' option takes only values in range 0 to 4, including. Terminating...\n" );
                    exit ( -1 );
                }
                    
                continue;
            }
                
            //================================ Print version info
            case 'v':
            {
                std::cout << "ProSHADE " << __PROSHADE_VERSION__ << std::endl << std::endl;
                exit ( 0 );
            }
                
            //================================ User needs help
            case 'h':
            {
                ProSHADE_internal_misc::printHelp ( );
                exit ( 0 );
            }
                
            //================================ Unknown option
            case '?':
            {
                //============================ This case is handled by getopt_long, nothing more needed.
                exit ( 0 );
            }
                
            //================================ Fallback option
            default:
            {
                ProSHADE_internal_misc::printHelp ( );
                exit ( 0 );
            }
        }
    }
    
    //======================================== Done
    return;
    
}

void ProSHADE::ProSHADE_settings::ignoreLsAddValuePy ( const int val )
{
    //======================================== Add the value
    this->ignoreLs.emplace_back ( val );
    
    //======================================== Done
    return;
}

ProSHADE::ProSHADE::ProSHADE ( ProSHADE_settings* setUp )
{
    //======================================== Initialise variables
    this->settings                            = setUp;
    
    this->distancesAvailable                  = false;
    this->symmetriesAvailable                 = false;
    this->fragmentsAvailable                  = false;
    
    //======================================== Report progress
    if ( setUp->verbose > 0 )
    {
        std::cout << "ProSHADE " << __PROSHADE_VERSION__ << " :" << std::endl;
        std::cout << "==========================" << std::endl << std::endl;
        
        setUp->printSettings ();
    }
    
    //======================================== If no task, report error
    if ( setUp->taskToPerform == NA )
    {
        std::cerr << "!!! ProSHADE ERROR !!! There is no task (functionality) selected for this run." << std::endl << std::endl;
        std::cerr << "Please select a functionality and supply it to the settings object (or using command line). Your options are:" << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -D or --distances                                                       " << std::endl;
        std::cerr << "             The shape distances will be computed between the first supplied " << std::endl;
        std::cerr << "             structure and all other structures. Requires at least two       " << std::endl;
        std::cerr << "             structures.                                                     " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -B or --buildDB                                                         " << std::endl;
        std::cerr << "             Pre-compute the spherical harmonics coefficients for all        " << std::endl;
        std::cerr << "             supplied structures for the given options. Save the results in  " << std::endl;
        std::cerr << "             a binary database file given by the \'x\' option.                 " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -F or --features                                                        " << std::endl;
        std::cerr << "             Return a table of statistics for input map. To save a clean and " << std::endl;
        std::cerr << "             resized map, use the \'m\' option to specify the file to save to. " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -S or --symmetry                                                        " << std::endl;
        std::cerr << "             Detect if any C, D, T or I symmetries are present in the first  " << std::endl;
        std::cerr << "             structure supplied. Also prints all symmetry elements for C and " << std::endl;
        std::cerr << "             D symmetries, but only two elements for T and I.                " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -M or --distancesFrag                                                   " << std::endl;
        std::cerr << "             Fragment the supplied map file and search each resulting frag-  " << std::endl;
        std::cerr << "             ment against the whole database supplied by the \'x\' option. Per " << std::endl;
        std::cerr << "             fragment distances are reported.                                " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -R or --rotate                                                          " << std::endl;
        std::cerr << "             Take a single density map and apply the rotation defined by the " << std::endl;
        std::cerr << "             \'--rotAng\', \'--rotAlp\', \'--rotBet\' and \'--rotGam\' options. The  " << std::endl;
        std::cerr << "             resulting rotated map is saved to location given by the         " << std::endl;
        std::cerr << "             \'--clearMap\' option. Many of the default settings are over-     " << std::endl;
        std::cerr << "             written and while can be modified using the command line        " << std::endl;
        std::cerr << "             options, this is not recommended.                               " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -O or --strOverlay                                                      " << std::endl;
        std::cerr << "             Given two structures, find the optimal overlay using the        " << std::endl;
        std::cerr << "             rotation and translation functions. The first structure is      " << std::endl;
        std::cerr << "             always unchanged, while a rotated and translated version of the " << std::endl;
        std::cerr << "             second structure will be written to the \'--clearMap\' option     " << std::endl;
        std::cerr << "             path or \'./rotStr\' file.                                        " << std::endl;
        std::cerr << "                                                                             " << std::endl;
        std::cerr << "     -E or --reBox                                                           " << std::endl;
        std::cerr << "             Takes a single map structure and computes the mask. Then, finds " << std::endl;
        std::cerr << "             the box around the mask, adds extra \'--cellBorderSpace\' angs.   " << std::endl;
        std::cerr << "             of space and outputs the re-sized map to \'--clearMap\'.          " << std::endl;
        std::cerr << "                                                                             " << std::endl;
    }

    if ( setUp->htmlReport )
    {
        //==================================== Create report directory
        mkdir                                 ( "proshade_report", 0777 );
        
        //==================================== Initialise empty HTML document
        std::stringstream hlpSS;
        hlpSS << __PROSHADE_RVAPI__ << "/jsrview";
        rvapi_init_document                   ( "ProSHADE Results Report",
                                                "./proshade_report",
                                                "ProSHADE Results Report",
                                                RVAPI_MODE_Html,
                                                RVAPI_LAYOUT_Plain,
                                                hlpSS.str().c_str(),
                                                NULL,
                                                "index.html",
                                                NULL,
                                                NULL );
        rvapi_flush                           ( );
    }
    
    //======================================== If overlaying two maps, proceed
    if ( setUp->taskToPerform == OverlayMap )
    {
        if ( setUp->verbose > 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                   MODE: Structure Overlay               |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Structure Overlay</h1>",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1 );
            setUp->htmlReportLine            += 1;
            
            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            rvapi_flush                       ( );
        }
        
        if ( setUp->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Overlay of structures started." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
        
        ProSHADE_internal::ProSHADE_symmetry symmetry ( setUp );
        
        //==================================== Done

    }
    
    //======================================== If rotating map, do so
    if ( setUp->taskToPerform == RotateMap )
    {
        if ( setUp->verbose > 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                      MODE: Map Rotation                 |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }

        ProSHADE_internal::ProSHADE_symmetry symmetry ( setUp );
    }
    
    //======================================== If creating database, do so
    if ( setUp->taskToPerform == BuildDB )
    {
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Building a database of structures</h1>",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1 );
            setUp->htmlReportLine            += 1;
            
            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Database building procedure started." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            
            rvapi_flush                       ( );
        }
        
        ProSHADE_internal::ProSHADE_distances *db = new ProSHADE_internal::ProSHADE_distances ( setUp );
        delete db;
        
        if ( setUp->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "COMPLETED." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            
            rvapi_flush                       ( );
        }
    }
    
    //======================================== If computing features, do so
    if ( setUp->taskToPerform == Features )
    {
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Density map features</h1>",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1 );
            setUp->htmlReportLine            += 1;
            
            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Density map features detection procedure started." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            
            rvapi_flush                       ( );
        }
        
        ProSHADE_internal::ProSHADE_mapFeatures features ( setUp );
        
        features.printInfo                    ( setUp->verbose );
        
        if ( setUp->clearMapFile != "" )
        {
            if ( setUp->verbose > 0 )
            {
                std::cout << "-----------------------------------------------------------" << std::endl;
                std::cout << "|                      SAVING CLEAN MAP                   |" << std::endl;
                std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
            }
            
            if ( setUp->verbose > 0 )
            {
                std::cout << "Clean map saved to file: " << setUp->clearMapFile << std::endl << std::endl;
            }
        }
        
        if ( setUp->mapFragBoxSize > 0.0 )
        {
            features.fragmentMap              ( setUp->axisOrder,
                                                setUp->verbose,
                                                setUp );
            this->fragmentList                = features.getFragmentsList ( );
            this->fragmentsAvailable          = true;
        }
        else
        {
            if ( setUp->verbose > 0 )
            {
                std::cout << "-----------------------------------------------------------" << std::endl;
                std::cout << "|                      MAP FRAGMENTATION                  |" << std::endl;
                std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
                
                std::cout << "Map fragmentation was not required." << std::endl << std::endl;
            }
        }
    }
    
    //======================================== If computing symmetry, do so
    if ( setUp->taskToPerform == Symmetry )
    {        
        ProSHADE_internal::ProSHADE_symmetry symmetry ( setUp );
        this->cyclicSymmetries                = symmetry.cnSymm;
        this->cyclicSymmetriesClear           = symmetry.cnSymmClear;
        this->dihedralSymmetries              = symmetry.dnSymm;
        this->dihedralSymmetriesClear         = symmetry.dnSymmClear;
        this->tetrahedralSymmetry             = symmetry.tetrAxes;
        this->octahedralSymmetry              = symmetry.octaAxes;
        this->icosahedralSymmetry             = symmetry.icosAxes;
        this->symmetriesAvailable             = true;
        
        //==================================== Print results to standard output, if need be
        if ( setUp->symmetryType == "" )
        {
            if ( setUp->verbose > -1 )
            {
                symmetry.printResultsClear    ( setUp->verbose );
            }
            if ( setUp->htmlReport )
            {
                symmetry.printResultsClearHTML( setUp );
            }
        }
        else
        {
            if ( setUp->verbose > -1 )
            {
                symmetry.printResultsRequest  ( setUp->symmetryType,
                                                setUp->symmetryFold,
                                                setUp->verbose );
            }
            if ( setUp->htmlReport )
            {
                symmetry.printResultsRequestHTML ( setUp->symmetryType,
                                                   setUp->symmetryFold,
                                                   setUp );
            }
        }
        
        //==================================== Align symmetry axes to coordinate axes, if required
        if ( setUp->clearMapFile != "" )
        {
            if ( symmetry.inputStructureDataType )
            {
                std::cerr << "!!! ProSHADE ERROR !!! The symmetry output map can only be used for map input, the functionality to align PDB symmetry axes (which are detected) to the system axes is not yet implemented. Please report this if you are interested in a swift solution." << std::endl;
                exit ( -1 );
            }
            
            if ( setUp->verbose != -1 )
            {
                std::cout << "-----------------------------------------------------------" << std::endl;
                std::cout << "|                 Structure Axis Alignment                |" << std::endl;
                std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
            }
            
            if ( setUp->htmlReport )
            {
                std::stringstream hlpSS;
                hlpSS << "<font color=\"green\">" << "Symmetry axes alignment started." << "</font>";
                rvapi_set_text                ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                setUp->htmlReportLineProgress+= 1;
                rvapi_flush                           ( );
            }
            
            if ( static_cast<unsigned int> ( this->icosahedralSymmetry.size() ) > 0 )
            {
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Applying the ICOSAHEDRAL axis alignment convention of Heymann, Chagoyen and Belnap (2005) I1 (2-fold axes on X, Y and Z)." << std::endl;
                }
                
                if ( setUp->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"green\">" << "Applying axes alignment using the convention of Heymann, Chagoyen and Belnap (2005) <b>I1</b> (2-fold axes on X, Y and Z)." << "</font>";
                    rvapi_set_text            ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                    setUp->htmlReportLineProgress     += 1;
                    rvapi_flush                           ( );
                }
                
                //============================ Take the highest symmetry and consequently highest peak axes
                std::sort                     ( this->icosahedralSymmetry.begin(), this->icosahedralSymmetry.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[0] > b[0]; });
                
                double maxFold                = this->icosahedralSymmetry.at(0)[0];
                int noMaxFold                 = 0;
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->icosahedralSymmetry.size() ); iter++ )
                {
                    if ( this->icosahedralSymmetry.at(iter)[0] == maxFold )
                    {
                        noMaxFold            += 1;
                    }
                }
                std::sort                     ( this->icosahedralSymmetry.begin(), this->icosahedralSymmetry.begin() + noMaxFold, [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
                
                maxFold                       = this->icosahedralSymmetry.at(noMaxFold)[0];
                int noMaxFold2                = 0;
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->icosahedralSymmetry.size() ); iter++ )
                {
                    if ( this->icosahedralSymmetry.at(iter)[0] == maxFold )
                    {
                        noMaxFold2           += 1;
                    }
                }
                
                std::sort                     ( this->icosahedralSymmetry.begin() + noMaxFold, this->icosahedralSymmetry.begin() + noMaxFold + noMaxFold2, [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
                
                //============================ Initialise axes
                std::vector<std::vector<double>> origAxis;
                std::vector<std::vector<double>> basisAxis;
                std::vector<double> hlpVec;
                hlpVec.emplace_back ( 0.0 ); hlpVec.emplace_back ( 0.0 );
                origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec );
                
                //============================ Find the three (two) perpendicular C2 axes
                int c2It                      = 0;
                while ( static_cast<int> ( noMaxFold + noMaxFold2 + c2It ) < static_cast<int> ( this->icosahedralSymmetry.size() ) )
                {
                    bool firstAxis            = true;
                    for ( unsigned int iter = noMaxFold + noMaxFold2 + c2It; iter < static_cast<unsigned int> ( this->icosahedralSymmetry.size() ); iter++ )
                    {
                        if ( firstAxis )
                        {
                            firstAxis         = false;
                            origAxis.at(0).at(0) = this->icosahedralSymmetry.at(iter)[1];
                            origAxis.at(1).at(0) = this->icosahedralSymmetry.at(iter)[2];
                            origAxis.at(2).at(0) = this->icosahedralSymmetry.at(iter)[3];
                        }
                        else
                        {
                            double dotProd    = origAxis.at(0).at(0) * this->icosahedralSymmetry.at(iter)[1] +
                                                origAxis.at(1).at(0) * this->icosahedralSymmetry.at(iter)[2] +
                                                origAxis.at(2).at(0) * this->icosahedralSymmetry.at(iter)[3];
                            
                            if ( ( dotProd < 0.05 ) && ( dotProd > -0.05 ) )
                            {
                                origAxis.at(0).at(1) = this->icosahedralSymmetry.at(iter)[1];
                                origAxis.at(1).at(1) = this->icosahedralSymmetry.at(iter)[2];
                                origAxis.at(2).at(1) = this->icosahedralSymmetry.at(iter)[3];
                                
                                break;
                            }
                        }
                    }
                    
                    if ( ( std::abs( origAxis.at(0).at(1) ) < 0.05 ) && ( std::abs( origAxis.at(1).at(1) ) < 0.05 ) && ( std::abs( origAxis.at(2).at(1) ) < 0.05 ) )
                    {
                        c2It             += 1;
                    }
                    else
                    {
                        break;
                    }
                }
                
                if ( ( std::abs( origAxis.at(0).at(1) ) < 0.05 ) && ( std::abs( origAxis.at(1).at(1) ) < 0.05 ) && ( std::abs( origAxis.at(2).at(1) ) < 0.05 ) )
                {
                    std::cerr << "!!! ProSHADE ERROR !!! Cannot find three perpendicular C2 axes in the Icosahedral symmetry. This looks like a bug, please report this case. Also, you can try increasing the resolution to overcome this problem." << std::endl;
                    if ( setUp->htmlReport )
                    {
                        std::stringstream hlpSS;
                        hlpSS << "<font color=\"red\">" << "Cannot find three perpendicular C2 axes in the Icosahedral symmetry. This looks like a bug, please report this case. Also, you can try increasing the resolution to overcome this problem." << "</font>";
                        rvapi_set_text        ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                        setUp->htmlReportLineProgress += 1;
                        rvapi_flush           ( );
                    }
                    exit ( -1 );
                }
                
                //============================ Basis vectors for alignment
                hlpVec.emplace_back ( 0.0 );
                basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec );
                
                basisAxis.at(0).at(0)         = 1.0;
                basisAxis.at(1).at(0)         = 0.0;
                basisAxis.at(2).at(0)         = 0.0;
                
                basisAxis.at(0).at(1)         = 0.0;
                basisAxis.at(1).at(1)         = 1.0;
                basisAxis.at(2).at(1)         = 0.0;
                
                basisAxis.at(0).at(2)         = 0.0;
                basisAxis.at(1).at(2)         = 0.0;
                basisAxis.at(2).at(2)         = 1.0;
                
                //============================ Get angle axis representation of the matrix
                std::array<double,5> joinAA   = ProSHADE_internal_misc::getAxisAngleFromSymmetryAxes ( origAxis, basisAxis, 0, 1 );
                
                //============================ Set up rotation ProSHADE run
                ProSHADE_settings* rotSet     = new ProSHADE_settings ( );
                
                // ... Settings regarding resolutions
                rotSet->mapResolution         = setUp->mapResolution;
                
                // ... Settings regarding space around the structure in lattice
                rotSet->extraSpace            = setUp->extraSpace;
                
                // ... Settings regarding rotation module, these should to be set this way!
                rotSet->clearMapData          = false;
                rotSet->useCOM                = false;
                rotSet->rotChangeDefault      = true;
                rotSet->ignoreLs              = std::vector<int> ();
                rotSet->maskBlurFactor        = 500.0;
                rotSet->maskBlurFactorGiven   = false;
                
                // ... Settings regarding the task
                rotSet->taskToPerform         = RotateMap;
                
                // ... Settings regarding where and if to save the clear map
                rotSet->clearMapFile          = setUp->clearMapFile;
                
                // ... Settings regarding the map rotation mode
                rotSet->rotXAxis              =  joinAA[1];
                rotSet->rotYAxis              =  joinAA[2];
                rotSet->rotZAxis              =  joinAA[3];
                rotSet->rotAngle              =  joinAA[0];
                
                // ... Settings regarding the map saving mode
                rotSet->axisOrder             = setUp->axisOrder;
                
                // ... Settings regarding the standard output amount
                rotSet->verbose               = -1;
                
                // ... Settings regarding the file to rotate
                rotSet->structFiles.emplace_back ( setUp->structFiles.at(0) );
                
                //============================ Rotate the structure as required
                if ( setUp->verbose > 0 )
                {
                    std::cout << ">> Structure rotation initiated." << std::endl;
                }
                
                ProSHADE *run                 = new ProSHADE ( rotSet );
                delete run;
                
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Structure with symmetry axis to box axes alignment writted in " << setUp->clearMapFile << std::endl;
                }
            }
            
            else if ( static_cast<unsigned int> ( this->octahedralSymmetry.size() ) > 0 )
            {
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Applying the OCTAHEDRAL axis alignment convention of Heymann, Chagoyen and Belnap (2005) (4-fold axes on X, Y and Z)." << std::endl;
                }
                
                if ( setUp->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"green\">" << "Applying axes alignment using the convention of Heymann, Chagoyen and Belnap (2005) (4-fold axes on X, Y and Z)." << "</font>";
                    rvapi_set_text            ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                    setUp->htmlReportLineProgress += 1;
                    rvapi_flush               ( );
                }
                
                //============================ Take the highest symmetry and consequently highest peak axes
                std::sort                     ( this->octahedralSymmetry.begin(), this->octahedralSymmetry.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[0] > b[0]; });
                
                double maxFold                = this->octahedralSymmetry.at(0)[0];
                int noMaxFold                 = 0;
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->octahedralSymmetry.size() ); iter++ )
                {
                    if ( this->octahedralSymmetry.at(iter)[0] == maxFold )
                    {
                        noMaxFold            += 1;
                    }
                }
                std::sort                     ( this->octahedralSymmetry.begin(), this->octahedralSymmetry.begin() + noMaxFold, [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
                
                maxFold                       = this->octahedralSymmetry.at(noMaxFold)[0];
                int noMaxFold2                = 0;
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->octahedralSymmetry.size() ); iter++ )
                {
                    if ( this->octahedralSymmetry.at(iter)[0] == maxFold )
                    {
                        noMaxFold2           += 1;
                    }
                }
                std::sort                     ( this->octahedralSymmetry.begin() + noMaxFold, this->octahedralSymmetry.begin() + noMaxFold2 + noMaxFold, [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
                
                //============================ Original symmetry axes for Z-alignment
                std::vector<std::vector<double>> origAxis;
                std::vector<std::vector<double>> basisAxis;
                std::vector<double> hlpVec;
                hlpVec.emplace_back ( 0.0 ); hlpVec.emplace_back ( 0.0 );
                origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec );
                
                origAxis.at(0).at(0)          = this->octahedralSymmetry.at(0)[1];
                origAxis.at(1).at(0)          = this->octahedralSymmetry.at(0)[2];
                origAxis.at(2).at(0)          = this->octahedralSymmetry.at(0)[3];
                
                origAxis.at(0).at(1)          = this->octahedralSymmetry.at(1)[1];
                origAxis.at(1).at(1)          = this->octahedralSymmetry.at(1)[2];
                origAxis.at(2).at(1)          = this->octahedralSymmetry.at(1)[3];
                
                //============================ Basis vectors for alignment
                hlpVec.emplace_back ( 0.0 );
                basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec );
                
                basisAxis.at(0).at(0)         = 1.0;
                basisAxis.at(1).at(0)         = 0.0;
                basisAxis.at(2).at(0)         = 0.0;
                
                basisAxis.at(0).at(1)         = 0.0;
                basisAxis.at(1).at(1)         = 1.0;
                basisAxis.at(2).at(1)         = 0.0;
                
                basisAxis.at(0).at(2)         = 0.0;
                basisAxis.at(1).at(2)         = 0.0;
                basisAxis.at(2).at(2)         = 1.0;
                
                //============================ Get angle axis representation of the matrix
                std::array<double,5> joinAA   = ProSHADE_internal_misc::getAxisAngleFromSymmetryAxes ( origAxis, basisAxis, 1, 2 );
                
                //============================ Set up rotation ProSHADE run
                ProSHADE_settings* rotSet     = new ProSHADE_settings ( );
                
                // ... Settings regarding resolutions
                rotSet->mapResolution         = setUp->mapResolution;
                
                // ... Settings regarding space around the structure in lattice
                rotSet->extraSpace            = setUp->extraSpace;
                
                // ... Settings regarding rotation module, these should to be set this way!
                rotSet->clearMapData          = false;
                rotSet->useCOM                = false;
                rotSet->rotChangeDefault      = true;
                rotSet->ignoreLs              = std::vector<int> ();
                rotSet->maskBlurFactor        = 500.0;
                rotSet->maskBlurFactorGiven   = false;
                
                // ... Settings regarding the task
                rotSet->taskToPerform         = RotateMap;
                
                // ... Settings regarding where and if to save the clear map
                rotSet->clearMapFile          = setUp->clearMapFile;
                
                // ... Settings regarding the map rotation mode
                rotSet->rotXAxis              =  joinAA[1];
                rotSet->rotYAxis              =  joinAA[2];
                rotSet->rotZAxis              =  joinAA[3];
                rotSet->rotAngle              =  joinAA[0];
                
                // ... Settings regarding the map saving mode
                rotSet->axisOrder             = setUp->axisOrder;
                
                // ... Settings regarding the standard output amount
                rotSet->verbose               = -1;
                
                // ... Settings regarding the file to rotate
                rotSet->structFiles.emplace_back ( setUp->structFiles.at(0) );
                
                //============================ Rotate the structure as required
                if ( setUp->verbose > 0 )
                {
                    std::cout << ">> Structure rotation initiated." << std::endl;
                }
                
                ProSHADE *run                 = new ProSHADE ( rotSet );
                delete run;
                
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Structure with symmetry axis to box axes alignment writted in " << setUp->clearMapFile << std::endl;
                }
            }
            
            else if ( static_cast<unsigned int> ( this->tetrahedralSymmetry.size() ) > 0 )
            {
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Applying the TETRAHEDRAL axis alignment convention from RELION (3-fold axis on Z)." << std::endl;
                }
                
                if ( setUp->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"green\">" << "Applying axis alignment using the RELION convention (3-fold axis on Z)." << "</font>";
                    rvapi_set_text            ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                    setUp->htmlReportLineProgress += 1;
                    rvapi_flush               ( );
                }
                
                //============================ Take the highest symmetry and consequently highest peak axes
                std::sort                     ( this->tetrahedralSymmetry.begin(), this->tetrahedralSymmetry.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[0] > b[0]; });
                
                double maxFold                = this->tetrahedralSymmetry.at(0)[0];
                int noMaxFold                 = 0;
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->tetrahedralSymmetry.size() ); iter++ )
                {
                    if ( this->tetrahedralSymmetry.at(iter)[0] == maxFold )
                    {
                        noMaxFold            += 1;
                    }
                }
                std::sort                     ( this->tetrahedralSymmetry.begin(), this->tetrahedralSymmetry.begin() + noMaxFold, [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
                std::sort                     ( this->tetrahedralSymmetry.begin() + noMaxFold, this->tetrahedralSymmetry.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[4] > b[4]; });
                
                //============================ Original symmetry axes for Z-alignment
                std::vector<std::vector<double>> origAxis;
                std::vector<std::vector<double>> basisAxis;
                std::vector<double> hlpVec;
                hlpVec.emplace_back ( 0.0 ); hlpVec.emplace_back ( 0.0 ); 
                origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec );
                
                origAxis.at(0).at(0)          = this->tetrahedralSymmetry.at(0)[1];
                origAxis.at(1).at(0)          = this->tetrahedralSymmetry.at(0)[2];
                origAxis.at(2).at(0)          = this->tetrahedralSymmetry.at(0)[3];
                
                origAxis.at(0).at(1)          = this->tetrahedralSymmetry.at(3)[1];
                origAxis.at(1).at(1)          = this->tetrahedralSymmetry.at(3)[2];
                origAxis.at(2).at(1)          = this->tetrahedralSymmetry.at(3)[3];
 
                //============================ Basis vectors for alignment
                hlpVec.emplace_back ( 0.0 );
                basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec );
                
                basisAxis.at(0).at(0)         = 1.0;
                basisAxis.at(1).at(0)         = 0.0;
                basisAxis.at(2).at(0)         = 0.0;
                
                basisAxis.at(0).at(1)         = 0.0;
                basisAxis.at(1).at(1)         = 1.0;
                basisAxis.at(2).at(1)         = 0.0;
                
                basisAxis.at(0).at(2)         = 0.0;
                basisAxis.at(1).at(2)         = 0.0;
                basisAxis.at(2).at(2)         = 1.0;
                
                //============================ Get angle axis representation of the matrix
                std::array<double,5> joinAA   = ProSHADE_internal_misc::getAxisAngleFromSymmetryAxes ( origAxis, basisAxis, 1, 2 );
                
                //============================ Set up rotation ProSHADE run
                ProSHADE_settings* rotSet     = new ProSHADE_settings ( );
                
                // ... Settings regarding resolutions
                rotSet->mapResolution         = setUp->mapResolution;
                
                // ... Settings regarding space around the structure in lattice
                rotSet->extraSpace            = setUp->extraSpace;
                
                // ... Settings regarding rotation module, these should to be set this way!
                rotSet->clearMapData          = false;
                rotSet->useCOM                = false;
                rotSet->rotChangeDefault      = true;
                rotSet->ignoreLs              = std::vector<int> ();
                rotSet->maskBlurFactor        = 500.0;
                rotSet->maskBlurFactorGiven   = false;
                
                // ... Settings regarding the task
                rotSet->taskToPerform         = RotateMap;
                
                // ... Settings regarding where and if to save the clear map
                rotSet->clearMapFile          = setUp->clearMapFile;
                
                // ... Settings regarding the map rotation mode
                rotSet->rotXAxis              =  joinAA[1];
                rotSet->rotYAxis              =  joinAA[2];
                rotSet->rotZAxis              =  joinAA[3];
                rotSet->rotAngle              =  joinAA[0];
                
                // ... Settings regarding the map saving mode
                rotSet->axisOrder             = setUp->axisOrder;
                
                // ... Settings regarding the standard output amount
                rotSet->verbose               = -1;
                
                // ... Settings regarding the file to rotate
                rotSet->structFiles.emplace_back ( setUp->structFiles.at(0) );
                
                //============================ Rotate the structure as required
                if ( setUp->verbose > 0 )
                {
                    std::cout << ">> Structure rotation initiated." << std::endl;
                }
                
                ProSHADE *run                 = new ProSHADE ( rotSet );
                delete run;
                
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Structure with symmetry axis to box axes alignment writted in " << setUp->clearMapFile << std::endl;
                }
            }
            
            //================================ Dihedral symmetry
            else if ( static_cast<unsigned int> ( this->dihedralSymmetries.size() ) > 0 )
            {
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Applying the DIHEDRAL axis alignment convention of Heymann, Chagoyen and Belnap (2005) (principle symm axis on Z, 2-fold on X)." << std::endl;
                }
                
                if ( setUp->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"green\">" << "Applying axes alignment using the convention of Heymann, Chagoyen and Belnap (2005) (principle symm axis on Z, 2-fold on X)." << "</font>";
                    rvapi_set_text            ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                    setUp->htmlReportLineProgress     += 1;
                    rvapi_flush                           ( );
                }
                
                //============================ Take the highest symmetry and consequently highest peak axes
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->dihedralSymmetries.size() ); iter++ )
                {
                    std::sort                 ( this->dihedralSymmetries.at(iter).begin(), this->dihedralSymmetries.at(iter).end(), [](const std::array<double,6>& a, const std::array<double,6>& b) { return a[0] > b[0]; });
                }
                std::sort                     ( this->dihedralSymmetries.begin(), this->dihedralSymmetries.end(), [](const std::vector< std::array<double,6> >& a, const std::vector< std::array<double,6> >& b) { return a.at(0)[0] > b.at(0)[0]; });
                
                double maxFold                = this->dihedralSymmetries.at(0).at(0)[0];
                int noMaxFold                 = 0;
                for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->dihedralSymmetries.size() ); iter++ )
                {
                    if ( this->dihedralSymmetries.at(iter).at(0)[0] == maxFold )
                    {
                        noMaxFold            += 1;
                    }
                }
                std::sort                     ( this->dihedralSymmetries.begin(), this->dihedralSymmetries.begin() + noMaxFold, [](const std::vector< std::array<double,6> >& a, const std::vector< std::array<double,6> >& b) { return a.at(0)[5] > b.at(0)[5]; });
                
                //============================ Original symmetry axes for Z-alignment
                std::vector<std::vector<double>> origAxis;
                std::vector<std::vector<double>> basisAxis;
                std::vector<double> hlpVec;
                hlpVec.emplace_back ( 0.0 ); hlpVec.emplace_back ( 0.0 );
                origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec );
                
                origAxis.at(0).at(0)          = this->dihedralSymmetries.at(0).at(0)[1];
                origAxis.at(1).at(0)          = this->dihedralSymmetries.at(0).at(0)[2];
                origAxis.at(2).at(0)          = this->dihedralSymmetries.at(0).at(0)[3];
                
                origAxis.at(0).at(1)          = this->dihedralSymmetries.at(0).at(1)[1];
                origAxis.at(1).at(1)          = this->dihedralSymmetries.at(0).at(1)[2];
                origAxis.at(2).at(1)          = this->dihedralSymmetries.at(0).at(1)[3];
                
                //============================ Basis vectors for alignment
                hlpVec.emplace_back ( 0.0 );
                basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec );
                
                basisAxis.at(0).at(0)         = 1.0;
                basisAxis.at(1).at(0)         = 0.0;
                basisAxis.at(2).at(0)         = 0.0;
                
                basisAxis.at(0).at(1)         = 0.0;
                basisAxis.at(1).at(1)         = 1.0;
                basisAxis.at(2).at(1)         = 0.0;
                
                basisAxis.at(0).at(2)         = 0.0;
                basisAxis.at(1).at(2)         = 0.0;
                basisAxis.at(2).at(2)         = 1.0;
                
                //============================ Get angle axis representation of the matrix
                std::array<double,5> joinAA   = ProSHADE_internal_misc::getAxisAngleFromSymmetryAxes ( origAxis, basisAxis, 1, 2 );
                
                //============================ Set up rotation ProSHADE run
                ProSHADE_settings* rotSet     = new ProSHADE_settings ( );
                
                // ... Settings regarding resolutions
                rotSet->mapResolution         = setUp->mapResolution;
                
                // ... Settings regarding space around the structure in lattice
                rotSet->extraSpace            = setUp->extraSpace;
                
                // ... Settings regarding rotation module, these should to be set this way!
                rotSet->clearMapData          = false;
                rotSet->useCOM                = false;
                rotSet->rotChangeDefault      = true;
                rotSet->ignoreLs              = std::vector<int> ();
                rotSet->maskBlurFactor        = 500.0;
                rotSet->maskBlurFactorGiven   = false;
                
                // ... Settings regarding the task
                rotSet->taskToPerform         = RotateMap;
                
                // ... Settings regarding where and if to save the clear map
                rotSet->clearMapFile          = setUp->clearMapFile;
                
                // ... Settings regarding the map rotation mode
                rotSet->rotXAxis              =  joinAA[1];
                rotSet->rotYAxis              =  joinAA[2];
                rotSet->rotZAxis              =  joinAA[3];
                rotSet->rotAngle              =  joinAA[0];
                
                // ... Settings regarding the map saving mode
                rotSet->axisOrder             = setUp->axisOrder;
                
                // ... Settings regarding the standard output amount
                rotSet->verbose               = -1;
                
                // ... Settings regarding the file to rotate
                rotSet->structFiles.emplace_back ( setUp->structFiles.at(0) );
                
                //============================ Rotate the structure as required
                if ( setUp->verbose > 0 )
                {
                    std::cout << ">> Structure rotation initiated." << std::endl;
                }
                
                ProSHADE *run                 = new ProSHADE ( rotSet );
                delete run;
                
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Structure with symmetry axis to box axes alignment writted in " << setUp->clearMapFile << std::endl;
                }
            }
            
            //================================ Cyclic symmetry
            else if ( static_cast<unsigned int> ( this->cyclicSymmetries.size() ) > 0 )
            {
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Applying the CYCLIC axis alignment convention of Heymann, Chagoyen and Belnap (2005) (symm axis on Z)." << std::endl;
                }
                
                if ( setUp->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"green\">" << "Applying axis alignment using the convention of Heymann, Chagoyen and Belnap (2005) (symm axis on Z)." << "</font>";
                    rvapi_set_text            ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                    setUp->htmlReportLineProgress += 1;
                    rvapi_flush               ( );
                }
                
                //============================ Take the highest symmetry axis
                std::sort                     ( this->cyclicSymmetries.begin(), this->cyclicSymmetries.end(), [](const std::array<double,5>& a, const std::array<double,5>& b) { return a[0] > b[0]; });
                
                //============================ Original symmetry axes for Z-alignment
                std::vector<std::vector<double>> origAxis;
                std::vector<std::vector<double>> basisAxis;
                std::vector<double> hlpVec;
                hlpVec.emplace_back ( 0.0 );
                origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec ); origAxis.emplace_back( hlpVec );
                
                origAxis.at(0).at(0)          = this->cyclicSymmetries.at(0)[1];
                origAxis.at(1).at(0)          = this->cyclicSymmetries.at(0)[2];
                origAxis.at(2).at(0)          = this->cyclicSymmetries.at(0)[3];
                
                //============================ Basis vectors for alignment
                hlpVec.emplace_back ( 0.0 ); hlpVec.emplace_back ( 0.0 );
                basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec ); basisAxis.emplace_back( hlpVec );
                
                basisAxis.at(0).at(0)         = 1.0;
                basisAxis.at(1).at(0)         = 0.0;
                basisAxis.at(2).at(0)         = 0.0;
                
                basisAxis.at(0).at(1)         = 0.0;
                basisAxis.at(1).at(1)         = 1.0;
                basisAxis.at(2).at(1)         = 0.0;
                
                basisAxis.at(0).at(2)         = 0.0;
                basisAxis.at(1).at(2)         = 0.0;
                basisAxis.at(2).at(2)         = 1.0;
                
                //============================ Get angle axis representation of the matrix
                std::array<double,5> joinAA   = ProSHADE_internal_misc::getAxisAngleFromSymmetryAxes ( origAxis, basisAxis, 1, 2 );
                
                //============================ Set up rotation ProSHADE run
                ProSHADE_settings* rotSet     = new ProSHADE_settings ( );
                
                // ... Settings regarding resolutions
                rotSet->mapResolution         = setUp->mapResolution;
                
                // ... Settings regarding space around the structure in lattice
                rotSet->extraSpace            = setUp->extraSpace;
                
                // ... Settings regarding rotation module, these should to be set this way!
                rotSet->clearMapData          = false;
                rotSet->useCOM                = false;
                rotSet->rotChangeDefault      = true;
                rotSet->ignoreLs              = std::vector<int> ();
                rotSet->maskBlurFactor        = 500.0;
                rotSet->maskBlurFactorGiven   = false;
                
                // ... Settings regarding the task
                rotSet->taskToPerform         = RotateMap;
                
                // ... Settings regarding where and if to save the clear map
                rotSet->clearMapFile          = setUp->clearMapFile;
                
                // ... Settings regarding the map rotation mode
                rotSet->rotXAxis              =  joinAA[1];
                rotSet->rotYAxis              =  joinAA[2];
                rotSet->rotZAxis              =  joinAA[3];
                rotSet->rotAngle              = -joinAA[0];
                
                // ... Settings regarding the map saving mode
                rotSet->axisOrder             = setUp->axisOrder;
                
                // ... Settings regarding the standard output amount
                rotSet->verbose               = -1;
                
                // ... Settings regarding the file to rotate
                rotSet->structFiles.emplace_back ( setUp->structFiles.at(0) );
                
                //============================ Rotate the structure as required
                if ( setUp->verbose > 0 )
                {
                    std::cout << ">> Structure rotation initiated." << std::endl;
                }
                
                ProSHADE *run                 = new ProSHADE ( rotSet );
                delete run;
                
                if ( setUp->verbose > -1 )
                {
                    std::cout << "Structure with symmetry axis to box axes alignment writted in " << setUp->clearMapFile << std::endl;
                }
            }
        }
        
        if ( setUp->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "COMPLETED." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
    }
    
    //======================================== If computing half maps, do so
    if ( setUp->taskToPerform == HalfMaps )
    {
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Half-maps re-boxing</h1>",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1 );
            setUp->htmlReportLine            += 1;
            
            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Starting the re-boxing of half-maps." << "</font>";
            rvapi_set_text                        ( hlpSS.str().c_str(),
                                                   "ProgressSection",
                                                   setUp->htmlReportLineProgress,
                                                   1,
                                                   1,
                                                   1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
        
        ProSHADE_internal::ProSHADE_mapFeatures features ( setUp );
        
    }
    
    //======================================== If simple re-boxing is to be done, do it
    if ( setUp->taskToPerform == SimpleRebox )
    {
        //==================================== Report progress
        if ( setUp->verbose > 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                     MODE: Re-boxing map                 |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Map Re-boxing</h1>",
                                               "body",
                                               setUp->htmlReportLine,
                                               0,
                                               1,
                                               1 );
            setUp->htmlReportLine            += 1;
            
            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            //================================ Sanity check
            if ( setUp->structFiles.size() != 1 )
            {
                std::cerr << "!!! ProSHADE ERROR  !!! Attampted to re-box a map, but supplied ";
                if ( setUp->structFiles.size() < 1 ) { std::cerr << "less "; }
                if ( setUp->structFiles.size() > 1 ) { std::cerr << "more "; }
                std::cerr << "than exactly one file. Please use the \'-f\' option to supply a single map file. Terminating ..." << std::endl << std::endl;
                
                if ( setUp->htmlReport )
                {
                    std::stringstream hlpSS;
                    hlpSS << "<font color=\"red\">" << "Supplied incorrect number (" << setUp->structFiles.size() << ") of structures - a single one expected." << "</font>";
                    rvapi_set_text            ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                settings->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
                    setUp->htmlReportLineProgress += 1;
                    rvapi_flush               ( );
                }
                
                exit ( -1 );
            }
            
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Starting computation of map re-boxing for the input structure " << setUp->structFiles.at(0) << " ." << "</font>";
            rvapi_set_text                        ( hlpSS.str().c_str(),
                                                   "ProgressSection",
                                                   setUp->htmlReportLineProgress,
                                                   1,
                                                   1,
                                                   1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
        
        //==================================== Run
        ProSHADE_internal::ProSHADE_mapFeatures features ( setUp );
    }
    
    //======================================== If computing distances, do so
    if ( setUp->taskToPerform == Distances )
    {
        //==================================== Report progress
        if ( setUp->verbose > 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                       MODE: Distances                   |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Distances</h1>",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1 );
            setUp->htmlReportLine            += 1;

            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Starting computation of distances between input structures." << "</font>";
            rvapi_set_text                        ( hlpSS.str().c_str(),
                                                   "ProgressSection",
                                                   setUp->htmlReportLineProgress,
                                                   1,
                                                   1,
                                                   1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
        
        //==================================== Compute the distances
        ProSHADE_internal::ProSHADE_distances* distObj = new ProSHADE_internal::ProSHADE_distances ( setUp );
        
        //==================================== Save the distances to this object
        if ( setUp->energyLevelDist ) { this->crossCorrDists   = distObj->getEnergyLevelsDistances ( ); }
        if ( setUp->traceSigmaDist )  { this->traceSigmaDists  = distObj->getTraceSigmaDistances ( );   }
        if ( setUp->fullRotFnDist )   { this->rotFunctionDists = distObj->getFullRotationDistances ( ); }
        
        //==================================== Free memory
        delete distObj;
        
        if ( setUp->verbose > 0 )
        {
            std::cout << std::endl << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                         COMPLETED                       |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            //================================ Record progress
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "COMPLETED." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                               "ProgressSection",
                                               setUp->htmlReportLineProgress,
                                               1,
                                               1,
                                               1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
        
        if ( setUp->verbose > 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                          RESULTS                        |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            //================================ Record progress
            rvapi_add_section                 ( "ResultsSection",
                                                "Distances",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            rvapi_add_grid                    ( "distGrid",
                                                true,
                                                "ResultsSection",
                                                2,
                                                0,
                                                1,
                                                1 );
            
            rvapi_add_table                   ( "DistancesTable",
                                                "Distances between structures",
                                                "distGrid",
                                                0,
                                                0,
                                                1,
                                                1,
                                                1 );
            
            // ... Column headers
            std::stringstream hlpSS;
            int columnIter                    = 0;
            
            if ( setUp->energyLevelDist )
            {
                hlpSS.str                     ( std::string ( ) );
                hlpSS << "This column contains the distances obtained by computing correlations between the spherical harmonics on different shells and then comparing these correlation tables between the two structures. These distances are all from structure " << setUp->structFiles.at(0) << " and the structure named in the row name.";
                rvapi_put_horz_theader        ( "DistancesTable", "Energy Levels Distances", hlpSS.str().c_str(), columnIter );
                columnIter                   += 1;
            }
            
            if ( setUp->energyLevelDist )
            {
                hlpSS.str                     ( std::string ( ) );
                hlpSS << "This column contains the distances obtained by calculating the distances in spherical harmonics coefficients integrated over all the shells, using Singular Value Decomposition to approximate the rotation difference. These distances are all from structure " << setUp->structFiles.at(0) << " and the structure named in the row name.";
                rvapi_put_horz_theader        ( "DistancesTable", "Trace Sigma Distances", hlpSS.str().c_str(), columnIter );
                columnIter                   += 1;
            }
            
            if ( setUp->energyLevelDist )
            {
                hlpSS.str                     ( std::string ( ) );
                hlpSS << "This column contains the distances obtained by calculating the distances in spherical harmonics coefficients integrated over all the shells, using the Rotation Function to find the optimal overlay rotation between the structures. These distances are all from structure " << setUp->structFiles.at(0) << " and the structure named in the row name.";
                rvapi_put_horz_theader        ( "DistancesTable", "Rotation Function Distances", hlpSS.str().c_str(), columnIter );
                columnIter                   += 1;
            }
            
            // ... Row headers
            for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( setUp->structFiles.size() ); iter++ )
            {
                hlpSS.str                     ( std::string ( ) );
                hlpSS << setUp->structFiles.at(iter);
                rvapi_put_vert_theader        ( "DistancesTable", hlpSS.str().c_str(), "", iter-1 );
            }
            
            rvapi_flush                       ( );
        }

        if ( setUp->energyLevelDist )
        {
            if ( setUp->verbose != -1 )
            {
                printf ( "Energy Level Descriptor distances      :           %+.4f", this->crossCorrDists.at(0) );
                for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( this->crossCorrDists.size() ); iter++ )
                {
                    printf ( " %+.4f", this->crossCorrDists.at(iter) );
                }
                std::cout << std::endl;
            }
            
        }

        if ( setUp->traceSigmaDist )
        {
            if ( setUp->verbose != -1 )
            {
                printf ( "Trace Sigma Descriptor distances       :           %+.4f", this->traceSigmaDists.at(0) );
                for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( this->traceSigmaDists.size() ); iter++ )
                {
                    printf ( " %+.4f", this->traceSigmaDists.at(iter) );
                }
                std::cout << std::endl;
            }
        }
        
        if ( setUp->fullRotFnDist )
        {
            if ( setUp->verbose != -1 )
            {
                printf ( "Rotation Function Descriptor distances :           %+.4f", this->rotFunctionDists.at(0) );
                for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( this->rotFunctionDists.size() ); iter++ )
                {
                    printf ( " %+.4f", this->rotFunctionDists.at(iter) );
                }
                std::cout << std::endl;
            }
        }
        
        if ( setUp->verbose != -1 )
        {
            for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( setUp->structFiles.size() ); iter++ )
            {
                std::cout << "Matching structure names               :   " << setUp->structFiles.at(iter) << std::endl;
            }
            std::cout << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            int colIter                       = 0;
            std::stringstream hlpSS;
            
            if ( setUp->energyLevelDist )
            {
                for ( unsigned int rIter = 0; rIter < static_cast<unsigned int> ( this->crossCorrDists.size() ); rIter++ )
                {
                    hlpSS.str                 ( std::string ( ) );
                    hlpSS << this->crossCorrDists.at(rIter);
                    rvapi_put_table_string    ( "DistancesTable", hlpSS.str().c_str(), rIter, colIter );
                }
                colIter                      += 1;
            }
            
            if ( setUp->traceSigmaDist )
            {
                for ( unsigned int rIter = 0; rIter < static_cast<unsigned int> ( this->traceSigmaDists.size() ); rIter++ )
                {
                    hlpSS.str                 ( std::string ( ) );
                    hlpSS << this->traceSigmaDists.at(rIter);
                    rvapi_put_table_string    ( "DistancesTable", hlpSS.str().c_str(), rIter, colIter );
                }
                colIter                      += 1;
            }
            
            if ( setUp->fullRotFnDist )
            {
                for ( unsigned int rIter = 0; rIter < static_cast<unsigned int> ( this->rotFunctionDists.size() ); rIter++ )
                {
                    hlpSS.str                 ( std::string ( ) );
                    hlpSS << this->rotFunctionDists.at(rIter);
                    rvapi_put_table_string    ( "DistancesTable", hlpSS.str().c_str(), rIter, colIter );
                }
                colIter                      += 1;
            }
            rvapi_flush                       ( );
        }
        
        if ( ( setUp->verbose > 0 ) && ( setUp->databaseName != "" ) )
        {
            std::cout << "NOTE: If you do not see a structure you were expecting and know to be in the database, it may be caused by the querry structure and the one you were expecting having too different dimmensions - see the help dialogue for parameter \'--dbSizeLim\'." << std::endl << std::endl;
        }
        
        this->distancesAvailable              = true;

    }
    
    //======================================== If computing fragments distances, do so
    if ( setUp->taskToPerform == DistancesFrag )
    {
        //==================================== Report progress
        if ( setUp->verbose != 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                     MODE: DistancesFrag                 |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            //================================ Report title
            rvapi_set_text                    ( "<h1>ProSHADE Results: Searching fragments in database</h1>",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1 );
            setUp->htmlReportLine            += 1;
            
            //==================================== Create section
            rvapi_add_section                 ( "ProgressSection",
                                                "Progress",
                                                "body",
                                                setUp->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            setUp->htmlReportLine            += 1;
            
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "Starting fragment search of file " << setUp->structFiles.at(0) << " against the database." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            rvapi_flush                       ( );
        }
        
        ProSHADE_internal::ProSHADE_distances* distObj = new ProSHADE_internal::ProSHADE_distances ( setUp );
        
        if ( setUp->verbose > 0 )
        {
            std::cout << std::endl << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                         COMPLETED                       |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        if ( setUp->htmlReport )
        {
            std::stringstream hlpSS;
            hlpSS << "<font color=\"green\">" << "COMPLETED." << "</font>";
            rvapi_set_text                    ( hlpSS.str().c_str(),
                                                "ProgressSection",
                                                setUp->htmlReportLineProgress,
                                                1,
                                                1,
                                                1 );
            setUp->htmlReportLineProgress    += 1;
            
            rvapi_flush                       ( );
        }
        
        if ( setUp->verbose > 0 )
        {
            std::cout << "-----------------------------------------------------------" << std::endl;
            std::cout << "|                          RESULTS                        |" << std::endl;
            std::cout << "-----------------------------------------------------------" << std::endl << std::endl;
        }
        
        std::vector< std::vector<double> > enLevDists;
        std::vector< std::vector<double> > trSigmaDists;
        std::vector< std::vector<double> > fullRotFnDists;
        if ( setUp->energyLevelDist )
        {
            enLevDists                        = distObj->getFragEnergyLevelsDistances ( );
        }
        
        if ( setUp->traceSigmaDist )
        {
            trSigmaDists                      = distObj->getFragTraceSigmaDistances ( );
        }
        
        if ( setUp->fullRotFnDist )
        {
            fullRotFnDists                    = distObj->getFragFullRotationDistances ( );
        }
        
        if ( setUp->verbose > -1 )
        {
            double maxFrag                    = std::max ( enLevDists.size(), std::max ( trSigmaDists.size(), fullRotFnDists.size() ) );
            for ( unsigned int frag = 0; frag < static_cast<unsigned int> ( maxFrag ); frag++ )
            {
                printf ( "Fragment %d Results:\n", frag );
                printf ( "--------------------\n" );
                if ( setUp->energyLevelDist )
                {
                    printf ( "Energy levels descriptor distances       :   %+.4f", enLevDists.at(frag).at(0) );
                    for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( enLevDists.at(frag).size() ); iter++ )
                    {
                        printf ( " %+.4f", enLevDists.at(frag).at(iter) );
                    }
                    std::cout << std::endl;
                }
                if ( setUp->traceSigmaDist )
                {
                    printf ( "Trace sigma descriptor distances         :   %+.4f", trSigmaDists.at(frag).at(0) );
                    for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( trSigmaDists.at(frag).size() ); iter++ )
                    {
                        printf ( " %+.4f", trSigmaDists.at(frag).at(iter) );
                    }
                    std::cout << std::endl;
                }
                if ( setUp->fullRotFnDist )
                {
                    printf ( "Full RF descriptor distances             :   %+.4f", fullRotFnDists.at(frag).at(0) );
                    for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( fullRotFnDists.at(frag).size() ); iter++ )
                    {
                        printf ( " %+.4f", fullRotFnDists.at(frag).at(iter) );
                    }
                    std::cout << std::endl;
                }
                std::cout << std::endl;
            }
            
            std::cout << std::endl;
            for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( setUp->structFiles.size() ); iter++ )
            {
                std::cout << "Matching structure names               :   " << setUp->structFiles.at(iter) << std::endl;
            }
        }
        
        if ( setUp->htmlReport )
        {
            int maxFrag                       = std::max ( enLevDists.size(), std::max ( trSigmaDists.size(), fullRotFnDists.size() ) );

            //================================ Create section
            rvapi_add_section                 ( "ResultsSection",
                                                "Fragment Distances",
                                                "body",
                                                settings->htmlReportLine,
                                                0,
                                                1,
                                                1,
                                                true );
            settings->htmlReportLine         += 1;
            
            //================================ Fragment distances table for Enetry Levels
            rvapi_add_table                   ( "FragmentDistancesELTable",
                                                "Fragment Energy Level distances to database entries",
                                                "ResultsSection",
                                                0,
                                                0,
                                                enLevDists.size(),
                                                enLevDists.at(0).size(),
                                                1 );
            
            // ... Column headers
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( enLevDists.at(0).size() ); iter++ )
            {
                std::stringstream hlpSS;
                hlpSS << "Database entry " << setUp->structFiles.at(iter);
                std::stringstream hlpSS2;
                hlpSS2 << "This is the Energy Levels distance for the fragment in the row to the database entry listed here (" << setUp->structFiles.at(iter) << ")";
                rvapi_put_horz_theader            ( "FragmentDistancesELTable", hlpSS.str().c_str(), hlpSS2.str().c_str(), iter );
            }
            
            // ... Row headers
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( maxFrag ); iter++ )
            {
                std::stringstream hlpSS;
                hlpSS << "Fragment Number " << iter;
                std::stringstream hlpSS2;
                hlpSS2 << "This is the Energy Levels distance for the fragment number " << iter << " to the database entry in the column header.";
                rvapi_put_vert_theader        ( "FragmentDistancesELTable", hlpSS.str().c_str(), hlpSS2.str().c_str(), iter );
            }
            
            // ... Fill in data
            for ( unsigned int rIter = 0; rIter < static_cast<unsigned int> ( enLevDists.size() ); rIter++ )
            {
                for ( unsigned int cIter = 0; cIter < static_cast<unsigned int> ( enLevDists.at(rIter).size() ); cIter++ )
                {
                    std::stringstream hlpSS;
                    hlpSS << std::showpos << ProSHADE_internal_misc::roundDouble ( enLevDists.at(rIter).at(cIter) * 1000.0 ) / 1000.0;
                    if ( hlpSS.str().length() != 6 ) { int hlp = 6 - hlpSS.str().length(); for ( int i = 0; i < hlp; i++ ) { hlpSS << " "; } }
                    if ( hlpSS.str().length() > 6 )  { hlpSS.str( hlpSS.str().substr( 0, 6 ) ); }
                    rvapi_put_table_string    ( "FragmentDistancesELTable", hlpSS.str().c_str(), rIter, cIter );
                }
            }
            
            //================================ Fragment distances table for Trace Sigma
            rvapi_add_table                   ( "FragmentDistancesTSTable",
                                                "Fragment Trace Sigma distances to database entries",
                                                "ResultsSection",
                                                enLevDists.size() + 1,
                                                0,
                                                trSigmaDists.size(),
                                                trSigmaDists.at(0).size(),
                                                1 );
            
            // ... Column headers
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( enLevDists.at(0).size() ); iter++ )
            {
                std::stringstream hlpSS;
                hlpSS << "Database entry " << setUp->structFiles.at(iter);
                std::stringstream hlpSS2;
                hlpSS2 << "This is the Trace Sigma distance for the fragment in the row to the database entry listed here (" << setUp->structFiles.at(iter) << ")";
                rvapi_put_horz_theader            ( "FragmentDistancesTSTable", hlpSS.str().c_str(), hlpSS2.str().c_str(), iter );
            }
            
            // ... Row headers
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( maxFrag ); iter++ )
            {
                std::stringstream hlpSS;
                hlpSS << "Fragment Number " << iter;
                std::stringstream hlpSS2;
                hlpSS2 << "This is the Trace Sigma distance for the fragment number " << iter << " to the database entry in the column header.";
                rvapi_put_vert_theader        ( "FragmentDistancesTSTable", hlpSS.str().c_str(), hlpSS2.str().c_str(), iter );
            }
            
            // ... Fill in data
            for ( unsigned int rIter = 0; rIter < static_cast<unsigned int> ( trSigmaDists.size() ); rIter++ )
            {
                for ( unsigned int cIter = 0; cIter < static_cast<unsigned int> ( trSigmaDists.at(rIter).size() ); cIter++ )
                {
                    std::stringstream hlpSS;
                    hlpSS << std::showpos << ProSHADE_internal_misc::roundDouble ( trSigmaDists.at(rIter).at(cIter) * 1000.0 ) / 1000.0;
                    if ( hlpSS.str().length() != 6 ) { int hlp = 6 - hlpSS.str().length(); for ( int i = 0; i < hlp; i++ ) { hlpSS << " "; } }
                    if ( hlpSS.str().length() > 6 )  { hlpSS.str( hlpSS.str().substr( 0, 6 ) ); }
                    rvapi_put_table_string    ( "FragmentDistancesTSTable", hlpSS.str().c_str(), rIter, cIter );
                }
            }
            
            //================================ Fragment distances table for Rotation Function
            rvapi_add_table                   ( "FragmentDistancesRFTable",
                                                "Fragment Rotation Function distances to database entries",
                                                "ResultsSection",
                                                enLevDists.size() + trSigmaDists.size() + 1,
                                                0,
                                                fullRotFnDists.size(),
                                                fullRotFnDists.at(0).size(),
                                                1 );
            
            // ... Column headers
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( fullRotFnDists.at(0).size() ); iter++ )
            {
                std::stringstream hlpSS;
                hlpSS << "Database entry " << setUp->structFiles.at(iter);
                std::stringstream hlpSS2;
                hlpSS2 << "This is the Rotation Function distance for the fragment in the row to the database entry listed here (" << setUp->structFiles.at(iter) << ")";
                rvapi_put_horz_theader            ( "FragmentDistancesRFTable", hlpSS.str().c_str(), hlpSS2.str().c_str(), iter );
            }
            
            // ... Row headers
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( maxFrag ); iter++ )
            {
                std::stringstream hlpSS;
                hlpSS << "Fragment Number " << iter;
                std::stringstream hlpSS2;
                hlpSS2 << "This is the Rotation Function distance for the fragment number " << iter << " to the database entry in the column header.";
                rvapi_put_vert_theader        ( "FragmentDistancesRFTable", hlpSS.str().c_str(), hlpSS2.str().c_str(), iter );
            }
            
            // ... Fill in data
            for ( unsigned int rIter = 0; rIter < static_cast<unsigned int> ( fullRotFnDists.size() ); rIter++ )
            {
                for ( unsigned int cIter = 0; cIter < static_cast<unsigned int> ( fullRotFnDists.at(rIter).size() ); cIter++ )
                {
                    std::stringstream hlpSS;
                    hlpSS << std::showpos << ProSHADE_internal_misc::roundDouble ( fullRotFnDists.at(rIter).at(cIter) * 1000.0 ) / 1000.0;
                    if ( hlpSS.str().length() != 6 ) { int hlp = 6 - hlpSS.str().length(); for ( int i = 0; i < hlp; i++ ) { hlpSS << " "; } }
                    if ( hlpSS.str().length() > 6 )  { hlpSS.str( hlpSS.str().substr( 0, 6 ) ); }
                    rvapi_put_table_string    ( "FragmentDistancesRFTable", hlpSS.str().c_str(), rIter, cIter );
                }
            }
            
            rvapi_flush                       ( );
        }
        
        delete distObj;
    }
    
    //======================================== Done
    
}

ProSHADE::ProSHADE::~ProSHADE ( )
{
    //======================================== Clear vectors
    this->crossCorrDists.clear                ( );
    this->traceSigmaDists.clear               ( );
    this->rotFunctionDists.clear              ( );
    this->cyclicSymmetries.clear              ( );
    this->cyclicSymmetriesClear.clear         ( );
    this->dihedralSymmetries.clear            ( );
    this->tetrahedralSymmetry.clear           ( );
    this->octahedralSymmetry.clear            ( );
    this->icosahedralSymmetry.clear           ( );
    this->fragmentList.clear                  ( );
    
    //======================================== Done
    
}

std::string ProSHADE::ProSHADE::getProSHADEVersion ( )
{
    //======================================== Check for sanity
    return ( std::string ( __PROSHADE_VERSION__ ) );
    
    //======================================== Done
    
}

void ProSHADE::ProSHADE_settings::appendStructure ( std::string str )
{
    //======================================== Check for sanity
    this->structFiles.emplace_back ( str );
    
    //======================================== Done
    
}

std::vector<std::string> ProSHADE::ProSHADE::getMapFragments ( )
{
    //======================================== Check for sanity
    if ( !this->fragmentsAvailable )
    {
        std::cerr << "!!! ProSHADE ERROR !!! Attempted to obtain list of map fragments, but did not compute the fragmentation beforehand. This is an undefined behaviour - terminating..." << std::endl;
        exit                                  ( -1 );
    }
    
    //======================================== Done
    return                                    ( this->fragmentList );
}

std::vector<double> ProSHADE::ProSHADE::getCrossCorrDists ( )
{
    //======================================== Check for sanity
    if ( this->distancesAvailable )
    {
        return ( this->crossCorrDists );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the cross-correlation distances without requesting their computation first. Returning empty vector." << std::endl;
        return ( std::vector<double> () );
    }
        
    //======================================== Done
    
}

std::vector<double> ProSHADE::ProSHADE::getTraceSigmaDists ( )
{
    //======================================== Check for sanity
    if ( this->distancesAvailable )
    {
        return ( this->traceSigmaDists );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the trace sigma distances without requesting their computation first. Returning empty vector." << std::endl;
        return ( std::vector<double> () );
    }
    
    //======================================== Done
    
}

std::vector<double> ProSHADE::ProSHADE::getRotFunctionDists ( )
{
    //======================================== Check for sanity
    if ( this->distancesAvailable )
    {
        return ( this->rotFunctionDists );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the rotation function based distances without requesting their computation first. Returning empty vector." << std::endl;
        return ( std::vector<double> () );
    }
    
    //======================================== Done
    
}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getCyclicSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->cyclicSymmetries );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the cyclic symmetries list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    
}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getClearCyclicSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->cyclicSymmetriesClear );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the clear cyclic symmetries list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    
}

std::vector< double > ProSHADE::ProSHADE::getCyclicSymmetriesPy ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< double > ret;
        
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->cyclicSymmetries.size() ); iter++ )
        {
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( this->cyclicSymmetries.at(iter).size() ); it++ )
            {
                ret.emplace_back              ( this->cyclicSymmetries.at(iter).at(it) );
            }
            ret.emplace_back                  ( -999.999 );
        }
        
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the cyclic symmetries list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< double > () );
    }
    
    //======================================== Done
    
}

std::vector< double > ProSHADE::ProSHADE::getDihedralSymmetriesPy ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< double > ret;
        
        for ( unsigned int i = 0; i < static_cast<unsigned int> ( this->dihedralSymmetries.size() ); i++ )
        {
            for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->dihedralSymmetries.at(i).size() ); iter++ )
            {
                for ( unsigned int it = 0; it < static_cast<unsigned int> ( this->dihedralSymmetries.at(i).at(iter).size() ); it++ )
                {
                    ret.emplace_back              ( this->dihedralSymmetries.at(i).at(iter).at(it) );
                }
                ret.emplace_back                  ( -999.999 );
            }
            ret.emplace_back                  ( -777.777 );
        }
        
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the dihedral symmetries list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< double > () );
    }
    
    //======================================== Done
    
}

std::vector< std::vector<std::array<double,6> > > ProSHADE::ProSHADE::getDihedralSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->dihedralSymmetries );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the dihedral symmetries list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::vector< std::array<double,6> > > () );
    }
    
    //======================================== Done
    
}

std::vector< std::vector<std::array<double,6> > > ProSHADE::ProSHADE::getClearDihedralSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->dihedralSymmetriesClear );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the clear dihedral symmetries list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::vector< std::array<double,6> > > () );
    }
    
    //======================================== Done
    
}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getTetrahedralSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->tetrahedralSymmetry );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the tetrahedral symmetry list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    
}

std::vector< double > ProSHADE::ProSHADE::getTetrahedralSymmetriesPy ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< double > ret;
        
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->tetrahedralSymmetry.size() ); iter++ )
        {
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( this->tetrahedralSymmetry.at(iter).size() ); it++ )
            {
                ret.emplace_back              ( this->tetrahedralSymmetry.at(iter).at(it) );
            }
            ret.emplace_back                  ( -999.999 );
        }
        
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the tetrahedral symmetry list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< double > () );
    }
    
    //======================================== Done
    
}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getOctahedralSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->octahedralSymmetry );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the octahedral symmetry list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    
}

std::vector< double > ProSHADE::ProSHADE::getOctahedralSymmetriesPy ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< double > ret;
        
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->octahedralSymmetry.size() ); iter++ )
        {
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( this->octahedralSymmetry.at(iter).size() ); it++ )
            {
                ret.emplace_back              ( this->octahedralSymmetry.at(iter).at(it) );
            }
            ret.emplace_back                  ( -999.999 );
        }
        
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the octahedral symmetry list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< double > () );
    }
    
    //======================================== Done
    
}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getIcosahedralSymmetries ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        return ( this->icosahedralSymmetry );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the icosahedral symmetry list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    
}

std::vector< double > ProSHADE::ProSHADE::getIcosahedralSymmetriesPy ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< double > ret;
        
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->icosahedralSymmetry.size() ); iter++ )
        {
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( this->icosahedralSymmetry.at(iter).size() ); it++ )
            {
                ret.emplace_back              ( this->icosahedralSymmetry.at(iter).at(it) );
            }
            ret.emplace_back                  ( -999.999 );
        }
        
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the icosahedral symmetry list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< double > () );
    }
    
    //======================================== Done
    
}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getRecommendedSymmetry ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        //==================================== Initialise
        ProSHADE_internal::ProSHADE_symmetry symObj;
        std::vector< std::array<double,5> > ret;
        std::array<double,5> hlpArr;
        
        //==================================== Identity element
        hlpArr[0]                             = 1.0;
        hlpArr[1]                             = 1.0;
        hlpArr[2]                             = 0.0;
        hlpArr[3]                             = 0.0;
        hlpArr[4]                             = 0.0;
        ret.emplace_back                      ( hlpArr );
        
        //==================================== If icosahedral exist, return that
        if ( this->icosahedralSymmetry.size() < 1 )
        {
            this->icosahedralSymmetry         = getIcosahedralSymmetries ( );
        }
        
        if ( static_cast<unsigned int> ( this->icosahedralSymmetry.size() ) > 0 )
        {
            std::vector< std::array<double,5> > hlp = symObj.generateIcosElements ( this->icosahedralSymmetry, this->settings, -1 );
            for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
            {
                ret.emplace_back              ( hlp.at(iter) );
            }
            return ( ret );
        }
        
        //==================================== If octahedral exist, return that
        if ( this->octahedralSymmetry.size() < 1 )
        {
            this->octahedralSymmetry          = getOctahedralSymmetries ( );
        }
        
        if ( static_cast<unsigned int> ( this->octahedralSymmetry.size() ) > 0 )
        {
            std::vector< std::array<double,5> > hlp = symObj.generateOctaElements ( this->octahedralSymmetry, this->settings, -1 );
            for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
            {
                ret.emplace_back              ( hlp.at(iter) );
            }
            return ( ret );
        }
        
        //==================================== If tetrahedral exist, return that
        if ( this->tetrahedralSymmetry.size() < 1 )
        {
            this->tetrahedralSymmetry         = getTetrahedralSymmetries ( );
        }
        
        if ( static_cast<unsigned int> ( this->tetrahedralSymmetry.size() ) > 0 )
        {
            std::vector< std::array<double,5> > hlp = symObj.generateTetrElements ( this->tetrahedralSymmetry, this->settings, -1 );
            for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
            {
                ret.emplace_back              ( hlp.at(iter) );
            }
            return ( ret );
        }
        
        //==================================== If dihedral exist, return that
        if ( this->dihedralSymmetries.size() < 1 )
        {
            this->dihedralSymmetries          = getDihedralSymmetries ( );
        }
        this->dihedralSymmetriesClear         = getClearDihedralSymmetries ( );
        
        if ( static_cast<unsigned int> ( this->dihedralSymmetriesClear.size() ) > 0 )
        {
            for ( unsigned int it = 0; it < 2; it++ )
            {
                if ( static_cast<int> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) % 2 == 0 )
                {
                    for ( int iter = -std::ceil( static_cast<int> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) / 2.0 ); iter <= std::ceil( static_cast<int> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) / 2.0 ); iter++ )
                    {
                        if ( iter == 0 ) { continue; }
                        if ( iter == -std::ceil( static_cast<int> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) / 2.0 ) ) { continue; }
                        hlpArr[0]             = this->dihedralSymmetriesClear.at(0).at(it)[0];
                        hlpArr[1]             = this->dihedralSymmetriesClear.at(0).at(it)[1];
                        hlpArr[2]             = this->dihedralSymmetriesClear.at(0).at(it)[2];
                        hlpArr[3]             = this->dihedralSymmetriesClear.at(0).at(it)[3];
                        hlpArr[4]             = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) );
                        ret.emplace_back      ( hlpArr );
                    }
                }
                else
                {
                    for ( int iter = -std::floor( static_cast<int> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) / 2.0 ); iter <= std::floor( static_cast<int> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) / 2.0 ); iter++ )
                    {
                        if ( iter == 0 ) { continue; }
                        hlpArr[0]             = this->dihedralSymmetriesClear.at(0).at(it)[0];
                        hlpArr[1]             = this->dihedralSymmetriesClear.at(0).at(it)[1];
                        hlpArr[2]             = this->dihedralSymmetriesClear.at(0).at(it)[2];
                        hlpArr[3]             = this->dihedralSymmetriesClear.at(0).at(it)[3];
                        hlpArr[4]             = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->dihedralSymmetriesClear.at(0).at(it)[0] ) );
                        ret.emplace_back      ( hlpArr );
                    }
                }
            }
            
            return ( ret );
        }
        
        //==================================== If cyclic exist, return that
        if ( this->cyclicSymmetries.size() < 1 )
        {
            this->cyclicSymmetries            = getCyclicSymmetries ( );
        }
        this->cyclicSymmetriesClear           = getClearCyclicSymmetries ( );
        
        if ( static_cast<unsigned int> ( this->cyclicSymmetriesClear.size() ) > 0 )
        {
            if ( static_cast<int> ( this->cyclicSymmetriesClear.at(0)[0] ) % 2 == 0 )
            {
                for ( int iter = -std::ceil( static_cast<int> ( this->cyclicSymmetriesClear.at(0)[0] ) / 2.0 ); iter <= std::ceil( static_cast<int> ( this->cyclicSymmetriesClear.at(0)[0] ) / 2.0 ); iter++ )
                {
                    if ( iter == 0 )  { continue; }
                    if ( iter == -std::ceil( static_cast<int> ( this->cyclicSymmetriesClear.at(0)[0] ) / 2.0 ) ) { continue; }
                    hlpArr[0]         = this->cyclicSymmetriesClear.at(0)[0];
                    hlpArr[1]         = this->cyclicSymmetriesClear.at(0)[1];
                    hlpArr[2]         = this->cyclicSymmetriesClear.at(0)[2];
                    hlpArr[3]         = this->cyclicSymmetriesClear.at(0)[3];
                    hlpArr[4]         = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->cyclicSymmetriesClear.at(0)[0] ) );
                    ret.emplace_back                      ( hlpArr );
                }
            }
            else
            {
                for ( int iter = -std::floor( static_cast<int> ( this->cyclicSymmetriesClear.at(0)[0] ) / 2.0 ); iter <= std::floor( static_cast<int> ( this->cyclicSymmetriesClear.at(0)[0] ) / 2.0 ); iter++ )
                {
                    if ( iter == 0 )  { continue; }
                    hlpArr[0]         = this->cyclicSymmetriesClear.at(0)[0];
                    hlpArr[1]         = this->cyclicSymmetriesClear.at(0)[1];
                    hlpArr[2]         = this->cyclicSymmetriesClear.at(0)[2];
                    hlpArr[3]         = this->cyclicSymmetriesClear.at(0)[3];
                    hlpArr[4]         = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->cyclicSymmetriesClear.at(0)[0] ) );
                    ret.emplace_back                      ( hlpArr );
                }
            }
            
            return ( ret );
        }
        
        std::cout << "!!! ProSHADE WARNING !!! Did not find any symmetries to report. Returning empty vector..." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the symmetry elements list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    std::cerr << "!!! ProSHADE ERROR !!! Something went wery wrong... Please report this. Returning empty vector..." << std::endl;
    return ( std::vector< std::array<double,5> > () );
}

std::vector< double > ProSHADE::ProSHADE::getRecommendedSymmetryElementsPy ( )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< std::array< double, 5 > > hlp = this->getRecommendedSymmetry ( );
        std::vector< double > ret;

        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
        {
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( hlp.at(iter).size() ); it++ )
            {
                ret.emplace_back              ( hlp.at(iter).at(it) );
            }
            ret.emplace_back                  ( -999.999 );
        }

        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the symmetry elements list without requesting their computation first. Returning empty vector." << std::endl;
        return ( std::vector< double > () );
    }

    //======================================== Done

}

std::vector< std::array<double,5> > ProSHADE::ProSHADE::getSpecificSymmetryElements ( std::string symType,
                                                                                      int symFold )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        //==================================== Initialise
        ProSHADE_internal::ProSHADE_symmetry symObj;
        std::vector< std::array<double,5> > ret;
        std::array<double,5> hlpArr;
        
        //==================================== Identity element
        hlpArr[0]                             = 1.0;
        hlpArr[1]                             = 1.0;
        hlpArr[2]                             = 0.0;
        hlpArr[3]                             = 0.0;
        hlpArr[4]                             = 0.0;
        ret.emplace_back                      ( hlpArr );
        
        //==================================== If no fold is given, use the highest scoring one
        if ( symFold == 0 )
        {
            if ( symType == "D" )
            {
                //============================ If not yet loaded, get them
                if ( this->dihedralSymmetries.size() < 1 )
                {
                    this->dihedralSymmetries  = getDihedralSymmetries ( );
                }
                
                //============================ Take highest peak
                symFold                       = this->dihedralSymmetries.at(0).at(0)[0];
            }
            else
            {
                //============================ If not yet loaded, get them
                if ( this->cyclicSymmetries.size() < 1 )
                {
                    this->cyclicSymmetries    = getCyclicSymmetries ( );
                }
                
                //============================ Take highest peak
                symFold                       = this->cyclicSymmetries.at(0)[0];
            }
        }
        
        //==================================== Cyclic symmetry
        if ( symType == "C" )
        {
            //================================ If not yet loaded, get them
            if ( this->cyclicSymmetries.size() < 1 )
            {
                this->cyclicSymmetries        = getCyclicSymmetries ( );
            }
            
            //================================ Locate the highest cyclic symmetry with required fold
            bool foundSym                     = false;
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( this->cyclicSymmetries.size() ); it++ )
            {
                if ( this->cyclicSymmetries.at(it)[0] == symFold )
                {
                    //======================== Init
                    foundSym                  = true;
                    
                    //======================== Find elements
                    if ( static_cast<int> ( this->cyclicSymmetries.at(it)[0] ) % 2 == 0 )
                    {
                        for ( int iter = -std::ceil( static_cast<int> ( this->cyclicSymmetries.at(it)[0] ) / 2.0 ); iter <= std::ceil( static_cast<int> ( this->cyclicSymmetries.at(it)[0] ) / 2.0 ); iter++ )
                        {
                            if ( iter == 0 )  { continue; }
                            if ( iter == -std::ceil( static_cast<int> ( this->cyclicSymmetries.at(it)[0] ) / 2.0 ) ) { continue; }
                            hlpArr[0]         = this->cyclicSymmetries.at(it)[0];
                            hlpArr[1]         = this->cyclicSymmetries.at(it)[1];
                            hlpArr[2]         = this->cyclicSymmetries.at(it)[2];
                            hlpArr[3]         = this->cyclicSymmetries.at(it)[3];
                            hlpArr[4]         = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->cyclicSymmetries.at(it)[0] ) );
                            ret.emplace_back                      ( hlpArr );
                        }
                    }
                    else
                    {
                        for ( int iter = -std::floor( static_cast<int> ( this->cyclicSymmetries.at(it)[0] ) / 2.0 ); iter <= std::floor( static_cast<int> ( this->cyclicSymmetries.at(it)[0] ) / 2.0 ); iter++ )
                        {
                            if ( iter == 0 )  { continue; }
                            hlpArr[0]         = this->cyclicSymmetries.at(it)[0];
                            hlpArr[1]         = this->cyclicSymmetries.at(it)[1];
                            hlpArr[2]         = this->cyclicSymmetries.at(it)[2];
                            hlpArr[3]         = this->cyclicSymmetries.at(it)[3];
                            hlpArr[4]         = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->cyclicSymmetries.at(it)[0] ) );
                            ret.emplace_back                      ( hlpArr );
                        }
                    }
                    
                    //======================== Done
                    break;
                }
            }
            
            //================================ Report if required fold was not found
            if ( !foundSym )
            {
                std::cout << "!!! ProSHADE WARNING !!! Requested symmetry elements for " << symType << "-" << symFold << " symmetry, but this symmetry was not found. Returning empty vector." << std::endl;
                return ( std::vector< std::array<double,5> > () );
            }
        }
        
        //==================================== Dihedral symmetry
        if ( symType == "D" )
        {
            //================================ If not yet loaded, get them
            if ( this->dihedralSymmetries.size() < 1 )
            {
                this->dihedralSymmetries      = getDihedralSymmetries ( );
            }
            
            //================================ Locate the highest cyclic symmetry with required fold
            bool foundSym                     = false;
            for ( unsigned int i = 0; i < static_cast<unsigned int> ( this->dihedralSymmetries.size() ); i++ )
            {
                if ( this->dihedralSymmetries.at(i).at(0)[0] == symFold )
                {
                    //======================== Init
                    foundSym                  = true;
                    
                    //======================== Find elements
                    for ( unsigned int it = 0; it < 2; it++ )
                    {
                        if ( static_cast<int> ( this->dihedralSymmetries.at(i).at(it)[0] ) % 2 == 0 )
                        {
                            for ( int iter = -std::ceil( static_cast<int> ( this->dihedralSymmetries.at(i).at(it)[0] ) / 2.0 ); iter <= std::ceil( static_cast<int> ( this->dihedralSymmetries.at(i).at(it)[0] ) / 2.0 ); iter++ )
                            {
                                if ( iter == 0 ) { continue; }
                                if ( iter == -std::ceil( static_cast<int> ( this->dihedralSymmetries.at(i).at(it)[0] ) / 2.0 ) ) { continue; }
                                hlpArr[0]     = this->dihedralSymmetries.at(i).at(it)[0];
                                hlpArr[1]     = this->dihedralSymmetries.at(i).at(it)[1];
                                hlpArr[2]     = this->dihedralSymmetries.at(i).at(it)[2];
                                hlpArr[3]     = this->dihedralSymmetries.at(i).at(it)[3];
                                hlpArr[4]     = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->dihedralSymmetries.at(i).at(it)[0] ) );
                                ret.emplace_back ( hlpArr );
                            }
                        }
                        else
                        {
                            for ( int iter = -std::floor( static_cast<int> ( this->dihedralSymmetries.at(i).at(it)[0] ) / 2.0 ); iter <= std::floor( static_cast<int> ( this->dihedralSymmetries.at(i).at(it)[0] ) / 2.0 ); iter++ )
                            {
                                if ( iter == 0 ) { continue; }
                                hlpArr[0]     = this->dihedralSymmetries.at(i).at(it)[0];
                                hlpArr[1]     = this->dihedralSymmetries.at(i).at(it)[1];
                                hlpArr[2]     = this->dihedralSymmetries.at(i).at(it)[2];
                                hlpArr[3]     = this->dihedralSymmetries.at(i).at(it)[3];
                                hlpArr[4]     = iter * ( ( 2.0 * M_PI ) / static_cast<double> ( this->dihedralSymmetries.at(i).at(it)[0] ) );
                                ret.emplace_back ( hlpArr );
                            }
                        }
                    }
                    
                    //======================== Done
                    break;
                }
            }
            
            //================================ Report if required fold was not found
            if ( !foundSym )
            {
                std::cout << "!!! ProSHADE WARNING !!! Requested symmetry elements for " << symType << "-" << symFold << " symmetry, but this symmetry was not found. Returning empty vector." << std::endl;
                return ( std::vector< std::array<double,5> > () );
            }
        }
        
        //==================================== Tetrahedral symmetry
        if ( symType == "T" )
        {
            //================================ If not yet loaded, get them
            if ( this->tetrahedralSymmetry.size() < 1 )
            {
                this->tetrahedralSymmetry     = getTetrahedralSymmetries ( );
            }
            
            //================================ Generate
            if ( this->tetrahedralSymmetry.size() > 1 )
            {
                std::vector< std::array<double,5 >> hlp = symObj.generateTetrElements ( this->tetrahedralSymmetry, settings, -1 );
                for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
                {
                    ret.emplace_back          ( hlp.at(iter) );
                }
            }
            
            //================================ Report if required fold was not found
            if ( this->tetrahedralSymmetry.size() < 1 )
            {
                std::cout << "!!! ProSHADE WARNING !!! Requested symmetry elements for tetrahedral symmetry, but this symmetry was not found. Returning empty vector." << std::endl;
                return ( std::vector< std::array<double,5> > () );
            }
        }
        
        //==================================== Octahedral symmetry
        if ( symType == "O" )
        {
            //================================ If not yet loaded, get them
            if ( this->octahedralSymmetry.size() < 1 )
            {
                this->octahedralSymmetry      = getOctahedralSymmetries ( );
            }
            
            //================================ Generate
            if ( this->octahedralSymmetry.size() > 1 )
            {
                std::vector< std::array<double,5> > hlp = symObj.generateOctaElements ( this->octahedralSymmetry, settings, -1 );
                for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
                {
                    ret.emplace_back          ( hlp.at(iter) );
                }
            }
            
            //================================ Report if required fold was not found
            if ( this->octahedralSymmetry.size() < 1 )
            {
                std::cout << "!!! ProSHADE WARNING !!! Requested symmetry elements for octahedral symmetry, but this symmetry was not found. Returning empty vector." << std::endl;
                return ( std::vector< std::array<double,5> > () );
            }
        }
        
        //==================================== Icosahedral symmetry
        if ( symType == "I" )
        {
            //================================ If not yet loaded, get them
            if ( this->icosahedralSymmetry.size() < 1 )
            {
                this->icosahedralSymmetry     = getIcosahedralSymmetries ( );
            }
            
            //================================ Generate
            if ( this->icosahedralSymmetry.size() > 1 )
            {
                std::vector< std::array<double,5> > hlp = symObj.generateIcosElements ( this->icosahedralSymmetry, settings, -1 );
                for ( unsigned int iter = 1; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
                {
                    ret.emplace_back          ( hlp.at(iter) );
                }
            }
            
            //================================ Report if required fold was not found
            if ( this->icosahedralSymmetry.size() < 1 )
            {
                std::cout << "!!! ProSHADE WARNING !!! Requested symmetry elements for icosahedral symmetry, but this symmetry was not found. Returning empty vector." << std::endl;
                return ( std::vector< std::array<double,5> > () );
            }
        }
        
        if ( ( symType != "I" ) && ( symType != "O" ) && ( symType != "T" ) && ( symType != "D" ) && ( symType != "C" ) )
        {
            std::cerr << "!!! ProSHADE ERROR !!! Unknown symmetry type supplied to the \"getSpecificSymmetryElements\" function. Allowed values are \"C\", \"D\", \"T\", \"O\" and \"I\". Terminating..." << std::endl;
            exit ( -1 );
        }
        
        //==================================== Symmetry element found, now returning
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the symmetry elements list without requesting their computation first. Returning empty vector of arrays." << std::endl;
        return ( std::vector< std::array<double,5> > () );
    }
    
    //======================================== Done
    std::cerr << "!!! ProSHADE ERROR !!! Something went wery wrong... Please report this. Returning empty vector..." << std::endl;
    return ( std::vector< std::array<double,5> > () );
    
}

std::vector< double > ProSHADE::ProSHADE::getSpecificSymmetryElementsPy ( std::string symType,
                                                                          int symFold )
{
    //======================================== Check for sanity
    if ( this->symmetriesAvailable )
    {
        std::vector< std::array< double, 5 > > hlp = this->getSpecificSymmetryElements ( symType, symFold );
        std::vector< double > ret;
        
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( hlp.size() ); iter++ )
        {
            for ( unsigned int it = 0; it < static_cast<unsigned int> ( hlp.at(iter).size() ); it++ )
            {
                ret.emplace_back              ( hlp.at(iter).at(it) );
            }
            ret.emplace_back                  ( -999.999 );
        }
        
        return ( ret );
    }
    else
    {
        std::cout << "!!! ProSHADE WARNING !!! Attempted to access the symmetry elements list without requesting their computation first. Returning empty vector." << std::endl;
        return ( std::vector< double > () );
    }
    
    //======================================== Done
    
}