/***************************************************************************
                          pdb.cpp  -  description
                             -------------------
    begin                : Mon Feb 11 14:09:48 2002
    copyright            : (C) 2002 by Cavalli Andrea
    email                : cavalli.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 <strstream>
#include <pdb.h>
#include <boost/iostreams/filtering_streambuf.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/copy.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#ifdef HAVE_BZLIB_H
#include <boost/iostreams/filter/bzip2.hpp>
#endif
#include<boost/regex.hpp>
using namespace Almost;
void rename__frag(PDBFragment & frag, string a1, string a2){
  for(size_t i=0;i<frag.size();i++){
    if(frag[i].name== a1) frag[i].name = a2;
  }
}

void rename__seg(PDBSegment & seg,string r1, string a1, string a2){
  bool num = false;
  if(r1.find("_")!=string::npos)
    num = 1;
  for(size_t i=0;i<seg.size();i++){
    if(!num){
      if((seg[i].name()==r1)||(r1=="*")){
	rename__frag(seg[i],a1,a2);
      }
    }
    else {
      if((seg[i].name(PDBFragment::LONG)==r1)||(r1=="*")){
	rename__frag(seg[i],a1,a2);
      }
    }
  }
}

void rename__mod(PDBModel & mod,string r1, string a1, string a2){
  for(size_t i=0;i<mod.size();i++)
    rename__seg(mod[i],r1,a1,a2);
}

void PDB::rename_atom(string r1, string a1, string a2){

  for(size_t i=0;i<model.size();i++)
    rename__mod(*model[i],r1,a1,a2);

}



void PDB::rename_frag(string r1, string r2){
  for(size_t i=0;i<model.size();i++){
    PDBModel & mod = *model[i];
    for(size_t j=0;j<mod.size();j++){
      PDBSegment & seg = mod[j];
      for(size_t k=0;k<seg.size();k++){
	if(seg[k].name()==r1) seg[k].set_name(r2);
      }
    }
  }
}


namespace PDBUtils {

  string simplifyWS(string str){
    size_t pos;
    while((pos= str.find("  "))!=string::npos){
      str.replace(str.begin()+pos,
		  str.begin()+pos+2," ");
    }
    return str;
  }


  vector<string> split(string s,char c=' '){
    string t;
    vector<string> r;
    size_t i=0;
    while(i<s.size()){
      if(s[i]!=c){ t += s[i]; ++i;}
      else { if(t.size()) r.push_back(t);t="";i++;}
    }
    if(t.size()) r.push_back(t);
    return r;
  }
}


PDB::PDB(string file){
  keep_hetatm_ = false;
  currentPDBModel = NULL;
  currentPDBSegment = NULL;
  segName = "";
  chainID = "";
  multiPDBModel = false;
  
  int comp = 0;
  int len = file.size()-1;
  if(file[len]=='z'&&file[len-1]=='g')
    comp =1;
  if(file[len]=='2'&&file[len-1]=='z'&&file[len-2]=='b')
    comp=2;
  if(comp==0){
    ifstream in;
    in.open(file.c_str());
    if(!in){
      cout<<"FATAL ERROR: no such file: "<<file<<endl;
      exit(-1);
    }
    read(in);
    in.close();
  }
  if(comp==1){
    ifstream in;
    in.open(file.c_str(),ios_base::in | ios_base::binary);
    boost::iostreams::filtering_istream in_gz;
    in_gz.push(boost::iostreams::gzip_decompressor());
    in_gz.push(in);
    read(in_gz);
    in.close();
  }
#ifdef HAVE_BZLIB_H
  if(comp==2){
    ifstream in;
    in.open(file.c_str(),ios_base::in | ios_base::binary);
    boost::iostreams::filtering_istream in_bz2;
    in_bz2.push(boost::iostreams::bzip2_decompressor());
    in_bz2.push(in);
    read(in_bz2);
    in.close();
  }
#endif
  //test if all models are consistent.
  //Check 1
  bool nsegm_check = true;
  //  bool nfrag_check = true;
  //all models contain the same number of segments ....
  if(model.size()>1){
    size_t nsegm = model[0]->size(); //insegments
    for(size_t i=1;i<model[i]->size();++i){
      if(model[i]->size()!= nsegm ){ nsegm_check=false; break;}
    }
    if(!nsegm_check){
      cout<<"PDB WARNING: the number of segment is not the same in all MODELS!"<<endl;
    }
  }
}

PDB::PDB(string file, bool keep){
  keep_hetatm_ = keep;
  currentPDBModel = NULL;
  currentPDBSegment = NULL;
  segName = "";
  chainID = "";
  multiPDBModel = false;
  
  int comp = 0;
  int len = file.size()-1;
  if(file[len]=='z'&&file[len-1]=='g')
    comp =1;
  if(file[len]=='2'&&file[len-1]=='z'&&file[len-2]=='b')
    comp=2;
  if(comp==0){
    ifstream in;
    in.open(file.c_str());
    if(!in){
      cout<<"FATAL ERROR: no such file: "<<file<<endl;
      exit(-1);
    }
    read(in);
    in.close();
  }
  if(comp==1){
    ifstream in;
    in.open(file.c_str(),ios_base::in | ios_base::binary);
    boost::iostreams::filtering_istream in_gz;
    in_gz.push(boost::iostreams::gzip_decompressor());
    in_gz.push(in);
    read(in_gz);
    in.close();
  }
#ifdef HAVE_BZLIB_H
  if(comp==2){
    ifstream in;
    in.open(file.c_str(),ios_base::in | ios_base::binary);
    boost::iostreams::filtering_istream in_bz2;
    in_bz2.push(boost::iostreams::bzip2_decompressor());
    in_bz2.push(in);
    read(in_bz2);
    in.close();
  }
#endif
  //test if all models are consistent.
  //Check 1
  bool nsegm_check = true;
  //  bool nfrag_check = true;
  //all models contain the same number of segments ....
  if(model.size()>1){
    size_t nsegm = model[0]->size(); //insegments
    for(size_t i=1;i<model[i]->size();++i){
      if(model[i]->size()!= nsegm ){ nsegm_check=false; break;}
    }
    if(!nsegm_check){
      cout<<"PDB WARNING: the number of segment is not the same in all MODELS!"<<endl;
    }
  }
}
      
PDB::~PDB(){
  for(size_t i=0;i<model.size();i++) delete model[i];
  model.clear();
}

//copy constructor
PDB::PDB(const PDB & pdb){
  multiPDBModel = pdb.multiPDBModel;
  for(size_t i=0;i<pdb.model.size();i++){
    model.push_back(new PDBModel(pdb[i]));
  }
  header_ = pdb.header_;
  title_ = pdb.title_;
  compound_ = pdb.compound_;
  source_ = pdb.source_;
  keywords_ = pdb.keywords_;
  expdata_ = pdb.expdata_;
  author_ = pdb.author_;
  journal_ = pdb.journal_;
  remark_ = pdb.remark_;
  seqres_ = pdb.seqres_;
}

string PDB::info() const {
  string str;
  strstream out;

  out<<model.size();
  string tmp;
  out>>tmp;
  str =  "PDB file contains "+tmp+" model(s).\n";
  for(size_t i=0;i<model.size();i++){
    {
      strstream out2;
      out2<<i;
      string tmp2;
      out2>>tmp2;
      str += "PDB model "+ tmp2 + " contains ";
    }
    {
      strstream out2;
      out2<< model[i]->size();
      string tmp2;
      out2>>tmp2;
      str += tmp2 + " segment(s)\n";
      for(size_t s = 0; s<model[i]->size();s++){
	str += (*model[i])[s].chainID()+" ";
      }
      str +="\n";
    }
  }
  return str;
}

int PDB::align_seqres(string c, int m, int i, string &r){
  using namespace boost;
  if(m>=(*this).size()) return -1;
  if(i>=(*this)[m].size()) return -1;
  PDBSegment & seg = (*this)[m][i];
  if(seg.size()==0) return -1;

  r = three2one(seg[0].name());
  int fseq = seg[0].fragSeq();
  for(int i=1;i<seg.size();i++){
    if(!(seg[i].fragSeq()>fseq)) return -1;
    for(int k=fseq+1;k<seg[i].fragSeq();k++){ r +=".";}
    fseq = seg[i].fragSeq();
    r += three2one(seg[i].name());
  }
  //  cout<<"DEBUG "<<r<<endl;
  //  vector<string> seq1 = (*this)[m][i].sequence(PDBSegment::SHORT);
  vector<string> seq_res = seqres(c);
  string seq;
  for(size_t i=0;i<seq_res.size();i++){
    seq += three2one(seq_res[i]);
  }
  match_results<std::string::const_iterator> what;
  match_flag_type flags = match_default;
  regex exp(r);
  string::const_iterator begin = seq.begin(),end = seq.end();
  if(regex_search(begin,end,what,exp,flags)){
    return what[0].first - seq.begin();
  }
  else return -1;
}

void PDB::read(istream &in){
  PDBAtomLine * atomLine;
  //   PDBHeaderLine * headerLine;
  while(in){
    PDBLine line;
    
    line.readLine(in);
    int type = line.type();
    
    switch(type){
    case PDBType::HEADER:
      header_.push_back(PDBHeaderLine(line.line()));
      break;
    case PDBType::TITLE:
      title_.push_back(PDBTitleLine(line.line()));
      break;
    case PDBType::COMPND:
      compound_.push_back(PDBCompndLine(line.line()));
      break;
    case PDBType::SOURCE:
      source_.push_back(PDBSourceLine(line.line()));
      break;
    case PDBType::KEYWDS:
      keywords_.push_back(PDBKeywordsLine(line.line()));
      break;
    case PDBType::EXPDTA:
      expdata_.push_back(PDBExpdataLine(line.line()));
      break;
    case PDBType::AUTHOR:
      author_.push_back(PDBAuthorLine(line.line()));
      break;
    case PDBType::JRNL:
      journal_.push_back(PDBJournalLine(line.line()));
      break;
    case PDBType::REMARK:
      remark_.push_back(PDBRemarkLine(line.line()));
      break;
    case PDBType::SEQRES:
      seqres_.push_back(PDBSeqResLine(line.line()));
      break;
    case PDBType::MODEL:
      currentPDBModel = new PDBModel();
      break;
    case PDBType::ENDMDL:
      model.push_back(currentPDBModel);
      segName="";
      chainID = "";
      currentPDBModel =NULL;
      currentPDBSegment = NULL;
      break;
    case PDBType::TER:
      //      currentPDBModel->segment.push_back(currentPDBSegment);
      currentPDBSegment = NULL;
      segName="";
      chainID = "";
      break;
    case PDBType::END:
      if(model.size()==0){
	if(currentPDBModel->size()>0)
	  model.push_back(currentPDBModel);
      }
      break;
    case PDBType::ATOM:
    case PDBType::HETATM:
      if(type==PDBType::HETATM&&!keep_hetatm_){
	//	cout<<"DEBUG discarding HETATM"<<endl;
	break;
      }
      if(type==PDBType::ATOM)
	atomLine = new PDBAtomLine(line.line());
      else
	atomLine = new PDBAtomLine(line.line(),PDBAtomLine::HETATM);
      //CHARMM ...
      //       if((atomLine->altLoc !=' ')&&(atomLine->altLoc !='A')){
      // 	cout <<"Alt location.....ignored"<<endl;
      // 	delete atomLine;
      // 	break;
      //       }
      //       cout<<"DEBUG "<<atomLine->serial<<" "
      // 	  <<atomLine->name<<" "
      // 	  <<atomLine->resName<<"\n";
      if(currentPDBModel == NULL)
	currentPDBModel = new PDBModel();
      if(currentPDBSegment == NULL){
	currentPDBSegment = new PDBSegment();
	currentPDBModel->segment.push_back(currentPDBSegment);
	segName = atomLine->segment();
	chainID = atomLine->chainID;
	currentPDBSegment->set_name(segName);
	currentPDBSegment->set_chainID(chainID);

	//	cout<<"Segment "<<segName<<endl;

	currentPDBFragment = new PDBFragment();
	currentPDBSegment->fragment.push_back(currentPDBFragment);
	fragName = atomLine->resName;
	fragSeq = atomLine->resSeq;

	//	cout<<"Framgment "<<fragName<<endl;


	currentPDBFragment->set_name(fragName);
	currentPDBFragment->set_fragSeq(fragSeq);
	currentPDBFragment->set_iCode(atomLine->iCode);
      }
      else if((segName!=atomLine->segment())){
	//if RESNAME and REDNUM unchanged ... ignore
	string old_res_name = currentPDBSegment->fragment.back()->name();
	string current_res_name = atomLine->resName;
	if((old_res_name!=current_res_name)&&
	   (currentPDBSegment->fragment.back()->fragSeq()
	    >atomLine->resSeq)){
	  currentPDBSegment = new PDBSegment();
	  currentPDBModel->segment.push_back(currentPDBSegment);
	  segName = atomLine->segment();
	  chainID = atomLine->chainID;
	  currentPDBSegment->set_name(segName);
	  currentPDBSegment->set_chainID(chainID);
	  
	  //	  cout<<"Segment "<<segName<<endl;

	  currentPDBFragment = new PDBFragment();
	  currentPDBSegment->fragment.push_back(currentPDBFragment);
	  fragName = atomLine->resName;
	  fragSeq = atomLine->resSeq;

	  //	  cout<<"Framgment "<<fragName<<endl;

	  currentPDBFragment->set_name(fragName);
	  currentPDBFragment->set_fragSeq(fragSeq);	
	  currentPDBFragment->set_iCode(atomLine->iCode);  
	}
      }
      
      if((fragName != atomLine->resName)||(fragSeq != atomLine->resSeq)){
	currentPDBFragment = new PDBFragment();
	currentPDBSegment->fragment.push_back(currentPDBFragment);
	fragName = atomLine->resName;
	
	//	cout<<"Framgment "<<fragName<<endl;


	fragSeq = atomLine->resSeq;
	currentPDBFragment->set_name(fragName);
	currentPDBFragment->set_fragSeq(fragSeq);
	currentPDBFragment->set_iCode(atomLine->iCode);
      }
      currentPDBFragment->atom.push_back(atomLine);
      //cout<<atomLine->name<<endl;

      break;
    default:;
    }
  }

  if(model.size()>0) multiPDBModel = true;
  if(model.size()==0){
    if(currentPDBModel->size()>0)
      model.push_back(currentPDBModel);
    else{
      cout <<"no model found"<<endl;
    }
  }
}


PDBModel & PDB::operator[](int i){
  return *model[i];
}

const PDBModel & PDB::operator[](int i) const {
  return *model[i];
}

struct resLess {
  bool operator()(string a, string b){
    string a1,a2;
    for(size_t i=0;i<a.size();i++){
      if(!isdigit(a[i])) continue;
      a1 +=a[i];
    }
    for(size_t i=0;i<b.size();i++){
      if(!isdigit(b[i])) continue;
      a2 +=b[i];
    }
    
    return (atoi(a1.c_str())< atoi(a2.c_str()));
  }
};

//vector<string> PDBSegment::sequence(int TYPE=LONG) const {



//content access

string PDB::header() const {
  string s;
  for(size_t i=0;i<header_.size();i++){
    s += header_[i].classification+" "+header_[i].date+" "
      +header_[i].idCode;
  }
  return PDBUtils::simplifyWS(s);
}

string PDB::title() const {
  string s;
  for(size_t i=0;i<title_.size();i++){
    s += title_[i].text;
  }
  return PDBUtils::simplifyWS(s);
}

string PDB::compound() const {
  string s;
  for(size_t i=0;i<compound_.size();i++){
    s += compound_[i].text;
  }
  return PDBUtils::simplifyWS(s);
}

string PDB::source() const {
  string s;
  for(size_t i=0;i<source_.size();i++){
    s += source_[i].text;
  }
  return PDBUtils::simplifyWS(s);
}

string PDB::author() const {
  string s;
  for(size_t i=0;i<author_.size();i++){
    vector<string> a = PDBUtils::split(author_[i].text,',');
    for(size_t j=0;j<a.size();j++)  s += a[j]+" ";
  }
  return PDBUtils::simplifyWS(s); 
}

string PDB::journal(string key) const {
  using namespace boost;
  string s;
  string res;  
  string rs = "(^"+key+")\\s+(\\d*\\s+)*(.*)";
  regex exp(rs);

  for(size_t i=0;i<journal_.size();i++){
    smatch what;
    s = PDBUtils::simplifyWS(journal_[i].text);
    if(regex_match(s,what,exp)){
      //      cout<<"REXM "<<what[0]<<endl;
      //      cout<<"REXM "<<what[1]<<endl;
      //      cout<<"REXM "<<what[2]<<endl;
      //      cout<<"REXM "<<what[3]<<endl;
      res += what[3];
    }
  }
  
  return res;

}

vector<string> PDB::seqres(string c) const {
  string s;
  if(seqres_.size()==0) return vector<string>();
  
  if(c==""||c==" "){
    for(size_t i=0;i<seqres_.size();i++){
      s += seqres_[i].text;
      if((i+1)!=seqres_.size()) s+=" ";
    }
    s = PDBUtils::simplifyWS(s);
    return PDBUtils::split(s);
  }

  for(size_t i=0;i<seqres_.size();i++){
    if(seqres_[i].chain==c[0]){
      s += seqres_[i].text;
      if((i+1)!=seqres_.size()) s+=" ";
    }
  }
  s = PDBUtils::simplifyWS(s);
  return PDBUtils::split(s);    
}


string PDB::remark(int i) const {
  string s;
  if(remark_.size())
    s = remark_[0].text;
  for(size_t i=1;i<remark_.size();i++){
    s += "\n"+remark_[i].text;
  }
  return s;
}

vector<string> PDB::remarks() const {
  vector<string> res;
  for(size_t i=1;i<remark_.size();i++){
    res.push_back(remark_[i].text);
  }
  return res;
}

void PDB::normalizeNames(string type){
  //roules ----
  
  if(type=="CHARMM22"){
    normalizeCh22();
  }
  else if(type=="CHARMM19"){
    normalizeCh19();
  }
  else{
    cout <<"Name normalization rule '"<<type<<"' unknown."<<endl;
  }
}

void PDB::normalizeCh19(){
  //     for(int i=0;i<(*this).size();++i){
  //     for(int j=0;j<(*this)[i].size();++j){
  //       bool patched = false;
  //       for(int k=0;k<(*this)[i][j].size();++k){
  // 	//loop on PDBAtomLine
  // 	string name = (*this)[i][j][k].name;
	
	
  // 	if(name =="H1") (*this)[i][j][k].name = "HT1";
  // 	else if(name =="H2") (*this)[i][j][k].name = "HT2";
  // 	else if(name =="H3") (*this)[i][j][k].name = "HT3";
  // 	else if(name =="OXT"){
  // 	  (*this)[i][j][k].name = "OT2";
  // 	  patched = true;
  // 	}

  // 	if((*this)[i][j][k].resName =="ILE"){
  // 	  if(name =="CD1") (*this)[i][j][k].name = "CD";
  // 	  else if(name =="HD11") (*this)[i][j][k].name = "HD1";
  // 	  else if(name =="HD12") (*this)[i][j][k].name = "HD2";
  // 	  else if(name =="HD13") (*this)[i][j][k].name = "HD3";
  // 	}
  //       }
      
  //       //change names in terminal ........
  //       if(patched){
  // 	int k = (*this)[i][j].size()-1;
  // 	string resName = (*this)[i][j][k].resName;
  // 	int resSeq = (*this)[i][j][k].resSeq;

  // 	while(k>=0){
  // 	  if((*this)[i][j][k].resSeq!= resSeq) break;
  // 	  else if((*this)[i][j][k].name =="O") {
  // 	    (*this)[i][j][k].name = "OT1";
  // 	    break;
  // 	  }
  // 	  else --k;
  // 	}
  //       }
  //     }
  //   }
}


void PDB::normalizeCh22(){
 
  //   for(int i=0;i<(*this).size();++i){
  //     for(int j=0;j<(*this)[i].size();++j){
  //       bool patched = false;
  //       for(int k=0;k<(*this)[i][j].size();++k){
  // 	//loop on PDBAtomLine
  // 	string name = (*this)[i][j][k].name;
	
  // 	if(name == "H") (*this)[i][j][k].name = "HN";
  // 	else if(name =="H1") (*this)[i][j][k].name = "HT1";
  // 	else if(name =="H2") (*this)[i][j][k].name = "HT2";
  // 	else if(name =="H3") (*this)[i][j][k].name = "HT3";
  // 	else if(name =="OXT"){
  // 	  (*this)[i][j][k].name = "OT2";
  // 	  patched = true;
  // 	}

  // 	if((*this)[i][j][k].resName =="ILE"){
  // 	  if(name =="CD1") (*this)[i][j][k].name = "CD";
  // 	  else if(name =="HD11") (*this)[i][j][k].name = "HD1";
  // 	  else if(name =="HD12") (*this)[i][j][k].name = "HD2";
  // 	  else if(name =="HD13") (*this)[i][j][k].name = "HD3";
  // 	}
  //       }
      
  //       //change names in terminal ........
  //       if(patched){
  // 	int k = (*this)[i][j].size()-1;
  // 	string resName = (*this)[i][j][k].resName;
  // 	int resSeq = (*this)[i][j][k].resSeq;

  // 	while(k>=0){
  // 	  if((*this)[i][j][k].resSeq!= resSeq) break;
  // 	  else if((*this)[i][j][k].name =="O") {
  // 	    (*this)[i][j][k].name = "OT1";
  // 	    break;
  // 	  }
  // 	  else --k;
  // 	}
  //       }
  //     }
  //   }

}
	  
bool comp(string a1,string a2){
  for(size_t i=0;i<a1.size();i++)
    a1[i] = tolower(a1[i]);
  if(a1==a2) return true;
  return false;
}

void info_ala(){
  aout<<"Aminoacid: Alanine ALA A\n";
  aout<<"     | \n";
  aout<<"  HN-N \n";
  aout<<"     |     HB1 \n";
  aout<<"     |    / \n";
  aout<<"  HA-CA--CB-HB2 \n";
  aout<<"     |    \\ \n";
  aout<<"     |     HB3 \n";
  aout<<"   O=C \n";
  aout<<"     |     \n";
}

void info_arg(){
  aout<<"Aminoacid: Arginine ARG R \n";
  aout<<"     |                      HH11 \n";
  aout<<"  HN-N                       | \n";
  aout<<"     |   HB1 HG1 HD1 HE     NH1-HH12 \n";
  aout<<"     |   |   |   |   |    //(+)   \n";
  aout<<"  HA-CA--CB--CG--CD--NE--CZ \n";
  aout<<"     |   |   |   |         \\ \n";
  aout<<"     |   HB2 HG2 HD2        NH2-HH22 \n";
  aout<<"   O=C                       | \n";
  aout<<"     |                      HH21 \n";
}

void info_asn(){
  aout<<"Aminoacid: Asparagine ASN N \n";
  aout<<"     |        \n";
  aout<<"  HN-N        \n";
  aout<<"     |   HB1 OD1    HD21 (cis to OD1) \n";
  aout<<"     |   |   ||    / \n";
  aout<<"  HA-CA--CB--CG--ND2 \n";
  aout<<"     |   |         \\ \n";
  aout<<"     |   HB2        HD22 (trans to OD1) \n";
  aout<<"   O=C            \n";
  aout<<"     |            \n";
}

void info_asp(){
  aout<<"Aminoacid: Aspartic acid ASP D \n";
  aout<<"     |        \n";
  aout<<"  HN-N        \n";
  aout<<"     |   HB1   OD1 \n";
  aout<<"     |   |    // \n";
  aout<<"  HA-CA--CB--CG \n";
  aout<<"     |   |    \\ \n";
  aout<<"     |   HB2   OD2(-) \n";
  aout<<"   O=C            \n";
  aout<<"     |    \n";
}

void info_cys(){
  aout<<"Aminoacid: Cysteine CYS C \n";
  aout<<"     |        \n";
  aout<<"  HN-N        \n";
  aout<<"     |   HB1 \n";
  aout<<"     |   |   \n";
  aout<<"  HA-CA--CB--SG \n";
  aout<<"     |   |     \\ \n";
  aout<<"     |   HB2    HG1 \n";
  aout<<"   O=C            \n";
  aout<<"     |            \n";
}

void info_gln(){
  aout<<"Aminoacid: Glutamine GLN Q \n";
  aout<<"     |           \n";
  aout<<"  HN-N           \n";
  aout<<"     |   HB1 HG1 OE1   HE21 (cis to OE1) \n";
  aout<<"     |   |   |   ||    / \n";
  aout<<"  HA-CA--CB--CG--CD--NE2 \n";
  aout<<"     |   |   |         \\ \n";
  aout<<"     |   HB2 HG2       HE22 (trans to OE1) \n";
  aout<<"   O=C           \n";
}

void info_glu(){
  aout<<"Aminoacid: Glutamic acid GLU E \n";
  aout<<"     |           \n";
  aout<<"  HN-N           \n";
  aout<<"     |   HB1 HG1   OE1 \n";
  aout<<"     |   |   |    // \n";
  aout<<"  HA-CA--CB--CG--CD \n";
  aout<<"     |   |   |    \\ \n";
  aout<<"     |   HB2 HG2   OE2(-) \n";
  aout<<"   O=C           \n";
  aout<<"     |           \n";
}

void info_gly(){
  aout<<"Aminoacid: Glycine GLY G \n";
  aout<<"     |           \n";
  aout<<"  HN-N           \n";
  aout<<"     |   HB1 HG1   OE1 \n";
  aout<<"     |   |   |    // \n";
  aout<<"  HA-CA--CB--CG--CD \n";
  aout<<"     |   |   |    \\ \n";
  aout<<"     |   HB2 HG2   OE2(-) \n";
  aout<<"   O=C           \n";
  aout<<"     |          \n";
}

void info_his(){
  aout<<"Aminoacid: Histidine HIS (HSD HSE HSP) H \n";
  aout<<"Neutral HIS, proton on ND1\n";
  aout<<"     |          HD1    HE1\n";
  aout<<"  HN-N           |     /\n";
  aout<<"     |   HB1    ND1--CE1\n";
  aout<<"     |   |     /      ||\n";
  aout<<"  HA-CA--CB--CG       ||\n";
  aout<<"     |   |     \\\\     ||\n";
  aout<<"     |   HB2    CD2--NE2\n";
  aout<<"   O=C           |\n";
  aout<<"     |          HD2\n";
  aout<<"Neutral His, proton on NE2\n";
  aout<<"     |                 HE1\n";
  aout<<"  HN-N             __  /\n";
  aout<<"     |   HB1    ND1--CE1\n";
  aout<<"     |   |     /      |\n";
  aout<<"  HA-CA--CB--CG       |\n";
  aout<<"     |   |     \\\\     |\n";
  aout<<"     |   HB2    CD2--NE2\n";
  aout<<"   O=C           |     \\\n";
  aout<<"     |          HD2    HE2\n";
  aout<<"Protonated His\n";
  aout<<"     |          HD1    HE1\n";
  aout<<"  HN-N           |     /\n";
  aout<<"     |   HB1    ND1--CE1\n";
  aout<<"     |   |     /      ||\n";
  aout<<"  HA-CA--CB--CG       ||\n";
  aout<<"     |   |     \\\\     ||\n";
  aout<<"     |   HB2    CD2--NE2(+)\n";
  aout<<"   O=C           |     \\\n";
  aout<<"     |          HD2    HE2\n";
}

void info_ile(){
  aout<<"Aminoacid: Isoleucine ILE I \n";
  aout<<"     |    HG21 HG22\n";
  aout<<"  HN-N      | / \n";
  aout<<"     |     CG2--HG23\n";
  aout<<"     |    /\n";
  aout<<"  HA-CA--CB-HB    HD1\n";
  aout<<"     |    \\\\       /\n";
  aout<<"     |     CG1--CD--HD2\n";
  aout<<"   O=C    / \\\\     \\\\	 \n";
  aout<<"     | HG11 HG12  HD3\n";
}

void info_leu(){
  aout<<"Aminoacid: Leucine LEU L \n";
  aout<<"     |        HD11 HD12\n";
  aout<<"  HN-N          | /\n";
  aout<<"     |   HB1   CD1--HD13\n";
  aout<<"     |   |    /\n";
  aout<<"  HA-CA--CB--CG-HG\n";
  aout<<"     |   |    \\ \n";
  aout<<"     |   HB2   CD2--HD23\n";
  aout<<"   O=C          | \\\n";
  aout<<"     |        HD21 HD22\n";
}


void info_lys(){
  aout<<"Aminoacid: Lysine LYS K \n";
  aout<<"     |                   \n";
  aout<<"  HN-N                   \n";
  aout<<"     |   HB1 HG1 HD1 HE1    HZ1\n";
  aout<<"     |   |   |   |   |     /   \n";
  aout<<"  HA-CA--CB--CG--CD--CE--NZ--HZ2\n";
  aout<<"     |   |   |   |   |     \\\n";
  aout<<"     |   HB2 HG2 HD2 HE2    HZ3\n";
  aout<<"   O=C                  \n";
  aout<<"     |                  \n";
}

void info_met(){
  aout<<"Aminoacid: Methionine MET M \n";
  aout<<"     |                   \n";
  aout<<"  HN-N                   \n";
  aout<<"     |   HB1 HG1     HE1 \n";
  aout<<"     |   |   |       |   \n";
  aout<<"  HA-CA--CB--CG--SD--CE--HE3\n";
  aout<<"     |   |   |       |   \n";
  aout<<"     |   HB2 HG2     HE2 \n";
  aout<<"   O=C                  \n";
  aout<<"     |                  \n";
}

void info_phe(){
  aout<<"Aminoacid: Phenylalanine PHE F \n";
  aout<<"     |        HD1  HE1    \n";
  aout<<"  HN-N         |    |   \n";
  aout<<"     |   HB1  CD1--CE1\n";
  aout<<"     |   |    //     \\\\\n";
  aout<<"  HA-CA--CB--CG      CZ--HZ\n";
  aout<<"     |   |    \\  __  /\n";
  aout<<"     |   HB2  CD2--CE2\n";
  aout<<"   O=C         |    |  \n";
  aout<<"     |        HD2  HE2   \n";
}

void info_pro(){
  aout<<"Aminoacid: Proline PRO P \n";
  aout<<"       HD1 HD2\n";
  aout<<"     |   \\ /\n";
  aout<<"     N---CD   HG1  ATOM CA   CP1     0.02\n";
  aout<<"     |     \\  /\n";
  aout<<"     |      CG\n";
  aout<<"     |     /  \\\n";
  aout<<"  HA-CA--CB   HG2\n";
  aout<<"     |   / \\\n";
  aout<<"     | HB1 HB2\n";
  aout<<"   O=C\n";
  aout<<"     |\n";
}

void info_ser(){
  aout<<"Aminoacid: Serine SER S \n";
  aout<<"     |       \n";
  aout<<"  HN-N       \n";
  aout<<"     |   HB1\n";
  aout<<"     |   |  \n";
  aout<<"  HA-CA--CB--OG\n";
  aout<<"     |   |     \\\n";
  aout<<"     |   HB2    HG1\n";
  aout<<"   O=C           \n";
  aout<<"     |           \n";
}

void info_thr(){
  aout<<"Aminoacid: Threonine THR T \n";
  aout<<"     |  \n";
  aout<<"  HN-N  \n";
  aout<<"     |     OG1--HG1\n";
  aout<<"     |    /\n";
  aout<<"  HA-CA--CB-HB  \n";
  aout<<"     |    \\     \n";
  aout<<"     |     CG2--HG21\n";
  aout<<"   O=C    / \\    \n";
  aout<<"     | HG21 HG22 \n";
}

void info_trp(){
  aout<<"Aminoacid: Tryptophan TRP W\n";
  aout<<"     |                  HE3\n";
  aout<<"  HN-N                   |\n";
  aout<<"     |   HB1            CE3\n";
  aout<<"     |   |             /  \\\\\n";
  aout<<"  HA-CA--CB---CG-----CD2   CZ3-HZ3\n";
  aout<<"     |   |    ||     ||     |\n";
  aout<<"     |   HB2  CD1    CE2   CH2-HH2\n";
  aout<<"   O=C       /   \\   / \\  //\n";
  aout<<"     |     HD1    NE1   CZ2\n";
  aout<<"                   |     |\n";
  aout<<"                  HE1   HZ2\n";
}


void info_tyr(){
  aout<<"Aminoacid: Tyrosine TYR Y\n";
  aout<<"     |        HD1  HE1    \n";
  aout<<"  HN-N         |    |   \n";
  aout<<"     |   HB1  CD1--CE1\n";
  aout<<"     |   |   //      \\\\\n";
  aout<<"  HA-CA--CB--CG      CZ--OH\n";
  aout<<"     |   |    \\  __  /     \\\n";
  aout<<"     |   HB2  CD2--CE2     HH\n";
  aout<<"   O=C         |    |  \n";
  aout<<"     |        HD2  HE2   \n";
}

void info_val(){
  aout<<"Aminoacid: Valine VAL V\n";
  aout<<"     |    HG11 HG12\n";
  aout<< "  HN-N      | / \n";
  aout<<"     |     CG1--HG13\n";
  aout<<"     |    /\n";
  aout<<"  HA-CA--CB-HB  \n";
  aout<<"     |    \\     \n";
  aout<<"     |     CG2--HG21\n";
  aout<<"   O=C    / \\   \n";
  aout<<"     | HG21 HG22\n";
}

void info_names(){
  aout<<"RESIDUE                     ABBREVIATION                SYNONYM                 \n";
  aout<<"-----------------------------------------------------------------------------\n";
  aout<<"Alanine                     ALA                         A                       \n";
  aout<<"Arginine                    ARG                         R                       \n";
  aout<<"Asparagine                  ASN                         N                       \n";
  aout<<"Aspartic acid               ASP                         D                       \n";
  aout<<"ASP/ASN ambiguous           ASX                         B                       \n";
  aout<<"Cysteine                    CYS                         C                       \n";
  aout<<"Glutamine                   GLN                         Q                       \n";
  aout<<"Glutamic acid               GLU                         E                       \n";
  aout<<"GLU/GLN ambiguous           GLX                         Z                       \n";
  aout<<"Glycine                     GLY                         G                       \n";
  aout<<"Histidine                   HIS                         H                       \n";
  aout<<"Isoleucine                  ILE                         I                       \n";
  aout<<"Leucine                     LEU                         L                       \n";
  aout<<"Lysine                      LYS                         K                       \n";
  aout<<"Methionine                  MET                         M                       \n";
  aout<<"Phenylalanine               PHE                         F                       \n";
  aout<<"Proline                     PRO                         P                       \n";
  aout<<"Serine                      SER                         S                       \n";
  aout<<"Threonine                   THR                         T                       \n";
  aout<<"Tryptophan                  TRP                         W                       \n";
  aout<<"Tyrosine                    TYR                         Y                       \n";
  aout<<"Unknown                     UNK                                                 \n";
  aout<<"Valine                      VAL                         V                       \n";

}

void info_formulas(){
  aout<<"NAME                    CODE           FORMULA                 MOL. WT.          \n";
  aout<<"-----------------------------------------------------------------------------\n";
  aout<<"Alanine                 ALA            C3 H7 N1 O2             89.09             \n";
  aout<<"Arginine                ARG            C6 H14 N4 O2            174.20            \n";
  aout<<"Asparagine              ASN            C4 H8 N2 O3             132.12            \n";
  aout<<"Aspartic acid           ASP            C4 H7 N1 O4             133.10            \n";
  aout<<"ASP/ASN ambiguous       ASX            C4 H71/2 N11/2 O31/2    132.61            \n";
  aout<<"Cysteine                CYS            C3 H7 N1 O2 S1          121.15            \n";
  aout<<"Glutamine               GLN            C5 H10 N2 O3            146.15            \n";
  aout<<"Glutamic acid           GLU            C5 H9 N1 O4             147.13            \n";
  aout<<"GLU/GLN ambiguous       GLX            C5 H91/2 N11/2 O31/2    146.64            \n";
  aout<<"Glycine                 GLY            C2 H5 N1 O2             75.07             \n";
  aout<<"Histidine               HIS            C6 H9 N3 O2             155.16            \n";
  aout<<"Isoleucine              ILE            C6 H13 N1 O2            131.17            \n";
  aout<<"Leucine                 LEU            C6 H13 N1 O2            131.17            \n";
  aout<<"Lysine                  LYS            C6 H14 N2 O2            146.19            \n";
  aout<<"Methionine              MET            C5 H11 N1 O2 S1         149.21            \n";
  aout<<"Phenylalanine           PHE            C9 H11 N1 O2            165.19            \n";
  aout<<"Proline                 PRO            C5 H9 N1 O2             115.13            \n";
  aout<<"Serine                  SER            C3 H7 N1 O3             105.09           \n";
  aout<<"Threonine               THR            C4 H9 N1 O3             119.12           \n";
  aout<<"Tryptophan              TRP            C11 H12 N2 O2           204.23            \n";
  aout<<"Tyrosine                TYR            C9 H11 N1 O3            181.19            \n";
  aout<<"Valine                  VAL            C5 H11 N1 O2            117.15            \n";
  aout<<"Undetermined            UNK            C5 H6 N1 O3             128.16            \n";
 
}

namespace Almost {

  void info(string aa){
    if(comp(aa,"ala")){ info_ala();     aout.flush();return;}
    if(comp(aa,"arg")){ info_arg();     aout.flush();return;}
    if(comp(aa,"asn")){ info_asn();     aout.flush();return;}
    if(comp(aa,"asp")){ info_asp();     aout.flush();return;}
    if(comp(aa,"cys")){ info_cys();     aout.flush();return;}
    if(comp(aa,"gln")){ info_gln();     aout.flush();return;}
    if(comp(aa,"glu")){ info_glu();     aout.flush();return;}
    if(comp(aa,"gly")){ info_gly();     aout.flush();return;}
    if(comp(aa,"his")){ info_his();     aout.flush();return;}
    if(comp(aa,"his")){ info_his();     aout.flush();return;}
    if(comp(aa,"ile")){ info_ile();     aout.flush();return;}
    if(comp(aa,"leu")){ info_leu();     aout.flush();return;}
    if(comp(aa,"lys")){ info_lys();     aout.flush();return;}
    if(comp(aa,"met")){ info_met();     aout.flush();return;}
    if(comp(aa,"phe")){ info_phe();     aout.flush();return;}
    if(comp(aa,"pro")){ info_pro();     aout.flush();return;}
    if(comp(aa,"ser")){ info_ser();     aout.flush();return;}
    if(comp(aa,"thr")){ info_thr();     aout.flush();return;}
    if(comp(aa,"trp")){ info_trp();     aout.flush();return;}
    if(comp(aa,"tyr")){ info_tyr();     aout.flush();return;}
    if(comp(aa,"val")){ info_val();     aout.flush();return;}
    if(comp(aa,"names")){ info_names();     aout.flush();return;}
    if(comp(aa,"formulas")){ info_formulas();     aout.flush();return;}
  }

  static map<string,string> three2one_map;
  int init_three2one(){
    three2one_map["ALA"] = "A";
    three2one_map["ASX"] = "B";
    three2one_map["CYS"] = "C"; 
    three2one_map["ASP"] = "D";
    three2one_map["GLU"] = "E";
    three2one_map["PHE"] = "F";
    three2one_map["GLY"] = "G";
    three2one_map["HIS"] = "H";
    three2one_map["HSD"] = "H";
    three2one_map["ILE"] = "I";
    three2one_map["LYS"] = "K";
    three2one_map["LEU"] = "L";
    three2one_map["MET"] = "M";
    three2one_map["ASN"] = "N";
    three2one_map["PRO"] = "P";
    three2one_map["GLN"] = "Q";
    three2one_map["ARG"] = "R";
    three2one_map["SER"] = "S";
    three2one_map["THR"] = "T";
    three2one_map["VAL"] = "V";
    three2one_map["TRP"] = "W";
    three2one_map["TYR"] = "Y";
    three2one_map["GLX"] = "Z";

    three2one_map["DAL"] = "a"; //
    three2one_map["DCY"] = "c"; //
    three2one_map["DAS"] = "d";//
    three2one_map["DGL"] = "e";//
    three2one_map["DPH"] = "f";//
    three2one_map["DGY"] = "g";//
    three2one_map["DHI"] = "h";//
    three2one_map["DIL"] = "i";//
    three2one_map["DLY"] = "k";//
    three2one_map["DLE"] = "l";//
    three2one_map["DME"] = "m";//
    three2one_map["DAN"] = "n";//
    three2one_map["DPR"] = "p";//
    three2one_map["DGN"] = "q";//
    three2one_map["DAR"] = "r";//
    three2one_map["DSE"] = "s";//
    three2one_map["DTH"] = "t";//
    three2one_map["DVA"] = "v";//
    three2one_map["DTR"] = "w";//
    three2one_map["DTY"] = "y";//
    return 0;
  }

  string four2three(string damin){

    if(damin=="DALA") return "DAL";
    if(damin=="DCYS") return "DCY";
    if(damin=="DASP") return "DAS";
    if(damin=="DGLU") return "DGL";
    if(damin=="DPHE") return "DPH";
    if(damin=="DGLY") return "DGY";
    if(damin=="DHIS") return "DHI";
    if(damin=="DILE") return "DIL";
    if(damin=="DLYS") return "DLY";
    if(damin=="DLEU") return "DLE";
    if(damin=="DMET") return "DME";
    if(damin=="DASN") return "DAN";
    if(damin=="DPRO") return "DPR";
    if(damin=="DGLN") return "DGN";
    if(damin=="DARG") return "DAR";
    if(damin=="DSER") return "DSE";
    if(damin=="DTHR") return "DTH";
    if(damin=="DVAL") return "DVA";
    if(damin=="DTRP") return "DTR";
    if(damin=="DTYR") return "DTY";
    else return "XXX";
  }

  static int init = init_three2one();


  string three2one(string aa){
    if(three2one_map.size()==0)
      init = init_three2one();

    map<string,string>::iterator pos = three2one_map.find(aa);
    if(pos == three2one_map.end()){
      return "X";
    }
    else return pos->second;
  }

  string one2three(string aa){
    if(three2one_map.size()==0)
      init = init_three2one();

    map<string,string>::iterator iter = three2one_map.begin(),
      end = three2one_map.end();
    while(iter!=end){
      if(iter->second==aa) return iter->first;
      ++iter;
    }
    return "UKN";
  }

  string one2three(char aa){ string s; s+=aa; return one2three(s);}

  int one2index(string aa){
    if(aa=="A") return 0;
    if(aa=="B") return 1;
    if(aa=="C") return 2;
    if(aa=="D") return 3;
    if(aa=="E") return 4;
    if(aa=="F") return 5;
    if(aa=="G") return 6;
    if(aa=="H") return 7;
    if(aa=="I") return 8;
    if(aa=="K") return 9;
    if(aa=="L") return 10;
    if(aa=="M") return 11;
    if(aa=="N") return 12;
    if(aa=="P") return 13;
    if(aa=="Q") return 14;
    if(aa=="R") return 15;
    if(aa=="S") return 16;
    if(aa=="T") return 17;
    if(aa=="V") return 18;
    if(aa=="W") return 19;
    if(aa=="Y") return 20;
    if(aa=="Z") return 21;
  
    return 99;
  }
  

  vector<string> read_fasta(string file){
    ifstream fasta;
    fasta.open(file.c_str());
    if(!fasta){
      throw "Unable to read fasta file!";
    }

    vector<string> res;
    string current;
    string line;
    while(getline(fasta,line)){
      if(line[0]=='>'){
	if(current.size()){
	  res.push_back(current);
	  current = "";
	}
	continue;
      }
      else {
	for(size_t i=0;i<line.size();i++){
	  char c = line[i];
	  if(isalpha(c)){
// 	    if(islower(c)){
// 	      c = toupper(c);
// 	    }
	    current += c;
	  }
	}
      }
    }
    res.push_back(current);
    return res;
  }

  vector<string> fasta2three(string s, bool frm_long){
    vector<string> res;
    for(int i=0;i<s.size();i++){
      string aa = one2three(s[i]);
      if(frm_long){
	char buff[256];
	sprintf(buff,"_%i",i+1);
	aa += buff;
      }
      res.push_back(aa);
    }
    return res;
  }

  bool is_protein(const PDBSegment & seg){
    if(three2one_map.size()==0)
      init = init_three2one();
    vector<string> seq = seg.sequence(1);
    for(size_t i=0;i<seq.size();i++){
      map<string,string>::iterator pos = three2one_map.find(seq[i]);
      if(pos==three2one_map.end()){
	//      cout<<"DEBUG NOT AN AMPINOACID "<<seq[i]<<endl;
	//      cout<<"DEBUG MAPSIZE "<<three2one_map.size()<<endl;
	return 0;
      }
    }
    return 1;
  }

  bool is_connected(const PDBSegment & seg){
    vector<double> X,Y,Z;
    vector<string> n;
    for(int i = 0; i< seg.size();i++){
      const PDBFragment & f = seg[i];
      for(int a = 0;a<f.size();a++){
	const PDBAtomLine & al = f[a];
	if(al.name=="N"||al.name=="CA"||al.name=="C"){
	  X.push_back(al.x);
	  Y.push_back(al.y);
	  Z.push_back(al.z);
	  char buff[256];
	  sprintf(buff,"%s_%i/%s",al.resName.c_str(),al.resSeq,al.name.c_str());
	  n.push_back(buff);
	}
      }
    }

    int s = X.size()-1;
    for(int i=0;i<s;i++){
      double dx = X[i]-X[i+1];
      double dy = Y[i]-Y[i+1];
      double dz = Z[i]-Z[i+1];
      double d2 = dx*dx+dy*dy+dz*dz;
      if(d2>9){
	//	aout<<"\t>>CHAIN BREAK BETWEEN "<<n[i]<<" "<<n[i+1]<<"\n";
	return false;
      }
    }
    return true;
  }

  PDBModel split(const PDBSegment & seg){
    //
    PDBModel mod;
    vector<double> X,Y,Z;
    vector<string> n;
    vector<int> frag;
    for(int i = 0; i< seg.size();i++){
      const PDBFragment & f = seg[i];
      for(int a = 0;a<f.size();a++){
	const PDBAtomLine & al = f[a];
	if(al.name=="N"||al.name=="CA"||al.name=="C"){
	  X.push_back(al.x);
	  Y.push_back(al.y);
	  Z.push_back(al.z);
	  char buff[256];
	  sprintf(buff,"%s",al.name.c_str());
	  n.push_back(buff);
	  frag.push_back(i);
	}
      }
    }

    PDBSegment * n_seg = new PDBSegment;
    int s = X.size()-1;
    for(int i=0;i<s;i++){
      if(n[i]!="C"||n[i+1]!="N") continue;
      double dx = X[i]-X[i+1];
      double dy = Y[i]-Y[i+1];
      double dz = Z[i]-Z[i+1];
      double d2 = dx*dx+dy*dy+dz*dz;
      if(d2>9){
	//	aout<<"\t>>CHAIN BREAK BETWEEN "<<n[i]<<" "<<n[i+1]<<"\n";
	n_seg->fragment.push_back(new PDBFragment(seg[frag[i]]));
	mod.segment.push_back(n_seg);
	n_seg = new PDBSegment;
      }
      else{
	n_seg->fragment.push_back(new PDBFragment(seg[frag[i]]));
      }
    }
    n_seg->fragment.push_back(new PDBFragment(seg[frag[s]]));
    mod.segment.push_back(n_seg);
    return mod;   
  }

}