buildable_monomer.h

Go to the documentation of this file.

// -*- mode: c++; -*-
//
//The Biomolecule Toolkit (BTK) is a C++ library for use in the
//modeling, analysis, and design of biological macromolecules.
//Copyright (C) 2004, Christopher Saunders <ctsa@users.sourceforge.net>
//
//This program is free software; you can redistribute it and/or modify
//it under the terms of the GNU Lesser General Public License as published
//by the Free Software Foundation; either version 2.1 of the License, or (at
//your option) any later version.
//
//This program is distributed in the hope that it will be useful,  but
//WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//Lesser General Public License for more details. 
//
//You should have received a copy of the GNU Lesser General Public License 
//along with this program; if not, write to the Free Software Foundation, 
//Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
//



#ifndef BTK_BUILDABLE_MONOMER_H
#define BTK_BUILDABLE_MONOMER_H


#include "buildable_atom.h"
#include "buildable_atom_algorithms.h"
#include "buildable_atom_group_concept.h"
#include "buildable_monomer_family_pattern.h"
#include "monomer.h"

#include <boost/concept_check.hpp>
#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>

#include <cassert>

#include <string>
#include <utility>


#include <iostream>

namespace BTK{


  template <typename AtomType>
  class BuildableMonomer : public Monomer<AtomType>, public BuildableAtomGroupConcept {
    BOOST_CLASS_REQUIRE2(AtomType, BuildableAtom, boost, ConvertibleConcept);

    typedef Monomer<AtomType> base_type;
    typedef BuildableMonomer<AtomType> self_type;

    typedef BuildableMonomerFamilyPattern pattern_t;
  public:
    typedef self_type buildable_monomer_types;


    typedef typename base_type::value_type value_type;
    typedef typename base_type::reference reference;
    typedef typename base_type::const_reference const_reference;
    typedef typename base_type::pointer pointer;
    typedef typename base_type::const_pointer const_pointer;
    typedef typename base_type::size_type size_type;
    typedef typename base_type::iterator iterator;
    typedef typename base_type::const_iterator const_iterator;


    typedef value_type atom_type;
    typedef iterator atom_iterator;
    typedef const_iterator const_atom_iterator;

    typedef std::map<ATOM::index_t,pointer> atom_index_map_t;

    typedef pattern_t::types::monomer_atom_idx_t monomer_atom_idx_t;

 
    BuildableMonomer(const BuildableMonomerFamilyPattern& family_pattern) 
      : _family_pattern(family_pattern) {}
    
    template <typename AtomIterator>
    BuildableMonomer(AtomIterator first,
                     AtomIterator last,
                     self_type* prev_monomer,
                     self_type* next_monomer,
                     const AtomGroupConcept::group_variant_container_t& modifiers,
                     const BuildableMonomerFamilyPattern& family_pattern) 
      : base_type(first,last), 
        _prev_monomer(prev_monomer), 
        _next_monomer(next_monomer), 
        _group_index(GROUP::name3_to_index(first->group_name3())),
        _group_variants(modifiers),
        _family_pattern(family_pattern)
    {
      set_buildable_group();
      initialize_atoms();
      initialize_bonds();
      if(_prev_monomer) _prev_monomer->set_next_monomer(this);
      if(_next_monomer) _next_monomer->set_prev_monomer(this);
    }
    
    BuildableMonomer(const self_type& source)
      : base_type(source), 
        _prev_monomer(0),
        _next_monomer(0),
        _group_index(source._group_index),
        _group_variants(source._group_variants),
        _family_pattern(source._family_pattern)
    {
      set_buildable_group();
      initialize_bonds();
    }
    
    template <typename T2>
    BuildableMonomer(const BuildableMonomer<T2>& source)
      : base_type(source),
        _prev_monomer(0),
        _next_monomer(0),
        _group_index(source._group_index),
        _group_variants(source._group_variants),
        _family_pattern(source._family_pattern)
    {
      set_buildable_group();
      initialze_bonds();
    }


    self_type& 
    operator=(const self_type& rhs) 
    {
      if (&rhs == this) return *this;
      base_type::operator=(rhs);
      _prev_monomer = 0;
      _next_monomer = 0;
      _group_index = rhs._group_index;
      _group_variants = rhs._group_variants;
      _family_pattern = rhs._family_pattern; 

      // buildable atoms should erase their own bonds when copied
      // so it should be as simple of as reinitializing bonds
      set_buildable_group();
      initialize_bonds();

      return *this;
    }

    template <typename T2>
    self_type& 
    operator=(const BuildableMonomer<T2>& rhs) 
    {
      base_type::operator=(rhs);
      _prev_monomer = 0;
      _next_monomer = 0;
      _group_index = rhs._group_index;
      _family_pattern = rhs._family_pattern;
 
      set_buildable_group();
      initialize_bonds();

      return *this;
    }

    // update atom group interface
    //
    virtual 
    GROUP::index_t 
    group_index() const 
    {
      return _group_index;
    }

    virtual 
    const AtomGroupConcept::group_variant_container_t& 
    group_variants() const 
    {
      return _group_variants;
    }


    virtual 
    const std::string& 
    group_name3() const 
    { 
      return GROUP::index_to_name3(group_index()); 
    }
    
    // initial definition from monomer
    //
    virtual 
    char 
    group_name1() const 
    { 
      return GROUP::index_to_name1(group_index()); 
    }

    //
    void
    insert_group_variant(GROUP_VARIANT::index_t v)
    {
      _group_variants.insert(v);
    }

    void
    erase_group_variant(GROUP_VARIANT::index_t v)
    {
      _group_variants.erase(v);
    }






    reference 
    atom(monomer_atom_idx_t ma) 
    {
      pointer a = atom_pointer(ma);
      assert(a);
      return *a;
    }

    const_reference 
    atom(monomer_atom_idx_t ma) const 
    {
      const_pointer a = atom_pointer(ma);
      assert(a);
      return *a;
    }



    reference 
    atom(ATOM::index_t atom_index) 
    {
      pointer a = atom_pointer(atom_index);
      assert(a);
      return *a;
    }

    const_reference 
    atom(ATOM::index_t atom_index) const 
    {
      const_pointer a = atom_pointer(atom_index);
      assert(a);
      return *a;
    }

    double 
    torsion_angle(TORSION::index_t t) const 
    {
      const pattern_t::types::torsion_atom_t& an = _family_pattern.torsion_atoms(group_index(),t);
      return torsion_angle_impl(an);
    }

    void 
    set_torsion_angle(TORSION::index_t t,
                      double angle) 
    {
      const pattern_t::types::torsion_atom_t& an = _family_pattern.torsion_atoms(group_index(),t);
      double delta_angle = angle - torsion_angle_impl(an);

      reference start    = atom(boost::get<2>(an));
      reference boundary = atom(boost::get<1>(an));

      restricted_rotate(start,boundary,delta_angle);
    }


    // this definition completes buildable atom group concept interface:
    // Doxygen comments are in the concept class
    //  
    virtual 
    bool 
    build_atom(const ATOM::index_t atom_index,
               bool is_bootstrap) 
    {
      // if source atom isn't assigned yet, then it can't be built
      //
      pointer build_atom_ptr = atom_pointer(atom_index);
      if(!build_atom_ptr) return false;

      // try to build from tetrads in this monomer
      //
      bool result = find_tetrads_and_build_atom(std::make_pair(0,atom_index),
                                                build_atom_ptr,
                                                is_bootstrap);
      if(result) return true;

      // if that fails, try to build from tetrads in neighboring monomers
      //
      if(_prev_monomer) {
        result = _prev_monomer->find_tetrads_and_build_atom(std::make_pair(1,atom_index),
                                                            build_atom_ptr,
                                                            is_bootstrap);
        if(result) return true;
      }

      if(_next_monomer) {
        result = _next_monomer->find_tetrads_and_build_atom(std::make_pair(-1,atom_index),
                                                            build_atom_ptr,
                                                            is_bootstrap);
        if(result) return true;
      }

      return false;
    }


  private:

    // internal atom lookup functions --
    //
    // what distinguishes these from the public atom functions is that lookup of an atom
    // that does not exist is not considered an error, this is because the buildable monomer class
    //  needs to lookup atoms that may not exist as a regular part of it's build routines.
    // 
    pointer 
    atom_pointer(monomer_atom_idx_t ma) 
    {
      int monomer_difference = ma.first;
      ATOM::index_t atom_index = ma.second;
      if(monomer_difference==0) {
        return atom_pointer(atom_index);
      } else if(monomer_difference<0) {
        if(! _prev_monomer) return 0;
        else return _prev_monomer->atom_pointer(std::make_pair(monomer_difference+1,atom_index));
      } else {
        if(! _next_monomer) return 0;
        else return _next_monomer->atom_pointer(std::make_pair(monomer_difference-1,atom_index));
      }
    }

    const_pointer 
    atom_pointer(monomer_atom_idx_t ma) const 
    {
      int monomer_difference = ma.first;
      ATOM::index_t atom_index = ma.second;
      if(monomer_difference==0) {
        return atom_pointer(atom_index);
      } else if(monomer_difference<0) {
        if(! _prev_monomer) return 0;
        else return _prev_monomer->atom_pointer(std::make_pair(monomer_difference+1,atom_index));
      } else {
        if(! _next_monomer) return 0;
        else return _next_monomer->atom_pointer(std::make_pair(monomer_difference-1,atom_index));
      }
    }

    pointer 
    atom_pointer(ATOM::index_t atom_index) 
    {
      typename atom_index_map_t::iterator a = _atom_index_to_atom_map.find(atom_index);
      if(a ==_atom_index_to_atom_map.end()) return 0;
      else return a->second;
    }

    const_pointer 
    atom_pointer(ATOM::index_t atom_index) const 
    {
      typename atom_index_map_t::const_iterator a = _atom_index_to_atom_map.find(atom_index);
      if(a ==_atom_index_to_atom_map.end()) return 0;
      else return a->second;
    }


    // given two atom indices, this function tries to find the bond length from the 
    // tetrads, if they are adjacent and if an appropriate tetrad exists. This function
    // was written for the convenience of the bootsrap build procedure
    //
    double
    get_bond_length_from_tetrads(monomer_atom_idx_t ma1,
                                 monomer_atom_idx_t ma2)
    {
      typedef pattern_t::types::tetrad_container_t tc_t;

      tc_t tc;
      double result = -1.;
      
      tc = _family_pattern.tetrads(group_index(),group_variants(),ma1);
      result = get_bond_length_from_tetrad_container(tc,ma2);
 
      if(result >= 0.) return result;

      tc = _family_pattern.tetrads(group_index(),group_variants(),ma2);
      result = get_bond_length_from_tetrad_container(tc,ma1);
      
      return result;
    }

    //  just grab a bond length from tetrads, used during bootstrap
    // build
    //
    double
    get_bond_length_from_tetrad_container(const pattern_t::types::tetrad_container_t& tc,
                                          const monomer_atom_idx_t& ma)
    {
      pattern_t::types::tetrad_container_t::const_iterator ti=tc.begin(),ti_end=tc.end();
      for(; ti != ti_end ; ++ti ){   
        if(ti->atom_index3()==ma) return ti->len34();
      }
      return -1.;
    }

    // convenient aggragate function for the build_atom procedure
    //
    bool
    find_tetrads_and_build_atom(monomer_atom_idx_t ma,
                                pointer build_atom_ptr,
                                bool is_bootstrap)
    {
      typedef pattern_t::types::tetrad_container_t tc_t;
      tc_t tc = _family_pattern.tetrads(group_index(),group_variants(), ma);

      return build_atom_from_tetrad_container(tc,build_atom_ptr,is_bootstrap);
    }
   
    // search through tetrads and try to build the build_atom from them
    //
    bool
    build_atom_from_tetrad_container(const pattern_t::types::tetrad_container_t& tc,
                                     pointer build_atom_ptr,
                                     bool is_bootstrap)
    {
      // for the build to work right, we must use buildable improper
      // torsion tetrads (pass 0) before buildable real torsion tetrads (pass 1)
      //
      for(unsigned pass=0;pass<2;++pass){
        pattern_t::types::tetrad_container_t::const_iterator ti=tc.begin(),ti_end=tc.end();
        for(; ti != ti_end ; ++ti ){
          if((pass == 0) && (ti->is_real_torsion())) continue;
          if((pass == 1) && (!ti->is_real_torsion())) continue;
          bool result = build_atom_from_tetrad(*ti,build_atom_ptr,is_bootstrap);
          if (result) return true;
        }
      }
      return false;
    }


    // try to build atom from a tetrad 
    //
    bool
    build_atom_from_tetrad(const AtomIndexTetrad& t,
                           pointer build_atom_ptr,
                           bool is_bootstrap)
    {
      // can a build be made from this tetrad?
      // check: 1) tetrad atoms are assigned and
      //        2) are built
      //
      const_pointer atom_ptr3 = atom_pointer(t.atom_index3());
      if(!atom_ptr3) return false;
      if(!is_bootstrap && !atom_ptr3->is_built()) return false;
      const_pointer atom_ptr2 = atom_pointer(t.atom_index2());
      if(!atom_ptr2) return false;
      if(!is_bootstrap && !atom_ptr2->is_built()) return false;
      
      if(!is_bootstrap){
        const_pointer atom_ptr1 = atom_pointer(t.atom_index1());
        if(!atom_ptr1) return false;
        if(!atom_ptr1->is_built()) return false;

        // regular tetrad build
        //
        build_atom_ptr->
          set_position(set_vector_from_dihedral(atom_ptr3->position(),
            atom_ptr2->position(),
            atom_ptr1->position(),
            t.len34(),
            t.ang234(),
            t.dih1234()));
        build_atom_ptr->set_is_built(true);

      } else {

        // build atoms 4,3,2 from scratch under length and angle constraints -
        // this should only occur to bootstrap a build of a completely new 
        // molecule
        //
        double len23 = get_bond_length_from_tetrads(t.atom_index3(),t.atom_index2());
        if( len23 <= 0.) return false;
        build_atom_ptr->set_position(new_vector(t.len34(),0.,0.));
        build_atom_ptr->set_is_built(true);
        atom_ptr3->set_position(new_vector(0.,0.,0.));
        atom_ptr3->set_is_built(true);
        atom_ptr2->set_position(new_vector(cos(t.ang234())*len23,sin(t.ang234())*len23,0.));
        atom_ptr2->set_is_built(true);
      }
      return true;
    }


    // update group pointers in buildable_atoms:
    void 
    set_buildable_group()
    {
      iterator i=begin(),i_end=end();
      for(;i!=i_end;++i) i->set_buildable_group(*this);
    }

    double 
    torsion_angle_impl(const pattern_t::types::torsion_atom_t& an) const 
    {
      return point_dihedral( atom(boost::get<0>(an)).position(),
                             atom(boost::get<1>(an)).position(),
                             atom(boost::get<2>(an)).position(),
                             atom(boost::get<3>(an)).position());
    }

    void 
    initialize_bonds() 
    {
      typedef pattern_t::types::bond_container_t bc_t;  
      bc_t bc = _family_pattern.bonds(group_index(),group_variants());

      std::cerr << "init bonds: group_num: " << group_num() << std::endl;

      bc_t::const_iterator bi=bc.begin(),bi_end=bc.end();
      for(; bi != bi_end; ++bi){
        pointer a1 = atom_pointer(bi->atom_index1());
        if(!a1) continue;
        pointer a2 = atom_pointer(bi->atom_index2());
        if(!a2) continue;
        std::cerr << "creating bond: atom1: " << *a1 << std::endl;
        std::cerr << "creating bond: atom2: " << *a2 << std::endl;

        add_bond(*a1,*a2);
      }
    }

    void 
    set_next_monomer(self_type* next_monomer)
    {
      _next_monomer = next_monomer;
    }

    void 
    set_prev_monomer(self_type* prev_monomer)
    {
      _prev_monomer = prev_monomer;
    }

    void 
    initialize_atoms() 
    {
      pattern_t::types::atom_container_t ac = _family_pattern.atoms(group_index(),group_variants());

      std::cerr << "init atoms: group_num: " << group_num() << std::endl;

      // assign the atoms that are expected
      for(atom_iterator i=atom_begin();i != atom_end(); ++i){
        ATOM::index_t atom_index = i->atom_index();

        if( ac.find(atom_index) == ac.end() ){
          // throwing away input atom -- add some notification options here
          std::cerr << "throwing away atom: " << i->atom_name() << std::endl;
          i = erase(i) - 1;
        } else {
          // allow atom in monomer
          std::cerr << "adding atom: " << i->atom_name() << std::endl;
          _atom_index_to_atom_map[ atom_index ] = &(*i);
          ac.erase(atom_index);

          i->set_is_built(true);
        }
      }

      // add missing atoms
      pattern_t::types::atom_container_t::const_iterator i=ac.begin(), i_end=ac.end();
      for(;i!=i_end;++i){
        shallow_push_back(new atom_type());
        back().set_atom_index(*i);
        _atom_index_to_atom_map[ *i ] = &(back());

        // new atoms are not built yet
        back().set_is_built(false);
        back().set_group(*this);
        back().set_buildable_group(*this);
        std::cerr << "adding new atom: " << back().atom_name() << std::endl;
      }
    }


  private:
    self_type* _prev_monomer;
    self_type* _next_monomer;

    atom_index_map_t _atom_index_to_atom_map;

    GROUP::index_t _group_index;

    AtomGroupConcept::group_variant_container_t _group_variants;

    const pattern_t& _family_pattern;

  };
} // namespace BTK



#endif

---