/***************************************************************************
                          molecules_mod.cpp  -  description
                             -------------------
    begin                : Sat Mar 22 01:02:37 2003
    copyright            : (C) 2002 by Cavalli Andrea
    email                : cavalli@bioc.unizh.ch
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <molecules/molecules.h>
#include <xmolecules.h>
#include <mdb/mdb.h>
#include <almostexception.h>
#include <almosttypes.h>
#include <almosttemplate.h>

#include <zeta/module.h>

#include <coor.h>
#include <molecules/coor3d.h>
#include <select.h>
#include <kernel/rmsd.h>
#include <moveset/bondrotation.h>
#include <secstruct.h>
#include <algorithm/maxsub.h>
#include <algorithm/molalgo.h>
#include <algorithm/mkbetapair.h>

#include <molecules/atomz.h>
#include <molecules/coor3dz.h>
#include <molecules/proteinz.h>
#include <molecules/nucleicacidz.h>
#include <molecules/lipidz.h>

using namespace Almost;

void auto_connect(Molecules & mol){
  //Find d(i,j)<CovRadoous;
  //  double cov_radious[ =
  double C_cov = 0.75;
  double N_cov = 0.7;
  double O_cov = 0.65;
  double H_cov = 0.45;
  double P_cov = 0.9;
  double S_cov = 0.95;
  double tol = 0.4;
  for(int i=0;i<mol.atom_size();i++){
    double a1_cov = 0;
    if(mol[i].name()[0] == 'C') a1_cov = C_cov;
    else if(mol[i].name()[0] == 'N') a1_cov = N_cov;
    else if(mol[i].name()[0] == 'O') a1_cov = O_cov;
    else if(mol[i].name()[0] == 'H') a1_cov = H_cov;
    else if(mol[i].name()[0] == 'P') a1_cov = P_cov;
    else if(mol[i].name()[0] == 'S') a1_cov = S_cov;
    else { 
      cout<<"Unable to find covalent radious for : "<<mol[i].name()<<endl;
      throw;
    }
    for(int j=i+1;j<mol.atom_size();j++){
      double a2_cov = 0;
      if(mol[j].name()[0] == 'C') a2_cov = C_cov;
      else if(mol[j].name()[0] == 'N') a2_cov = N_cov;
      else if(mol[j].name()[0] == 'O') a2_cov = O_cov;
      else if(mol[j].name()[0] == 'H') a2_cov = H_cov;
      else if(mol[j].name()[0] == 'P') a2_cov = P_cov;
      else if(mol[j].name()[0] == 'S') a2_cov = S_cov;
      else { 
	cout<<"Unable to find covalent radious for : "<<mol[j].name()<<endl;
	throw;
      }

      double d = mol.distance(i,j);
      if(d<(a1_cov+a2_cov)+tol){
	cout<<">> Putative bond "<<mol.atom_name(i,Molecules::BASE)<<" "<<mol.atom_name(j,Molecules::BASE)
	    <<" "<<d<<" "<<a1_cov+a2_cov<<"\n";
      }
    }
  }

}


Protein auto_connect_pdb(PDBSegment & segment){
  Protein p("PDB");
  {
    vector<string> sequ = segment.sequence();
    for(int i=0;i<sequ.size();i++){
     int fragd =  p.add_fragment(sequ[i]);
     for(int j=0;j<segment[i].size();j++){
       Atom atom(segment[i][j].name,"UKN",
		 segment[i][j].x,segment[i][j].y,segment[i][j].z,
		 -1,
		 -1);
       p.add_atom(atom,fragd);
     }
    }
  }
  //Add segment frags
  
  //Find d(i,j)<CovRadoous;
  //  double cov_radious[ =
  double C_cov = 0.75;
  double N_cov = 0.7;
  double O_cov = 0.65;
  double H_cov = 0.45;
  double P_cov = 0.9;
  double S_cov = 0.95;
  double tol = 0.4;
  for(int i=0;i<p.atom_size();i++){
    double a1_cov = 0;
    if(p[i].name()[0] == 'C') a1_cov = C_cov;
    else if(p[i].name()[0] == 'N') a1_cov = N_cov;
    else if(p[i].name()[0] == 'O') a1_cov = O_cov;
    else if(p[i].name()[0] == 'H') a1_cov = H_cov;
    else if(p[i].name()[0] == 'P') a1_cov = P_cov;
    else if(p[i].name()[0] == 'S') a1_cov = S_cov;
    else { 
      cout<<"Unable to find covalent radious for : "<<p[i].name()<<endl;
      throw;
    }
    for(int j=i+1;j<p.atom_size();j++){
      double a2_cov = 0;
      if(p[j].name()[0] == 'C') a2_cov = C_cov;
      else if(p[j].name()[0] == 'N') a2_cov = N_cov;
      else if(p[j].name()[0] == 'O') a2_cov = O_cov;
      else if(p[j].name()[0] == 'H') a2_cov = H_cov;
      else if(p[j].name()[0] == 'P') a2_cov = P_cov;
      else if(p[j].name()[0] == 'S') a2_cov = S_cov;
      else { 
	cout<<"Unable to find covalent radious for : "<<p[j].name()<<endl;
	throw;
      }

      double d = p.distance(i,j);
      if(d<(a1_cov+a2_cov)+tol){
	cout<<">> Putative bond "<<p.atom_name(i,Protein::BASE)<<" "<<p.atom_name(j,Protein::BASE)
	    <<" "<<d<<" "<<a1_cov+a2_cov<<"\n";
	if(!p.add_bond(i,j,0)){
	  aout<<"\t>> Unable to add bond ....\n";
	}
      }
    }
  }
  return p;
}

template<>
inline string to_string(const Almost::DCDTrj & ){
  return "<Almost::DCDTrj>";
}


template<>
inline string to_string(const Almost::Atom & ){
  return "<Almost::Atom>";
}


template<>
inline string to_string(const Kernel::Molecule<Almost::Atom, int> & ){
  return "<Kernel::Molecule<Almost::Atom, int>>";
}

template<>
inline string to_string(const Kernel::Molecule<Almost::Atom, Almost::MoleculeBond> & ){
  return "<Kernel::Molecule<Almost::Atom, Almost::MoleculeBond>>";
}


template<>
inline string to_string(const Protein & protein){
  return "<Protein>";
}

template<>
inline string to_string(const NucleicAcid & nucleicacid){
  return "<NucleicAcid>";
}

template<>
inline string to_string(const Lipid & lipid){
  return "<Lipid>";
}

template<>
inline string to_string(const Solvent & solvent){
  return "<Solvent>";
}

template<>
inline string to_string(const Molecules & molecules){
  return "<Molecules>";
}

template<>
inline string to_string(const XMolecules & molecules){
  return "<XMolecules>";
}

template<>
inline string to_string(const PDBSegment & pdb){
  return "<PDB>";
}

template<>
inline string to_string(const MDB & mdb){
  return "<MDB>";
}

template<>
inline string to_string(const DSSPSecStruct & ){
  return "<DSSPSecStruct>";
}

template<>
inline string to_string(const AlmSecStruct & ){
  return "<AlmSecStruct>";
}

template<>
inline string to_string(const MakeBetaPair & ){
  return "<MakeBetaPair>";
}


#define ALMTOSTRING(s) \
template<> \
inline string to_string(const s & ){ \
  return "<" #s ">"; \
}

ALMTOSTRING(Almost::BondRotation<double>);
ALMTOSTRING(Almost::NameSelection);

template<>
inline string to_string( const Coor3D & coor){
  string str;
  str = "[ ";
  for(int i=0;i<3;i++){
    str += to_string(coor[i]);
    if(i!=2) str +=", ";
  }

  str += " ]";
  return str;
}


void prot2alm(const Protein & prot, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  FormOStream fout(file);
  fout<<prot<<endl;
}

void na2alm(const NucleicAcid & na, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  FormOStream fout(file);
  fout<<na<<endl;
}

void lipid2alm(const Lipid & prot, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  FormOStream fout(file);
  fout<<prot<<endl;
}

void solv2alm(const Solvent & solv, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  FormOStream fout(file);
  fout<<solv<<endl;
}


void prot2pdb(const Protein & prot, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  PDBOStream  fout(file);
  fout<<prot;;
}

void na2pdb(const NucleicAcid & prot, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  PDBOStream  fout(file);
  fout<<prot;;
}

void lipid2pdb(const Lipid & lip, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  PDBOStream  fout(file);
  fout<<lip;;
}

void solv2pdb(const Solvent & solv, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  PDBOStream fout(file);
  fout<<solv;
}

void mol2pdb(const Molecules & mols, string ofile){
  ofstream file;
  file.open(ofile.c_str());
  PDBOStream  fout(file);
  fout<<mols;;
}

Coor<double> mols_coor(const Molecules & molecules){
  Coor<double> coor;
  coor.resize(molecules.atom_size());
  for(int i=0; i<molecules.atom_size();i++){
    int pos = Coor<double>::DIM*i;
    coor.coor[pos] = molecules[i][0];
    coor.coor[pos+1] = molecules[i][1];
    coor.coor[pos+2] = molecules[i][2];
  }
  return coor;
}

template<typename float_type>
inline string to_string(const Coor<float_type> & ){
  return "<Coor>"; 
}

ostream & operator<<(ostream & out,const Coor<double> & ){
  out<<"<Coor>";
  return out;
}

template<>
inline string to_string(const Coor<double> & ){
  return "<Coor>"; 
}

Coor<double> mol_coor(const Molecules & m){
  return Coor<double>(m);
}



void check_bonds(const Molecules & molecules, const MDB & mdb){
      vector<int> bond_first_;
      vector<int> bond_second_;
      vector<double> b0__;
      for(int i=0;i<molecules.atom_size();i++){
	vector<int> bonds = molecules.bonds(i);
	string name1 = molecules[i].type_name();
	for(int j=0;j<bonds.size();j++){
	  if(bonds[j]<i) continue;
	  string name2 = molecules[bonds[j]].type_name();
	  BondT bt = mdb.bondT(name1,name2);
	  bond_first_.push_back(i);
	  bond_second_.push_back(bonds[j]);
	  b0__.push_back(bt.equ_length());
	}
      }

      for(int i=0;i<bond_first_.size();i++){
	double d = molecules.distance(bond_first_[i],bond_second_[i]);
	if(d>2*b0__[i]){
	  cout<<"\t>> Bond "<<bond_first_[i]<<" "<<molecules[bond_first_[i]].name()<<" "
	      <<bond_second_[i]<<" "<<molecules[bond_second_[i]].name()<<" "
	      <<"to long "<<d<<" "<<b0__[i]<<"\n";
	}
	
      }

}

template<>
inline string to_string(const MaxSub & ){
  return "<MaxSub>"; 
}

namespace Almost {
  void maxsub(const Molecules & m1, const Molecules & m2, double d,
	      double & rmsd,
	      int & size){
    MaxSub ms(m1,m2,d);
    rmsd = ms.rmsd();
    size = ms.size();
  }

  MaxSub maxsub_const(const Molecules & m1,  double d){
    return MaxSub(m1,d);
  }
}


Protein mutate(Protein & p,int aa, string f,MDB & mdb,string fp, string lp){
  return p.mutate(aa,f,mdb,fp,lp);
}


void mol2psf(const Molecules & mol, MDB & mdb, string out_file){
  ofstream out;
  out.open(out_file.c_str());
  if(!out){
    throw (">> Unable to open file: "+out_file).c_str();
  }

  char buff[1024];
  int nrem = 3;
  sprintf(buff,"PSF\n\n%8d !NTITLE\n",nrem);
  out<<buff;
  //Remarks
  out<<" REMARKS PSF generated with almost, what else?\n";
  out<<" REMARKS (c) 1999-2009 Andrea Cavalli\n";
  out<<" REMARKS amc82@cam.ac.uk, http://www.open-almost.org\n\n";

  //Atoms
  int natoms = mol.atom_size();
  sprintf(buff,"%8d !NATOM\n",natoms);
  out<<buff;
  for(int i=0;i<natoms;i++){
    int fd = mol.find_fragment(i);
    sprintf(buff,"%8d %-4s %-4d %-4s %-4s %-4s %10.6f     %9.4f  %10d\n",
	    i+1,
	    "MAIN",
	    fd+1,
	    mol.fragment_name(fd,Molecules::SHORT).c_str(), 
	    mol[i].name().c_str(),
	    mol[i].type_name().c_str(),
	    mol[i].charge(),
	    mol[i].mass(),
	    0);
    out<<buff;
  }
  out<<"\n";
  {
    // NBOND
    map<int,int> bonds;
    for(int i=0;i<natoms;i++){
      const vector<int> & b = mol.bonds(i);
      for(int j=0;j<b.size();j++){
	if(i<b[j])
	  bonds[i] = b[j];
      }
    }
    sprintf(buff,"%8lu !NBOND: bonds\n",bonds.size());
    out<<buff;
    int count = 0;
    map<int,int>::iterator iter = bonds.begin(),end = bonds.end();
    while(iter!=end){
      if(count==4){ out<<"\n";count = 0;}
      sprintf(buff," %7d %7d",iter->first+1,iter->second+1);
      out<<buff;
      ++count;
      ++iter;
    }
    out<<"\n\n";
  }
  {
    //Angles
    vector<vector<int> > angles;
    for(int i=0;i<mol.protein_size();i++){
      vector<vector<int> > aa = mol.protein(i).angles();
      for(int a = 0; a<aa.size();a++)
	angles.push_back(aa[a]);
    }

    for(int i=0;i<mol.nucleicacid_size();i++){
      vector<vector<int> > aa = mol.nucleicacid(i).angles();
      for(int a = 0; a<aa.size();a++)
	angles.push_back(aa[a]);
    }

    for(int i=0;i<mol.lipid_size();i++){
      vector<vector<int> > aa = mol.lipid(i).angles();
      for(int a = 0; a<aa.size();a++)
	angles.push_back(aa[a]);
    }

    for(int i=0;i<mol.solvent_size();i++){
      vector<vector<int> > aa = mol.solvent(i).angles();
      for(int a = 0; a<aa.size();a++)
	angles.push_back(aa[a]);
    }

    sprintf(buff,"%8lu !NTHETA: angles\n",angles.size());
    out<<buff;
    int count = 0;
    for(int i=0;i<angles.size();i++){
      if(count==3){ out<<"\n"; count=0;}
      sprintf(buff," %7d %7d %7d",
	      angles[i][0]+1,
	      angles[i][1]+1,
	      angles[i][2]+1);
      out<<buff;
      ++count;
    }
    out<<"\n\n";
  }

  {
    //Dihedrals
    vector<vector<int> > dihedrals;
    for(int i=0;i<mol.protein_size();i++){
      vector<vector<int> > aa = mol.protein(i).dihedrals();
      for(int a = 0; a<aa.size();a++)
	dihedrals.push_back(aa[a]);
    }

    for(int i=0;i<mol.nucleicacid_size();i++){
      vector<vector<int> > aa = mol.nucleicacid(i).dihedrals();
      for(int a = 0; a<aa.size();a++)
	dihedrals.push_back(aa[a]);
    }

    for(int i=0;i<mol.lipid_size();i++){
      vector<vector<int> > aa = mol.lipid(i).dihedrals();
      for(int a = 0; a<aa.size();a++)
	dihedrals.push_back(aa[a]);
    }

    for(int i=0;i<mol.solvent_size();i++){
      vector<vector<int> > aa = mol.solvent(i).dihedrals();
      for(int a = 0; a<aa.size();a++)
	dihedrals.push_back(aa[a]);
    }

    sprintf(buff,"%8lu !NPHI: dihedrals\n",dihedrals.size());
    out<<buff;
    int count = 0;
    for(int i=0;i<dihedrals.size();i++){
      if(count==3){ out<<"\n"; count=0;}
      sprintf(buff," %7d %7d %7d %7d",
	      dihedrals[i][0]+1,
	      dihedrals[i][1]+1,
	      dihedrals[i][2]+1,
	      dihedrals[i][3]+1);
      out<<buff;
      ++count;
    }
    out<<"\n\n";
  }

  {
    //Impropers
    vector<vector<int> > impropers;
    for(int i=0;i<mol.protein_size();i++){
      vector<vector<int> > aa = mol.protein(i).impropers();
      for(int a = 0; a<aa.size();a++)
	impropers.push_back(aa[a]);
    }

    for(int i=0;i<mol.nucleicacid_size();i++){
      vector<vector<int> > aa = mol.nucleicacid(i).impropers();
      for(int a = 0; a<aa.size();a++)
	impropers.push_back(aa[a]);
    }

    for(int i=0;i<mol.lipid_size();i++){
      vector<vector<int> > aa = mol.lipid(i).impropers();
      for(int a = 0; a<aa.size();a++)
	impropers.push_back(aa[a]);
    }

    for(int i=0;i<mol.solvent_size();i++){
      vector<vector<int> > aa = mol.solvent(i).impropers();
      for(int a = 0; a<aa.size();a++)
	impropers.push_back(aa[a]);
    }

    sprintf(buff,"%8lu !NIMPHI: impropers\n",impropers.size());
    out<<buff;
    int count = 0;
    for(int i=0;i<impropers.size();i++){
      if(count==3){ out<<"\n"; count=0;}
      sprintf(buff," %7d %7d %7d %7d",
	      impropers[i][0]+1,
	      impropers[i][1]+1,
	      impropers[i][2]+1,
	      impropers[i][3]+1);
      out<<buff;
      ++count;
    }
    out<<"\n\n";
  }

  sprintf(buff,"\n\n");
  out<<buff;
  sprintf(buff,"%8d !NDON: donors\n",0);
  out<<buff;
  sprintf(buff,"\n\n");
  out<<buff;
  sprintf(buff,"%8d !NACC: acceptors\n",0);
  out<<buff;
  sprintf(buff,"\n\n");
  out<<buff;


  {
    sprintf(buff,"%8d !NNB\n\n",0);
    out<<buff;
    int i, fullrows;
    fullrows = natoms/8;
    for (i=0; i<fullrows; ++i) {
      sprintf(buff, "%8d%8d%8d%8d%8d%8d%8d%8d\n",0,0,0,0,0,0,0,0);
      out<<buff;
    }
    for (i=natoms - fullrows*8; i; --i){
      sprintf(buff, "%8d",0);
      out<<buff;
    }
  } 
  sprintf(buff,"\n\n");
  out<<buff;

  sprintf(buff,"%8d %7d !NGRP\n%8d%8d%8d\n",1,0,0,0,0);
  out<<buff;
  sprintf(buff,"\n");
  out<<buff;
}



// template< typename Base, typename Derived>

// string (Derived::*to_derived(string (Base::*f)() const))() const{
//   return f;
// }



void init_molecules_types(ZModule & module){
  init_atom_class(module);
  init_coor3d_class(module);
  init_protein_class(module);
  init_nucleicacid_class(module);
  init_lipid_class(module);

  typedef void (Coor<double>::*p_copy)(const Coor<double> &);
  Class<Coor<double> >(module,"coor")
    .def_method("size",&Coor<double>::size)
    .def_method("resize",&Coor<double>::resize)
    .def_method("clear",&Coor<double>::clear)
    .def_method("X",&Coor<double>::X_coor)
    .def_method("Y",&Coor<double>::Y_coor)
    .def_method("Z",&Coor<double>::Z_coor)
    .def_method<p_copy>("copy",&Coor<double>::copy);




  Class<MaxSub,ArgList3<const Molecules & ,
    const Molecules & , double> >(module,"maxsub")
    .def_method("rmsd",&MaxSub::rmsd)
    .def_method("size",&MaxSub::size)
    .def_method("alignment",&MaxSub::alignment)
    .def_method("score",&MaxSub::score)
    .def_method("apply_dcd",&MaxSub::apply_dcd)
    ;



  Class<Solvent,ArgList1<string> >(module,"solvent")
    //    .def_method("name",&Solvent::name)
    //    .def_method("set_name",&Solvent::set_name)
    .def_method("chainID",&Solvent::chainID)
    .def_method("set_chainID",&Solvent::set_chainID)
    .def_method("find_atom",&Solvent::find_atom)
    .def_method("build",&Solvent::build)
    .def_method("build_sequ",&Solvent::build_sequ)
    .def_method("translate",&Solvent::translate)
    .def_method("rotate",&Solvent::rotate)
    .def_method("cog",&Solvent::cog)
    .def_method("com",&Solvent::com)
    .def_method<string (Solvent::*)() const>("name",&Solvent::name)
    .def_method<int (Solvent::*)() const>("atom_size",&Solvent::atom_size)
    .def_doc("atom_size","Returns the number of atoms in the solvent.")
    ;


 Class<Molecules>(module,"molecules")
    .def_method("add_solvent",&Molecules::add_solvent)
    .def_doc("add_solvent","Adds the solvent to the molecules collection.")
    .def_method("add_lipid",&Molecules::add_lipid)
    .def_doc("add_lipid","Adds the lipids to the molecules collection.")
    .def_method("add_protein",&Molecules::add_protein)
    .def_doc("add_protein","Adds the protein to the molecules collection.")
    .def_method("add_nucleicacid",&Molecules::add_nucleicacid)
    .def_doc("add_nucleicacid","Adds the nucleicacid to the molecules collection.")
    .def_method("atom_size",&Molecules::atom_size)
    .def_doc("atom_size","Returns the number of atoms in the molecules collection.")
    .def_method("atom_name",&Molecules::atom_name)
    .def_doc("atom_name","Return the name of the atom with atom descriptor specified by the first parameter. The second sets the mode.")
    .def_method("fragment_size",&Molecules::fragment_size)
    .def_doc("fragment_size","Returns the number of fragments in the molecules collection.")
    .def_method("protein_size",&Molecules::protein_size)
    .def_doc("protein_size","Returns the number of proteins in the molecules collection.")
    .def_method("nucleicacid_size",&Molecules::nucleicacid_size)
    .def_doc("nucleicacid_size","Returns the number of nucleicacids in the molecules collection.")
    .def_method("lipid_size",&Molecules::lipid_size)
    .def_doc("lipid_size","Returns the number of lipids in the molecules collection.")
    .def_method("solvent_size",&Molecules::solvent_size)
    .def_doc("solvent_size","Returns the number of solvents in the molecules collection.")
    .def_method<Protein & (Molecules::*)(int)>("protein",&Molecules::protein)
    .def_doc("protein","Returns a reference to the i-th protein.")
    .def_method<NucleicAcid & (Molecules::*)(int)>("nucleicacid",&Molecules::nucleicacid)
    .def_doc("nucleicacid","Returns a reference to the i-th nucleicacid.")
    .def_method<Lipid & (Molecules::*)(int)>("lipid",&Molecules::lipid)
    .def_doc("lipid","Returns a reference to the i-th lipid.")
    .def_method<Solvent & (Molecules::*)(int)>("solvent",&Molecules::solvent)
    .def_doc("solvent","Returns a reference to the i-th solvent.")
    .def_method("charge",&Molecules::charge)
    .def_doc("charge","Returns the charge of the i-th atom.")
    .def_method("cog",&Molecules::cog)
    .def_method("com",&Molecules::com)
    .def_method("translate",&Molecules::translate)
    .def_method("sequence",&Molecules::sequence)
    .def_method("find_atom",&Molecules::find_atom)
    .def_method<int  (Molecules::*)(int) const>("find_fragment",&Molecules::find_fragment)
    .def_method("distance",&Molecules::distance)
    //    .def_method("distanceP2",&Molecules::distanceP2)
    ;
    

  Class<XMolecules,ArgList2<const Molecules &,int> >(module,"xmolecules")
    .def_method("set",&XMolecules::set)
    .def_method("clear",&XMolecules::clear)
    .def_method("size",&XMolecules::size);

  ZObj::BinaryOperator<INDEX>::add_ifndef<XMolecules,int>
      (index_operator<XMolecules,int,Molecules &>);








  Class<DSSPSecStruct,ArgList1<const Molecules &> >(module,"dssp")
    .def_method("secstruct",&DSSPSecStruct::secstruc_hbonds)
    ;

  Class<AlmSecStruct,ArgList1<const Molecules &> >(module,"almsecstruct")
    .def_method("secstruct",&AlmSecStruct::secstruct)
    .def_method("secstructf",&AlmSecStruct::secstructf)
    ;

  Class<MakeBetaPair>(module,"makebeta_pair")
    .def_method("apply",&MakeBetaPair::apply)
    .def_method("parallel",&MakeBetaPair::parallel)
    .def_method("up",&MakeBetaPair::up)
    .def_method("betabb",&MakeBetaPair::betabb)
    .def_method("offset",&MakeBetaPair::offset)
    .def_method("set_parallel",&MakeBetaPair::set_parallel)
    .def_method("set_up",&MakeBetaPair::set_up)
    .def_method("set_betabb",&MakeBetaPair::set_betabb)
    .def_method("set_offset",&MakeBetaPair::set_offset)
    ;
}

int f____(){ return ZType<float>::type_id_+ZType<float&>::type_id_;}
extern "C" {
  void init_molecules(){
    typedef double (*d_fun_)(const Molecules & mols, const NameSelection & s);
    Module module("molecules","Build and modify molecules");

    module.def_function("mols_coor","Return a Coor structure with the coordinates of the molecule.",mols_coor);  
    init_molecules_types(module.self());

    module.def_function("protein2alm","Prints the protein to the file specified  in almost format.",prot2alm)
      .def_function("protein2pdb","Prints the protein to the file specified  in PDB format.",prot2pdb)
      .def_function("solvent2alm","Prints the solvent to the file specified  in almost format.",solv2alm)
      .def_function("solvent2pdb","Prints the solvent to the file specified  in PDB format.",solv2pdb)
      .def_function("lipid2alm","Prints the lipid to the file specified  in almost format.",lipid2alm)
      .def_function("lipid2pdb","Prints the lipid to the file specified  in PDB format.",lipid2pdb)
      .def_function("nucleicacid2alm","Prints the na to the file specified  in almost format.",na2alm)
      .def_function("nucleicacid2pdb","Prints the na to the file specified  in PDB format.",na2pdb)
      .def_function("molecules2pdb","Prints the molecules to the file specified  in PDB format.",mol2pdb)
      .def_function("boxsolvate","Creates a water box around a protein.",boxsolvate)
      .def_function("spheresolvate","Creates a water sphere around a protein.",spheresolvate)
      .def_function("molspheresolvate","Creates a water sphere around molecules.",molspheresolvate)
      .def_function("shellsolvate","Creates a water shell around a protein.",shellsolvate)
      .def_function("molshellsolvate","Creates a water shell around molecules.",molshellsolvate)
      .def_function("fit","Fits second molecules to first using selected atoms.",fit)
      .def_function("subfit","Fits second molecules to the first one using selection 1 and returns the rmsd computed on the atoms of selection 2.",subfit)
      .def_function("fit2","Fits second molecules to first using the two sets of selected atoms.",fit2)
      .def_function("mfit","Fits second molecules to fitrs using selected atoms and weighting each contribution  with the atom mass.",mfit)
      .def_function("set_torsion","Sets torsion between desired atom to the specified value.",set_torsion)
      .def_function("set_torsion2","Sets torsion between desired atom to the specified value.",set_torsion2)
      .def_function("randomize_backbone","Randomizes backbone torsion angles.",randomize_backbone)
      .def_function("coordinate","Returns coordinates of the molecules",mol_coor)
      .def_function("ca_secstruct","Computes secondary structure using CA-atoms",secstruc_ca)
      .def_function("dssp_secstruct","Computes secondary structure using CA-atoms",secstruc_hbonds)
      .def_function("check_bonds","Checks bond lengths",check_bonds)
      .def_function("rmsd_prof","Compute amino-acid rmsd profile",rmsd_prof)
      .def_function("mutate","Wrapper",mutate)
      .def_function("ramarmsd","Computes mean distance on Ramachandram plot",ramarmsd)
      .def_function<d_fun_>("rgyr","Computes radious of gyration of the selected atoms.",Almost::rgyr)
      .def_function("dot","Computed dot product",coor3D_scal)
      .def_function("auto_connect","Find bonded structire of molecules",auto_connect)
      .def_function("auto_connect_pdb","Find bonded structire of molecules",auto_connect_pdb)
      .def_function("maxsub_const","",Almost::maxsub_const)
      .def_function("render_tachyon","",Almost::render_tachyon)
      .def_function("find_clashes","",Almost::find_clashes)
      .def_function("psf","Writes charmm PSF (protein structure file).",mol2psf)
      .def_function("find_disu","Cpmputes the number of disulfide bonds with ditstance lower than the cutoff",Almost::find_disu)
      ;

    module
      .def_const("LONG",(int)Molecules::LONG)
      .def_const("SHORT",(int)Molecules::SHORT)
      .def_const("BASE",(int)Molecules::BASE)
      ;

  }
}