Containers.h

Go to the documentation of this file.

 
#ifndef PNAB_CONTAINERS_H
#define PNAB_CONTAINERS_H
 
#include <string>
#include <random>
#include <array>
#include <openbabel/mol.h>
#include <openbabel/atom.h>
#include <openbabel/bond.h>
#include <openbabel/obconversion.h>
#include <openbabel/math/matrix3x3.h>
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
#ifndef DEG_TO_RAD
#define DEG_TO_RAD (M_PI/180.0)
#endif
 
namespace PNAB {
    class RuntimeParameters {
 
    public:
        RuntimeParameters() : energy_filter{}, max_distance(), ff_type(), glycosidic_bond_distance(0.0),
                              num_steps(0), seed(0), weighting_temperature(298.0), monte_carlo_temperature(298.0),
                              mutation_rate(0.75), crossover_rate(0.75), population_size(1000), strand{}, is_hexad(false),
                              build_strand(std::vector<bool> {true, false, false, false, false, false}), num_candidates(10),
                              strand_orientation(std::vector<bool> {true, true, true, true, true, true}){};
 
        // Thresholds
        std::vector<double> energy_filter;      
        double max_distance;                    
        // Force Field Parameter
        std::string ff_type;                    
        // Algorithm parameters
        std::string search_algorithm;           
        std::size_t num_steps;                  
        unsigned int seed;                      /*< @brief The seed for the random number generator. Use the same value to get reproducible
                                                * results.
                                                *
                                                *  @sa ConformationSearch::MonteCarloSearch
                                                *  @sa ConformationSearch::RandomSearch
                                                *  @sa ConformationSearch::GeneticAlgorithmSearch
                                                */
 
        double dihedral_step;                   
        double weighting_temperature;           
        double monte_carlo_temperature;         
        double mutation_rate;                   
        double crossover_rate;                  
        int population_size;                    
        //Strand parameters
        std::vector<std::string> strand;        
        std::vector<bool> build_strand;         
        std::vector<bool> strand_orientation;   
        bool is_hexad;                          
        // Glycosidic bond
        double glycosidic_bond_distance;        
        unsigned int num_candidates;            
    };
 
    class HelicalParameters {
 
    public:
        HelicalParameters() : h_twist{0}, h_rise{0}, inclination{0}, tip{0}, x_displacement{0}, y_displacement{0},
                              shift{0}, slide{0}, rise{0}, tilt{0}, roll{0}, twist{0},
                              buckle{0}, propeller{0}, opening{0}, shear{0}, stretch{0}, stagger{0}, is_helical{true} {};
 
        //Helical parameters
        double     inclination,                         
                   tip,                                 
                   h_twist,                             
                   x_displacement,                      
                   y_displacement,                      
                   h_rise,                              
                   shift,                               
                   slide,                               
                   rise,                                
                   tilt,                                
                   roll,                                
                   twist,                               
                   buckle,                              
                   propeller,                           
                   opening,                             
                   shear,                               
                   stretch,                             
                   stagger;                             
        bool is_helical;                                
 
        void computeHelicalParameters();
 
        OpenBabel::vector3 getGlobalTranslationVec(bool is_base_pair=false, bool is_second_strand=false);
 
        OpenBabel::vector3 getStepTranslationVec(unsigned n = 0, bool is_base_pair=false, bool is_second_strand=false);
 
        // Lu, X. J., El Hassan, M. A., & Hunter, C. A. (1997). Structure and conformation of helical nucleic acids:
        // rebuilding program (SCHNArP). Journal of molecular biology, 273(3), 681-691.
        OpenBabel::matrix3x3 getGlobalRotationMatrix(bool is_base_pair=false, bool is_second_strand=false);
 
        OpenBabel::matrix3x3 getStepRotationMatrix(unsigned n = 0, bool is_base_pair=false, bool is_second_strand=false);
 
    private:
        OpenBabel::matrix3x3 rodrigues_formula(OpenBabel::vector3 axis_vector, double theta);
 
        std::vector<OpenBabel::vector3> StepParametersToReferenceFrame();
 
        void ReferenceFrameToHelicalParameters(OpenBabel::vector3 origin2, OpenBabel::vector3 x2, OpenBabel::vector3 y2, OpenBabel::vector3 z2);
 
    };
 
    class Backbone {
    public:
 
        Backbone() : backbone{}, file_path{}, interconnects{}, linker{}, vector_atom_deleted{}, fixed_bonds{} {}
 
        Backbone(std::string file_path, std::array<unsigned, 2> interconnects, std::array<unsigned,2> linker, std::vector<std::vector<unsigned>> fixed_bonds = {});
 
        OpenBabel::OBAtom* getHead() {
            return backbone.GetAtom(interconnects[0]);
        }
 
        OpenBabel::OBAtom* getTail() {
            return backbone.GetAtom(interconnects[1]);
        }
 
        OpenBabel::OBAtom* getLinker() {
            return backbone.GetAtom(linker[0]);
        }
 
        OpenBabel::OBAtom* getVector() {
            if (!vector_atom_deleted)
                return backbone.GetAtom(linker[1]);
            else {
                std::cerr << "Called getVector() for backbone with no vector atom." << std::endl;
                return nullptr;
            }
        }
 
        void center() {
            backbone.Center();
        }
 
        void rotate(double* rot) {
            backbone.Rotate(rot);
        }
 
        void translate(OpenBabel::vector3 vec) {
            backbone.Translate(vec);
        }
 
        OpenBabel::OBMol getMolecule() {
            return backbone;
        }
 
        void deleteVectorAtom();
 
        std::array<unsigned , 2> interconnects,          
                                 linker;                 
        std::vector<std::vector<unsigned>> fixed_bonds;  
        OpenBabel::OBMol backbone;                       
        std::string file_path;                           
 
    private:
        void validate();
 
        bool vector_atom_deleted;                        
    };
 
    class Base {
 
    public:
 
        Base() : name{}, code{}, linker{}, base{}, vector_atom_deleted{}, pair_name{} {};
 
        Base(std::string name, std::string code, std::string file_path, std::array<std::size_t, 2> linker,
             std::string pair_name = "");
 
        OpenBabel::OBAtom* getLinker() {
            return base.GetAtom(static_cast<unsigned>(linker[0]));
        }
 
        OpenBabel::OBAtom* getVector() {
            if (!vector_atom_deleted)
                return base.GetAtom(static_cast<unsigned>(linker[1]));
            else
                return nullptr;
        }
 
        OpenBabel::OBMol getMolecule() {
            return base;
        }
 
        std::string getCode() {
            return code;
        }
 
        std::string getName() {
            return name;
        }
 
        void deleteVectorAtom() {
            if (!vector_atom_deleted) {
                size_t id = getLinker()->GetId();
                base.DeleteAtom(getVector());
                linker = {base.GetAtomById(id)->GetIdx(), 0};
                vector_atom_deleted = true;
            }
        }
 
        std::string getBasePairName() {
            return pair_name;
        }
 
        std::string name,                               
                    code,                               
                    pair_name,                          
                    file_path;                          
        OpenBabel::OBMol base;                          
        std::array<std::size_t, 2 > linker;             
 
    private:
        void validate();
 
        bool vector_atom_deleted;                       
    };
 
    class Bases {
    public:
        Bases(std::vector<Base> input_bases);
 
        Bases() {};
 
        PNAB::Base getBaseFromName(std::string name) {
            for (auto v : bases) {
                if (v.getName().find(name) != std::string::npos)
                    return v;
            }
            std::cerr << "Base \"" << name << "\" does not exists in list of bases. Please check input file."
                      << std::endl;
            throw 1;
        }
 
        std::vector<Base> getBasesFromStrand(std::vector<std::string> strand);
 
        std::vector<Base> getComplimentBasesFromStrand(std::vector<std::string> strand);
 
    private:
        std::vector<Base> bases;                          
        bool all_bases_pair;                              
        std::map<std::string, PNAB::Base> name_base_map;  
    };
 
    class BaseUnit {
    public:
        BaseUnit(Base b, Backbone backbone, double glycosidic_bond_distance = 0.0);
 
        BaseUnit() {};
 
        const OpenBabel::OBMol getMol() {
            return unit;
        }
 
        const std::array< std::size_t, 2 > getBaseIndexRange() {
            return base_index_range;
        };
 
        const std::array< std::size_t, 2 > getBackboneIndexRange() {
            return backbone_index_range;
        };
 
        const std::array< std::size_t, 2 > getBackboneLinkers() {
            return backbone_interconnects;
        };
 
        std::size_t getBaseConnectIndex() {
            return base_connect_index;
        }
 
        std::vector<std::vector<unsigned>> getFixedBonds() {
            return fixed_bonds;
        }
 
    private:
        OpenBabel::OBMol unit;                                                       
        std::array< std::size_t, 2 > base_index_range,                               
                                     backbone_index_range;                           
        std::size_t base_connect_index;                                              
        std::array< std::size_t, 2 > backbone_interconnects;                         
        std::vector<std::vector<unsigned>> fixed_bonds;                              
    };
 
    struct ConformerData {
        OpenBabel::OBMol molecule;                  
        double *monomer_coord,                      
                distance,                           
                bondE,                              
                angleE,                             
                torsionE,                           
                VDWE,                               
                total_energy,                       
                rmsd;                               
        std::size_t index;                          
        bool accepted;                              
        std::vector<double> dihedral_angles;        
 
        bool operator<(const ConformerData &cd) const {
            return (total_energy < cd.total_energy);
        }
    };
 
}
 
#endif //PNAB_CONTAINERS_H