/***************************************************************************
                          protein.cpp  -  description
                             -------------------
    begin                : Mon Jan 27 06:38:48 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 <almostexception.h>
#include <almosttypes.h>
#include <almosttemplate.h>
#include <pdb.h>
#include <mdb/mdb.h>
#include <moveset/bondrotation.h>
#include <algorithm>

#include <molecules/protein.h>
#include <buildutils.h>






Almost::Protein Almost::Protein::mutate(int aa, string frgm, MDB & mdb,
					string fpatch, string lpatch,
					bool keep_bb) const {
  vector<string> sequ = sequence(Protein::LONG);
  --aa;
  if((aa<fragment_offset())||(aa>=fragment_.size()+fragment_offset())){
    throw "Almost::Protein::mutate invalid amino-acid number\n";
  }
  //  aout<<"DEBUG  "<<sequ[aa]<<endl;
  sequ[aa-fragment_offset()] = frgm;
  Almost::Protein p("MUT");
  p.build_sequ(sequ,mdb,fpatch,lpatch);
  p.set_fragment_offset(fragment_offset());
  p.set_atom_offset(atom_offset());
  int n_wt = atom_offset();
  int n_mut = p.atom_offset();
  for(int i=0;i<sequ.size();i++){
    int fd_wt = i + fragment_offset();
    int fd_mut = i + p.fragment_offset(); 
    int wt_size = fragment_atoms(fd_wt).size();
    int mut_size = p.fragment_atoms(fd_mut).size();
    //     vector<int>::const_iterator 
    //       wt_iter = frgment_atoms(fd_wt).begin(),
    //       wt_end  = frgment_atoms(fd_wt).end(),
    //       mut_iter = p.frgment_atoms(fd_wt).begin(),
    //       mut__end  = p.frgment_atoms(fd_wt).end();
    if(i!=aa-fragment_offset()){
      for(int k=0;k<mut_size;k++){
	p[n_mut][0]=(*this)[n_wt][0];
	p[n_mut][1]=(*this)[n_wt][1];
	p[n_mut][2]=(*this)[n_wt][2];
	++n_mut;++n_wt;
      }
    }
    else {
      for(int k=0;k<mut_size;k++){
	p[n_mut][0]= 9999;
	p[n_mut][1]= 9999;
	p[n_mut][2]= 9999;
	++n_mut;
      }
      n_wt += wt_size;
    }
  }

  if(keep_bb){
    vector<int> s1,s2;
    s1 = fragment_atoms(aa);
    s2 = p.fragment_atoms(aa);
    int N1=-1,N2=-1,CA1=-1,CA2=-1,C1=-1,C2=-1,O1=-1,O2=-1;
    for(int q=0;q<s1.size();q++){
      if((*this)[s1[q]].name()=="N") N1 = s1[q];
      if((*this)[s1[q]].name()=="CA") CA1 = s1[q];
      if((*this)[s1[q]].name()=="C") C1 = s1[q];
      if((*this)[s1[q]].name()=="O") O1 = s1[q];
    }

    for(int q=0;q<s2.size();q++){
      if(p[s2[q]].name()=="N") N2 = s2[q];
      if(p[s2[q]].name()=="CA") CA2 = s2[q];
      if(p[s2[q]].name()=="C") C2 = s2[q];
      if(p[s2[q]].name()=="O") O2 = s2[q];
    }

    //
    if(N1!=-1&&N2!=-1){
      p[N2][0] = (*this)[N1][0];
      p[N2][1] = (*this)[N1][1];
      p[N2][2] = (*this)[N1][2];
    }


    if(CA1!=-1&&CA2!=-1){
      p[CA2][0] = (*this)[CA1][0];
      p[CA2][1] = (*this)[CA1][1];
      p[CA2][2] = (*this)[CA1][2];
    }

    if(C1!=-1&&C2!=-1){
      p[C2][0] = (*this)[C1][0];
      p[C2][1] = (*this)[C1][1];
      p[C2][2] = (*this)[C1][2];
    }

    if(O1!=-1&&O2!=-1){
      p[O2][0] = (*this)[O1][0];
      p[O2][1] = (*this)[O1][1];
      p[O2][2] = (*this)[O1][2];
    }

  }
  p.generate_missing(mdb,fpatch,lpatch);
  return p;
}



void Almost::Protein::disu_bond(string frag1, string frag2, MDB & mdb){

  //
  // FRAG 1 and 2 are CYS
  
  logMsg(Almost::ProteinPr,"Adding disulfite bond between "+frag1+" "+frag2);

  TemplatePatch disu = mdb.patch("DISU");

  //delete
  
  //atoms
  {
    vector<string> datom = disu.del_atoms();
  
    for(int i=0;i<datom.size();i++){
      string a = datom[i];
      int pos =0;
      while((pos<a.size())&&isdigit(a[pos]))  ++pos;
      int idx = -1;
      idx = atoi(string(a.begin(),a.begin()+pos).c_str());
      if(idx!=-1){
	if(idx==1){
	  int na = find_atom(frag1+"/"+string(a.begin()+pos,a.end()));
	  if(na!=1){
	    //	    aout<<"Deleting atom "<<na<<endl;
	    remove_atom(na);
	  }
	  else {
	    throw "Unable to find: "+frag1+"/"+string(a.begin()+pos,a.end());
	  }
	}
	else if(idx==2){
	  int na = find_atom(frag2+"/"+string(a.begin()+pos,a.end()));
	  if(na!=1){
	    //	    aout<<"Deleting atom "<<na<<endl;
	    remove_atom(na);
	  }
	  else {
	    throw "Unable to find: "+frag2+"/"+string(a.begin()+pos,a.end());
	  }
	}
      }
    }
  }
  //and so on ....

  //add-modify

  //atoms
  {
    vector<TAtom> atom = disu.atoms();
    for(int i=0;i<atom.size();i++){
      string a = atom[i].name();
      int idx = atoi(a.c_str());
      if(idx==0){
	throw "Wrong disulfide-patch";
      }
      if(idx==1){
	string a1 = frag1+"/"+string(a.begin()+1,a.end());
	int na = find_atom(a1);
	if(na==-1){
	  throw "Unable  to find atom: "+a1;
	}
	//	aout<<"Modifying "<<na<<" : "<<a1<<endl;
	//	atom_[na].set_name(atom[i].name());
	atom_[na].set_type_name(atom[i].type_name());
	atom_[na].set_charge(atom[i].charge());
      }
      else if(idx==2){
	string a2 = frag2+"/"+string(a.begin()+1,a.end());
	int na = find_atom(a2);
	if(na==-1){
	  throw "Unable  to find atom:"+a2;
	}
	//	aout<<"Modifying  "<<na<<" : "<<a2<<endl;
	//	atom_[na].set_name(atom[i].name());
	atom_[na].set_type_name(atom[i].type_name());
	atom_[na].set_charge(atom[i].charge());
      }
      else { 	throw "Wrong disulfide-patch"; }
    }
    
    //dihedrals
    {
      //auto
    }

    //dihedrals
    {
      vector<TDihedral> imph = disu.dihedrals();
      for(int i=0;i<imph.size();i++){
	int idx[4];
	string a[4];
	int na[4];
	vector<int> iv(4);

	for(int j=0;j<4;j++){
	  a[j] = imph[i].atom(j);
	  idx[j] = atoi(a[j].c_str());
	  if(idx[j]<=0) throw "Wrong disulfide-patch";
	  else if(idx[j]==1){
	    a[j] = frag1+"/"+string(a[j].begin()+1,a[j].end());
	  }
	  else if(idx[j]==2){
	    a[j] = frag2+"/"+string(a[j].begin()+1,a[j].end());
	  }
	  
	  na[j] = find_atom(a[j]);
	  if(na[j]==-1){
	    throw "Unable  to find atom: " +a[j];
	  }
	  
	  iv[j] = na[j];
	}

	if(!dihedral_.add_dihedral(iv))
	  aout<<"Warning can't add dihedral"<<endl;
      }
    }



    //impropers
    {
      vector<TImproper> imph = disu.impropers();
      for(int i=0;i<imph.size();i++){
	int idx[4];
	string a[4];
	int na[4];
	vector<int> iv(4);

	for(int j=0;j<4;j++){
	  a[j] = imph[i].atom(j);
	  idx[j] = atoi(a[j].c_str());
	  if(idx[j]<=0) throw "Wrong disulfide-patch";
	  else if(idx[j]==1){
	    a[j] = frag1+"/"+string(a[j].begin()+1,a[j].end());
	  }
	  else if(idx[j]==2){
	    a[j] = frag2+"/"+string(a[j].begin()+1,a[j].end());
	  }
	  
	  na[j] = find_atom(a[j]);
	  if(na[j]==-1){
	    throw "Unable  to find atom: " +a[j];
	  }
	  
	  iv[j] = na[j];
	}

	if(!improper_.add_improper(iv))
	  aout<<"Warning can't add improper"<<endl;
      }
    }
	
  }

  //bonds
  {
    vector<TBond> bond = disu.bonds();
    for(int i=0;i<bond.size();i++){
      string a1,a2;
      a1 = bond[i].first();
      a2 = bond[i].second();
      int idx1 = atoi(a1.c_str());
      int idx2 = atoi(a2.c_str());
      
      if((idx1<=0)||(idx2<=0)) throw "Wrong disulfide-patch";
      
      int na1,na2;
      
      if(idx1==1){
	a1 = frag1+"/"+string(a1.begin()+1,a1.end());
	na1 = find_atom(a1);
	if(na1==-1){
	  throw "Unable  to find atom: "+a1;
	}
      }
      else if(idx1==2){
	a1 = frag2+"/"+string(a1.begin()+1,a1.end());
	na1 = find_atom(a1);
	if(na1==-1){
	  throw "Unable  to find atom: "+a1;
	}
      }

      if(idx2==1){
	a2 = frag1+"/"+string(a2.begin()+1,a2.end());
	na2 = find_atom(a2);
	if(na2==-1){
	  throw "Unable  to find atom: "+a2;
	}
      }
      else if(idx2==2){
	a2 = frag2+"/"+string(a2.begin()+1,a2.end());
	na2 = find_atom(a2);
	if(na2==-1){
	  throw "Unable  to find atom: "+a2;
	}
      }

      int kind = bond[i].kind();

      if(bonded(na1,na2)){
	//aout<<"Deleting old bond between "<<na1<<" "<<na2<<endl;
	remove_bond(na1,na2);
      }
      if(!add_bond(na1,na2,kind)){
	aout<<"Couldn't add bond"<<endl;
      }
      else { 
	//aout<<"DISU bond added"<<endl;
      }
    }
  }


  //regenerate angles
  angle_generate();
  if(disu.dihedrals().size()==0) dihedral_generate();
  logMsg(Almost::ProteinPr,"Done");
}


void Almost::Protein::auto_disu_bonds(double max_dist,MDB & mdb){
  logMsg(Almost::ProteinPr,
	 "autogenerating disulfiede bonds for protein "+name_);
  vector<int> cysd;
  for(int i=0;i<fragment_.size();i++){
    int pos = fragment_[i].find("_");
    string base;
    if(pos!=string::npos){
      base = string(fragment_[i].begin(),fragment_[i].begin()+pos);
    }
    if(base=="CYS"){
      cysd.push_back(i);
    }
  }


  //compute distance between S atoms ....
  for(int i=0;i<cysd.size();i++){
    for(int j=i+1;j<cysd.size();j++){
      int s1 = find_atom(fragment_[cysd[i]]+"/"+"SG")-atom_offset_;
      int s2 = find_atom(fragment_[cysd[j]]+"/"+"SG")-atom_offset_;
      double d2 = 0;
      for(int k=0;k<3;k++){
	d2 += (atom_[s1][k]-atom_[s2][k])*(atom_[s1][k]-atom_[s2][k]);
      }
      if(d2<(max_dist*max_dist)){
	disu_bond(fragment_[cysd[i]],fragment_[cysd[j]],mdb);
      }
    }
  }
}



 
void Almost::Protein::linearize(){
  BondRotation<double> bond_rotation(*this);
  vector<pair<int,string> > atoms;
  vector<double*> coor;
  for(int i=0;i<atom_.size();i++)
    coor.push_back(const_cast<double*>(atom_[i].coor()));

  for(int i=0;i<atom_.size();++i){
    string name = atom_[i].name();
    if((name=="N")||
       (name=="C")||
       (name=="CA"))
      atoms.push_back(pair<int,string>(i,name));
  }
  
  for(int i=0;i<atoms.size()-3;i++){
    string dihe_name = atoms[i].second;
    for(int j=0;j<3;++j)
      dihe_name +="-"+atoms[i+j+1].second;
    if((dihe_name=="N-CA-C-N")||(dihe_name=="C-N-CA-C")){
      bond_rotation.set_dihedral(coor,
				 atoms[i].first,
				 atoms[i+1].first,
				 atoms[i+2].first,
				 atoms[i+3].first,
				 M_PI);
      bond_rotation.set_dihedral(coor,
				 atoms[i].first,
				 atoms[i+1].first,
				 atoms[i+2].first,
				 atoms[i+3].first,
				 M_PI);
    }
  }

  //Side chains ...
  //   for(int i=0;i<fragment_.size();++i){
  //     string fragment = fragment_name_[i];
  //     vector<int> atoms = fragment_atoms_[i];
  //     int N = -1;
  //     int CA= -1;
  //     int CB= -1;
  //     int SC= -1; //"first" atom in side chain;
      
  //     for(int j=0;j<atoms.size();++j){
  //       if(atom_[atoms[j]].name()=="N") N=atoms[j];
  //       if(atom_[atoms[j]].name()=="CA") CA=atoms[j];
  //       if(atom_[atoms[j]].name()=="CB") CB=atoms[j];
  //     }
      
  //     if((N==-1)||
  //        (CA==-1)||
  //        (CB==-1)) continue;

  //     //find "first" atom in side chain;
  //     vector<int> out_bonds = this->bonds(CB,Molecules::OUT);
      
  //     int j=0;
  //     while(j<out_bonds.size()){
  //       int out_atom = out_bonds[j++];
  //       if((out_atom != N)&&
  //      	 (out_atom != CA)&&
  //      	 (out_atom != CB) ){
  //      	SC = out_atom;
  //      	break;
  //       }
  //     }

  //     if(SC!=-1){
  //       bond_rotation.set_dihedral(coor_manager.atomCoor,N,CA,CB,SC,0);
  //       bond_rotation.set_dihedral(coor_manager.atomCoorBuffer,N,CA,CB,SC,0);
  //     }
  //   }

}


vector<double> Almost::Protein::phi() const {
  vector<int> N;
  vector<int> C;
  vector<int> CA;
  int b = atom_offset();
  int e = atom_size()+atom_offset();
  for(int i = b; i<e;i++){
    if((*this)[i].name()=="N")  N.push_back(i);
    if((*this)[i].name()=="C")  C.push_back(i);
    if((*this)[i].name()=="CA") CA.push_back(i);
  }

  vector<double> t;
  //first is undef == 999
  t.push_back(999);
  for(int i=1;i<fragment_size();i++){
    const  double * c[4];
    c[0] = (*this)[C[i-1]].coor();
    c[1] = (*this)[N[i]].coor();
    c[2] = (*this)[CA[i]].coor();
    c[3] = (*this)[C[i]].coor();
    double a = Kernel::Geometry::dihedral(c[0],c[1],c[2],c[3]);
    t.push_back(a);
  }
  return t;
}

vector<vector<int> > Almost::Protein::phi_atoms() const {
  vector<int> N;
  vector<int> C;
  vector<int> CA;
  int b = atom_offset();
  int e = atom_size()+atom_offset();
  for(int i = b; i<e;i++){
    if((*this)[i].name()=="N")  N.push_back(i);
    if((*this)[i].name()=="C")  C.push_back(i);
    if((*this)[i].name()=="CA") CA.push_back(i);
  }

  vector<vector<int> > t;
  t.resize(fragment_size());
  //first is undef == 999
  for(int i=1;i<fragment_size();i++){
    int c[4];
    c[0] = C[i-1];
    c[1] = N[i];
    c[2] = CA[i];
    c[3] = C[i];
    vector<int> v;
    for(int q=0;q<4;q++) v.push_back(c[q]);
    t[i] = v;
  }
  return t;
}


vector<double> Almost::Protein::psi() const {
  vector<int> N;
  vector<int> C;
  vector<int> CA;
  int b = atom_offset();
  int e = atom_size()+atom_offset();
  for(int i = b; i<e;i++){
    if((*this)[i].name()=="N")  N.push_back(i);
    if((*this)[i].name()=="C")  C.push_back(i);
    if((*this)[i].name()=="CA") CA.push_back(i);
  }

  vector<double> t;
  for(int i=0;i<fragment_size()-1;i++){
    const  double * c[4];
    c[0] = (*this)[N[i]].coor();
    c[1] = (*this)[CA[i]].coor();
    c[2] = (*this)[C[i]].coor();
    c[3] = (*this)[N[i+1]].coor();
    double a = Kernel::Geometry::dihedral(c[0],c[1],c[2],c[3]);
    t.push_back(a);
  }
  t.push_back(999);
  return t;
}

vector<vector<int> > Almost::Protein::psi_atoms() const {
  vector<int> N;
  vector<int> C;
  vector<int> CA;
  int b = atom_offset();
  int e = atom_size()+atom_offset();
  for(int i = b; i<e;i++){
    if((*this)[i].name()=="N")  N.push_back(i);
    if((*this)[i].name()=="C")  C.push_back(i);
    if((*this)[i].name()=="CA") CA.push_back(i);
  }

  vector<vector<int> > t;
  t.resize(fragment_size());
  for(int i=0;i<fragment_size()-1;i++){
    int c[4];
    c[0] = C[i-1];
    c[1] = N[i];
    c[2] = CA[i];
    c[3] = C[i];
    vector<int> v;
    for(int q=0;q<4;q++) v.push_back(c[q]);
    t[i] = v;    
  }

  return t;
}

vector<double> Almost::Protein::omega() const {
  vector<int> N;
  vector<int> C;
  vector<int> CA;
  int b = atom_offset();
  int e = atom_size()+atom_offset();
  for(int i = b; i<e;i++){
    if((*this)[i].name()=="N")  N.push_back(i);
    if((*this)[i].name()=="C")  C.push_back(i);
    if((*this)[i].name()=="CA") CA.push_back(i);
  }

  vector<double> t;
  for(int i=0;i<fragment_size()-1;i++){
    const  double * c[4];
    c[0] = (*this)[CA[i]].coor();
    c[1] = (*this)[C[i]].coor();
    c[2] = (*this)[N[i+1]].coor();
    c[3] = (*this)[CA[i+1]].coor();
    double a = Kernel::Geometry::dihedral(c[0],c[1],c[2],c[3]);
    t.push_back(a);
  }
  t.push_back(999);
  return t;
}

vector<vector<int> > Almost::Protein::omega_atoms() const {
  vector<int> N;
  vector<int> C;
  vector<int> CA;
  int b = atom_offset();
  int e = atom_size()+atom_offset();
  for(int i = b; i<e;i++){
    if((*this)[i].name()=="N")  N.push_back(i);
    if((*this)[i].name()=="C")  C.push_back(i);
    if((*this)[i].name()=="CA") CA.push_back(i);
  }

  vector<vector<int> > t;
  t.resize(fragment_size());
  for(int i=0;i<fragment_size()-1;i++){
    int c[4];
    c[0] = C[i-1];
    c[1] = N[i];
    c[2] = CA[i];
    c[3] = C[i];
    vector<int> v;
    for(int q=0;q<4;q++) v.push_back(c[q]);
    t[i] = v;
  }

  return t;
}

inline bool check_index(int a, int b, int c, int d){
  if(a==-1) return false;
  if(b==-1) return false;
  if(c==-1) return false;
  if(d==-1) return false;
  return true;
}

#define CHI(a,b,c,d,res) {						\
    if(check_index(a,b,c,d)){						\
      double chi = Kernel::Geometry::dihedral(p[a].coor(),		\
					      p[b].coor(),		\
					      p[c].coor(),		\
					      p[d].coor());		\
      rot.push_back(chi);						\
    }									\
    else{ cout<<">> Missing atom "<<res<<"\n"; rot.push_back(999);}	\
  }



vector<double> rotamer_ARG(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD_ = -1, NE_ = -1,
    CZ_ = -1, NH1_ = -1, NH2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD")
      CD_= atom[i];
    if(p[atom[i]].name() =="NE")
      NE_ = atom[i];
    if(p[atom[i]].name() =="CZ")
      CZ_= atom[i];
    if(p[atom[i]].name() =="NH1")
      NH1_= atom[i];
    if(p[atom[i]].name() =="NH2")
      NH2_= atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"ARG");
  CHI(CA_,CB_,CG_,CD_,"ARG");
  CHI(CB_,CG_,CD_,NE_,"ARG");
  CHI(CG_,CD_,NE_,CZ_,"ARG");
  CHI(CD_,NE_,CZ_,NH1_,"ARG");
  CHI(CD_,NE_,CZ_,NH2_,"ARG"); 
  return rot;
}


vector<double> rotamer_ASN(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, OD1_ = -1, ND2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="OD1")
      OD1_= atom[i];
    if(p[atom[i]].name() =="ND2")
      ND2_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"ASN");
  CHI(CA_,CB_,CG_,OD1_,"ASN");
  CHI(CA_,CB_,CG_,ND2_,"ASN");
  return rot;
}

vector<double> rotamer_ASP(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, OD1_ = -1, OD2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="OD1")
      OD1_= atom[i];
    if(p[atom[i]].name() =="OD2")
      OD2_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"ASP");
  CHI(CA_,CB_,CG_,OD1_,"ASP");
  CHI(CA_,CB_,CG_,OD2_,"ASP");
  return rot;
}

vector<double> rotamer_CYS(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, SG_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="SG")
      SG_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,SG_,"CYS");
  return rot;
}


vector<double> rotamer_GLN(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD_ = -1, OE1_ = -1,
    NE2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD")
      CD_= atom[i];
    if(p[atom[i]].name() =="OE1")
      OE1_ = atom[i];
    if(p[atom[i]].name() =="NE2")
      NE2_= atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"GLN");
  CHI(CA_,CB_,CG_,CD_,"GLN");
  CHI(CB_,CG_,CD_,OE1_,"GLN");
  CHI(CB_,CG_,CD_,NE2_,"GLN");
  return rot;
}

vector<double> rotamer_GLU(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD_ = -1, OE1_ = -1,
    OE2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD")
      CD_= atom[i];
    if(p[atom[i]].name() =="OE1")
      OE1_ = atom[i];
    if(p[atom[i]].name() =="OE2")
      OE2_= atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"GLU");
  CHI(CA_,CB_,CG_,CD_,"GLU");
  CHI(CB_,CG_,CD_,OE1_,"GLU");
  CHI(CB_,CG_,CD_,OE2_,"GLU");
  return rot;
}

vector<double> rotamer_HIS(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, ND1_ = -1, CD2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="ND1")
      ND1_= atom[i];
    if(p[atom[i]].name() =="CD2")
      CD2_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"HIS");
  CHI(CA_,CB_,CG_,ND1_,"HIS");
  CHI(CA_,CB_,CG_,CD2_,"HIS");
  return rot;
}

vector<double> rotamer_ILE(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG1_ = -1, CG2_ = -1, CD1_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG1")
      CG1_ = atom[i];
    if(p[atom[i]].name() =="CG2")
      CG2_= atom[i];
    if(p[atom[i]].name() =="CD"||p[atom[i]].name() =="CD1")
      CD1_ = atom[i];

  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG1_,"ILE");
  CHI(N_,CA_,CB_,CG2_,"ILE");
  CHI(CA_,CB_,CG1_,CD1_,"ILE");
  return rot;
}

vector<double> rotamer_LEU(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD1_ = -1, CD2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD1")
      CD1_= atom[i];
    if(p[atom[i]].name() =="CD2")
      CD2_ = atom[i];

  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"LEU");
  CHI(CA_,CA_,CB_,CD1_,"LEU");
  CHI(CA_,CB_,CB_,CD2_,"LEU");
  return rot;
}

vector<double> rotamer_LYS(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD_ = -1, CE_ = -1,
    NZ_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD")
      CD_= atom[i];
    if(p[atom[i]].name() =="CE")
      CE_ = atom[i];
    if(p[atom[i]].name() =="NZ")
      NZ_= atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"LYS");
  CHI(CA_,CB_,CG_,CD_,"LYS");
  CHI(CB_,CG_,CD_,CE_,"LYS");
  CHI(CG_,CD_,CE_,NZ_,"LYS");
  return rot;
}

vector<double> rotamer_MET(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, SD_ = -1, CE_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="SD")
      SD_= atom[i];
    if(p[atom[i]].name() =="CE")
      CE_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"MET");
  CHI(CA_,CB_,CG_,SD_,"MET");
  CHI(CB_,CG_,SD_,CE_,"MET");
  return rot;
}


vector<double> rotamer_PHE(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD1_ = -1, CD2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD1")
      CD1_= atom[i];
    if(p[atom[i]].name() =="CD2")
      CD2_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"PHE");
  CHI(CA_,CA_,CB_,CD1_,"PHE");
  CHI(CA_,CB_,CB_,CD2_,"PHE");
  return rot;
}

//PRO???



vector<double> rotamer_SER(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, OG_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="OG")
      OG_ = atom[i];
  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,OG_,"SER");
  return rot;
}


vector<double> rotamer_THR(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, OG1_ = -1, CG2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="OG1")
      OG1_ = atom[i];
    if(p[atom[i]].name() =="CG2")
      CG2_= atom[i];

  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,OG1_,"THR");
  CHI(N_,CA_,CB_,CG2_,"THR");
  return rot;
}


vector<double> rotamer_TRP(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD1_ = -1, CD2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD1")
      CD1_= atom[i];
    if(p[atom[i]].name() =="CD2")
      CD2_ = atom[i];

  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"TRP");
  CHI(CA_,CA_,CB_,CD1_,"TRP");
  CHI(CA_,CB_,CB_,CD2_,"TRP");
  return rot;
}

vector<double> rotamer_TYR(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG_ = -1, CD1_ = -1, CD2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG")
      CG_ = atom[i];
    if(p[atom[i]].name() =="CD1")
      CD1_= atom[i];
    if(p[atom[i]].name() =="CD2")
      CD2_ = atom[i];

  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG_,"TYR");
  CHI(CA_,CA_,CB_,CD1_,"TYR");
  CHI(CA_,CB_,CB_,CD2_,"TYR");
  return rot;
}

vector<double> rotamer_VAL(const Almost::Protein & p,
			   const vector<int> & atom, bool verbose){
  int N_ = -1, CA_ = -1, CB_ = -1, CG1_ = -1, CG2_ = -1;


  for(int i=0;i<atom.size();i++){
    if(p[atom[i]].name() == "N")
      N_ = atom[i];
    if(p[atom[i]].name() =="CA")
      CA_ = atom[i];
    if(p[atom[i]].name() =="CB")
      CB_= atom[i];
    if(p[atom[i]].name() =="CG1")
      CG1_ = atom[i];
    if(p[atom[i]].name() =="CG2")
      CG2_= atom[i];


  }
  vector<double> rot;
  //CHI1
  CHI(N_,CA_,CB_,CG1_,"VAL");
  CHI(N_,CA_,CB_,CG2_,"VAL");
  return rot;
}

vector<double> Almost::Protein::rotamer(int aa, bool verbose) const{
  /*           chi1        chi2         chi3         chi4
	       ARG      N-CA-CB-CG  CA-CB-CG-CD  CB-CG-CD-NE  CG-CD-NE-CZ
	       ASN      N-CA-CB-CG  CA-CB-CG-ND2(?)(OD1?)
	       ASP      N-CA-CB-CG  CA-CB-CG-OD1
	       CYS      N-CA-CB-SG
	       GLN      N-CA-CB-CG  CA-CB-CG-CD  CB-CG-CD-NE2/OE1
	       GLU      N-CA-CB-CG  CA-CB-CG-CD  CB-CG-CD-OE1
	       HIS      N-CA-CB-CG  CA-CB-CG-ND1/CD2(?)
	       ILE      N-CA-CB-CG(*) ???
	       LEU      N-CA-CB-CG  CA-CB-CG-CD1/CD2
	       LYS      N-CA-CB-CG  CA-CB-CG-CD  CB-CG-CD-CE  CG-CD-CE-NZ
	       MET      N-CA-CB-CG  CA-CB-CG-SD  CB-CG-SD-CE
	       PHE      N-CA-CB-CG  CA-CB-CG-CD1/2
	       PRO      ??
	       SER      N-CA-CB-OG
	       THR      N-CA-CB-OG1(*)
	       TRP      N-CA-CB-CG  CA-CB-CG-CD1/2
	       TYR      N-CA-CB-CG  CA-CB-CG-CD1/CD2(?)
	       VAL      N-CA-CB-CG
  */

  string fname = fragment_name(aa,Protein::SHORT);
  vector<int> fatom = fragment_atoms(aa);
  //   {
  //     for(int i=0;i<fatom.size();i++) cout<<fatom[i]<<" ";
  //     cout<<"\n";
  //   }
  if(fname=="ARG") return rotamer_ARG(*this,fatom,verbose);
  if(fname=="ASN") return rotamer_ASN(*this,fatom,verbose); 
  if(fname=="ASP") return rotamer_ASP(*this,fatom,verbose); 
  if(fname=="CYS") return rotamer_CYS(*this,fatom,verbose); 
  if(fname=="GLN") return rotamer_GLN(*this,fatom,verbose); 
  if(fname=="GLU") return rotamer_GLU(*this,fatom,verbose); 
  if(fname=="HSP") return rotamer_HIS(*this,fatom,verbose); 
  if(fname=="ILE") return rotamer_ILE(*this,fatom,verbose); 
  if(fname=="LEU") return rotamer_LEU(*this,fatom,verbose); 
  if(fname=="LYS") return rotamer_LYS(*this,fatom,verbose); 
  if(fname=="MET") return rotamer_MET(*this,fatom,verbose); 
  if(fname=="PHE") return rotamer_PHE(*this,fatom,verbose); 
  //   if(fname=="PRO") return rotamer_PRO(*this,fatom,verbose); 
  if(fname=="SER") return rotamer_SER(*this,fatom,verbose); 
  if(fname=="THR") return rotamer_THR(*this,fatom,verbose); 
  if(fname=="TRP") return rotamer_TRP(*this,fatom,verbose); 
  if(fname=="TYR") return rotamer_TYR(*this,fatom,verbose); 
  if(fname=="VAL") return rotamer_VAL(*this,fatom,verbose); 

  return vector<double>();
}