/***************************************************************************
                          mft_mod.cpp  -  description
                             -------------------
    begin                : Tue May 15 08:10:32 2007
    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 <zeta/module.h>
#include <zeta/class.h>
#include <zeta/zextclass.h>
#include <zeta/zobj.hpp>
#include <molecules.h>
#include <select.h>
#include <iterator>
#include <secstruct.h>
#ifdef __MFT__SC__
#include <ff.h>
#include <moveset/bondrotation.h>
#include <rose/dunbracklib.h>
#include <forcefield/const/envconst.h>
#include <forcefield/const/envpair.h>
#endif
#include <pdb.h>
using namespace Almost;

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

ALMTOSTRING(Molecules);


extern vector<string> aminoacids();
extern void prot2pdb(const Protein & prot, string ofile);
extern void mol2pdb(const Molecules & prot, string ofile);
namespace mft {
#ifdef __MFT__SC__
  vector<int> get_atoms(const Molecules & molecules,
			string fname,
			string a1,
			string a2,
			string a3,
			string a4){

    vector<int> a; a.resize(4);
    for(int i=0;i<4;i++) a[i] = -1;

    {
      NameSelection s1(molecules,fname+"/"+a1);
      if(s1.result().size()==1) a[0] = *s1.result().begin();
      else cout<<"Error in selection of "<<fname+"/"+a1<<endl;
    }

    {
      NameSelection s1(molecules,fname+"/"+a2);
      if(s1.result().size()==1) a[1] = *s1.result().begin();
      else cout<<"Error in selection of "<<fname+"/"+a2<<endl;
    }

    {
      NameSelection s1(molecules,fname+"/"+a3);
      if(s1.result().size()==1) a[2] = *s1.result().begin();
      else cout<<"Error in selection of "<<fname+"/"+a3<<endl;
    }

    {
      NameSelection s1(molecules,fname+"/"+a4);
      if(s1.result().size()==1) a[3] = *s1.result().begin();
      else cout<<"Error in selection of "<<fname+"/"+a4<<endl;
    }

    return a;
  }
  
  bool is_ok(vector<int> a){
    for(int i=0;i<4;i++) if(a[i]<0) return false;
    return true;
  }


  vector<vector<int> > get_rotamers_ala(const Molecules & molecules, string fname){ return vector<vector<int> >();
  }

  vector<vector<int> > get_rotamers_gly(const Molecules & molecules, string fname){ return vector<vector<int> >();
  }

  vector<vector<int> > get_rotamers_arg(const Molecules & molecules, string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CB","CG","CD","NE");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CG","CD","NE","CZ");
      if(is_ok(a)) res.push_back(a);	
    }
    return res;
  }

  vector<vector<int> > get_rotamers_asn(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","OD1");//??
      if(is_ok(a)) res.push_back(a);	
    }
    return res;
  }

  vector<vector<int> > get_rotamers_asp(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","OD1");//??
      if(is_ok(a)) res.push_back(a);	
    }
    return res;
  }

  vector<vector<int> > get_rotamers_cys(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","SG");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_gln(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CB","CG","CD","OE1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_glu(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
    if(is_ok(a)) res.push_back(a);

    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CB","CG","CD","OE1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_his(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
    if(is_ok(a)) res.push_back(a);

    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","ND1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_ile(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG1");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG1","CD");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_leu(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_lys(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CB","CG","CD","CE");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CG","CD","CE","NZ");
      if(is_ok(a)) res.push_back(a);	
    }
    return res;
  }

  vector<vector<int> > get_rotamers_met(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","SD");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CB","CG","SD","CE");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_phe(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_pro(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    return res;
  }
    

  vector<vector<int> > get_rotamers_ser(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","OG");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_thr(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","OG1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }


  vector<vector<int> > get_rotamers_trp(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_tyr(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG");
      if(is_ok(a)) res.push_back(a);
    }
    {
      vector<int> a = get_atoms(molecules,fname,"CA","CB","CG","CD1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }

  vector<vector<int> > get_rotamers_val(const Molecules & molecules,string fname){
    vector<vector<int> > res;
    {
      vector<int> a = get_atoms(molecules,fname,"N","CA","CB","CG1");
      if(is_ok(a)) res.push_back(a);
    }
    return res;
  }


  vector<vector<int> > get_rotamers(const Molecules & m,
				    string lname, string sname){


    if(sname=="ALA")
      return get_rotamers_ala(m,lname);
    if(sname=="CYS")
      return get_rotamers_cys(m,lname);
    if(sname=="ASP")
      return get_rotamers_asp(m,lname);
    if(sname=="GLU")
      return get_rotamers_glu(m,lname);
    if(sname=="PHE")
      return get_rotamers_phe(m,lname);
    if(sname=="GLY")
      return get_rotamers_gly(m,lname);
    if(sname=="HIS")
      return get_rotamers_his(m,lname);
    if(sname=="ILE")
      return get_rotamers_ile(m,lname);
    if(sname=="LYS")
      return get_rotamers_lys(m,lname);
    if(sname=="LEU")
      return get_rotamers_leu(m,lname);
    if(sname=="MET")
      return get_rotamers_met(m,lname);
    if(sname=="ASN")
      return get_rotamers_asn(m,lname);
    if(sname=="PRO")
      return get_rotamers_pro(m,lname);
    if(sname=="GLN")
      return get_rotamers_gln(m,lname);
    if(sname=="ARG")
      return get_rotamers_arg(m,lname);
    if(sname=="SER")
      return get_rotamers_ser(m,lname);
    if(sname=="THR")
      return get_rotamers_thr(m,lname);
    if(sname=="VAL")
      return get_rotamers_val(m,lname);
    if(sname=="TRP")
      return get_rotamers_trp(m,lname);
    if(sname=="TYR")
      return get_rotamers_tyr(m,lname);

    return vector<vector<int> >();
  }

  void assign_backbone(Protein &trg, const Protein  & src){
    int size = trg.fragment_size();
    for(int aa=0;aa<size;aa++){
      vector<int> s1,s2;
      s1 = src.fragment_atoms(aa);
      s2 = trg.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(src[s1[q]].name()=="N") N1 = s1[q];
	if(src[s1[q]].name()=="CA") CA1 = s1[q];
	if(src[s1[q]].name()=="C") C1 = s1[q];
	if(src[s1[q]].name()=="O") O1 = s1[q];
      }

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

      //
      if(N1!=-1&&N2!=-1){
	trg[N2][0] = src[N1][0];
	trg[N2][1] = src[N1][1];
	trg[N2][2] = src[N1][2];
      }
      
      
      if(CA1!=-1&&CA2!=-1){
	trg[CA2][0] = src[CA1][0];
	trg[CA2][1] = src[CA1][1];
	trg[CA2][2] = src[CA1][2];
      }
      
      if(C1!=-1&&C2!=-1){
	trg[C2][0] = src[C1][0];
	trg[C2][1] = src[C1][1];
	trg[C2][2] = src[C1][2];
      }
      
      if(O1!=-1&&O2!=-1){
	trg[O2][0] = src[O1][0];
	trg[O2][1] = src[O1][1];
	trg[O2][2] = src[O1][2];
      }            
    }
  }

  void self_energy(const Molecules & m1,
		   int aa1,
		   const EnergyOptions<double> & opts,
		   DunbrackLib & db,
		   MDB & mdb, ostream &out){

    int phi1 = m1.protein(0).phi()[aa1]*180.0/M_PI;
    int psi1 = m1.protein(0).psi()[aa1]*180.0/M_PI;
    string s1 = m1.fragment_name(aa1,Molecules::SHORT);
    {
      //      stringstream str1;
      //      str1<<"//"<<m1.fragment_name(aa1,Molecules::SHORT)<<"_"<<aa1+1;
      int size = db.size(s1,phi1,psi1);
      out<<"Size "<<s1<<" "<<aa1<<" "<<max(int(1),size)<<"\n";
    }
    Molecules mols;
    {
      Protein p1("P1");
      string fpatch = "NONE";
      string lpatch = "NONE";
      if(aa1==0) fpatch = "DEFAULT";
      if(aa1==m1.fragment_size()-1) lpatch = "DEFAULT";
      vector<string> seq;
      stringstream str;
      str<<m1.fragment_name(aa1,Molecules::SHORT)<<"_"<<aa1+1;
      seq.push_back(str.str());
      p1.build_template(seq,mdb,fpatch,lpatch);
      vector<int> atm = m1.fragment_atoms(aa1);
      for(int i=0;i<atm.size();i++){
	p1[i][0] = m1[atm[i]][0];
	p1[i][1] = m1[atm[i]][1];
	p1[i][2] = m1[atm[i]][2];
      }
      mols.add_protein(p1);
    }
    map<double,DunbrackLib::Rotamer>::reverse_iterator i1,i2,e1,e2;
    Coor<double> coor(mols);
    stringstream str1;
    str1<<"//"<<m1.fragment_name(aa1,Molecules::SHORT)<<"_"<<aa1+1;
    vector<vector<int> > r1 = get_rotamers(mols,str1.str(),s1);
    BondRotation<double> br(mols,0);
    ForceField<double> ff(mols,opts,mdb);
    i1 = db.rbegin(s1,phi1,psi1);
    e1 = db.rend(s1,phi1,psi1);
    if(r1.size()){
      while(i1!=e1){
	//set the first
	double *chi1 =  i1->second.chi;
	for(int c = 0; c < r1.size();c++){
	  br.set_dihedral(coor,
			  r1[c][0],
			  r1[c][1],
			  r1[c][2],
			  r1[c][3],
			  chi1[c]*(M_PI/180.0));
	}
	EnergyData<double> e = ff.energy(coor); 
	out<<e["energy"]<<" "<<e["vdw"]<<" "<<e["elec"]<<" "<<e["solv"]<<"\n";
	++i1;
      }
    }
    if(r1.size()==0){
      EnergyData<double> e = ff.energy(coor); 
      out<<e["energy"]<<" "<<e["vdw"]<<" "<<e["elec"]<<" "<<e["solv"]<<"\n";
    }
  }

  void inter_energy(const Molecules & m1, const Molecules &m2,
		    int aa1, int aa2, const EnergyOptions<double> & opts,
		    DunbrackLib & db,
		    MDB & mdb, ostream &out){

    //Build 2 facke proteins;
    Molecules mols;
    {
      Protein p1("P1");
      string fpatch = "NONE";
      string lpatch = "NONE";
      if(aa1==0) fpatch = "DEFAULT";
      if(aa1==m1.fragment_size()-1) lpatch = "DEFAULT";
      vector<string> seq;
      stringstream str;
      str<<m1.fragment_name(aa1,Molecules::SHORT)<<"_"<<aa1+1;
      seq.push_back(str.str());
      p1.build_template(seq,mdb,fpatch,lpatch);
      vector<int> atm = m1.fragment_atoms(aa1);
      for(int i=0;i<atm.size();i++){
	p1[i][0] = m1[atm[i]][0];
	p1[i][1] = m1[atm[i]][1];
	p1[i][2] = m1[atm[i]][2];
      }
      mols.add_protein(p1);;
    }
    {
      Protein p2("P2");
      string fpatch = "NONE";
      string lpatch = "NONE";
      if(aa2==0) fpatch = "DEFAULT";
      if(aa2==m2.fragment_size()-1) lpatch = "CTER";
      //      cout<<"DEBUG PATCH "<<fpatch<<" "<<lpatch<<endl;
      vector<string> seq;
      stringstream str;
      str<<m2.fragment_name(aa2,Molecules::SHORT)<<"_"<<aa2+1;
      seq.push_back(str.str());
      p2.build_template(seq,mdb,fpatch,lpatch);
      vector<int> atm = m2.fragment_atoms(aa2);
      //      cout<<"DEBUG "<<atm.size()<<" "<<p2.atom_size()<<endl;
      for(int i=0;i<atm.size();i++){
	p2[i][0] = m2[atm[i]][0];
	p2[i][1] = m2[atm[i]][1]; 
	p2[i][2] = m2[atm[i]][2];
      }
      mols.add_protein(p2);;
    }

    //Build ff
    ForceField<double> ff(mols,opts,mdb);
    //    mol2pdb(mols,"db.pdb");
    //Iterate thru rotamers
    BondRotation<double> br(mols,0);
    int phi1 = m1.protein(0).phi()[aa1]*180.0/M_PI;
    int psi1 = m1.protein(0).psi()[aa1]*180.0/M_PI;
    string s1 = m1.fragment_name(aa1,Molecules::SHORT);
    int phi2 = m2.protein(0).phi()[aa2]*180.0/M_PI;
    int psi2 = m2.protein(0).psi()[aa2]*180.0/M_PI;
    string s2 = m2.fragment_name(aa2,Molecules::SHORT);
    //Remember to treat special cases with 0 rotames!!!!
    map<double,DunbrackLib::Rotamer>::reverse_iterator i1,i2,e1,e2;
    Coor<double> coor(mols);
    stringstream str1,str2;
    str1<<"//"<<m1.fragment_name(aa1,Molecules::SHORT)<<"_"<<aa1+1;
    str2<<"//"<<m2.fragment_name(aa2,Molecules::SHORT)<<"_"<<aa2+1;
    vector<vector<int> > r1 = get_rotamers(mols,str1.str(),s1);
    vector<vector<int> > r2 = get_rotamers(mols,str2.str(),s2);
    
    i1 = db.rbegin(s1,phi1,psi1);
    i2 = db.rbegin(s2,phi2,psi2);
    e1 = db.rend(s1,phi1,psi1);
    e2 = db.rend(s2,phi2,psi2);
    //    cout<<"DEBUG apply rots"<<endl;
    int count1 = 0;
    int count2 = 0;
    if(r1.size()&&r2.size()){
      while(i1!=e1){
	//set the first
	double *chi1 =  i1->second.chi;
	for(int c = 0; c < r1.size();c++){
	  br.set_dihedral(coor,
			  r1[c][0],
			  r1[c][1],
			  r1[c][2],
			  r1[c][3],
			  chi1[c]*(M_PI/180.0));
	}
	i2 = db.rbegin(s2,phi2,psi2);
	count2 = 0;
	while(i2!=e2){
	  //Apply
	  double * chi2 = i2->second.chi;
	  for(int c = 0; c < r2.size();c++){
	    br.set_dihedral(coor,
			    r2[c][0],
			    r2[c][1],
			    r2[c][2],
			    r2[c][3],
			    chi2[c]*(M_PI/180.0));
	  }
	  out<<s1<<" "<<aa1<<" "<<" "<<count1<<" "
	     <<s2<<" "<<aa2<<" "<<count2<<" ";
	  EnergyData<double> e = ff.energy(coor); 
	  out<<e["energy"]<<" "<<e["vdw"]<<" "<<e["elec"]<<" "<<e["solv"]<<"\n";
	  ++count2;
	  ++i2;
	}
	++count1;
	++i1;
      }
    }

    if(r1.size()==0&&r2.size()){
      i2 = db.rbegin(s2,phi2,psi2);
      while(i2!=e2){
	//Apply
	double * chi2 = i2->second.chi;
	for(int c = 0; c < r2.size();c++){
	  br.set_dihedral(coor,
			  r2[c][0],
			  r2[c][1],
			  r2[c][2],
			  r2[c][3],
			  chi2[c]*(M_PI/180.0));
	}
	out<<s1<<" "<<aa1<<" "<<" "<<count1<<" "
	   <<s2<<" "<<aa2<<" "<<count2<<" ";	
	EnergyData<double> e = ff.energy(coor); 
	out<<e["energy"]<<" "<<e["vdw"]<<" "<<e["elec"]<<" "<<e["solv"]<<"\n";
// 	for(int c = 0; c < r2.size();c++){
// 	  cout<<chi2[c]<<" ";
// 	}
	
// 	cout<<"\n";
	++count2;
	++i2;
      }
    }


    if(r1.size()&&r2.size()==0){
      while(i1!=e1){
	//set the first
	double *chi1 =  i1->second.chi;
	for(int c = 0; c < r1.size();c++){
	  br.set_dihedral(coor,
			  r1[c][0],
			  r1[c][1],
			  r1[c][2],
			  r1[c][3],
			  chi1[c]*(M_PI/180.0));
	}
	out<<s1<<" "<<aa1<<" "<<" "<<count1<<" "
	   <<s2<<" "<<aa2<<" "<<count2<<" ";
	EnergyData<double> e = ff.energy(coor); 
	out<<e["energy"]<<" "<<e["vdw"]<<" "<<e["elec"]<<" "<<e["solv"]<<"\n";
	++i1;
	++count1;
      }
    }


    if(r1.size()==0&&r2.size()==0){
      out<<s1<<" "<<aa1<<" "<<" "<<count1<<" "
	 <<s2<<" "<<aa2<<" "<<count2<<" ";
      EnergyData<double> e = ff.energy(coor); 
      out<<e["energy"]<<" "<<e["vdw"]<<" "<<e["elec"]<<" "<<e["solv"]<<"\n";
    }

  }

  void mkenergy_tables(const Molecules & molecules,
		       const EnergyOptions<double> options,
		       MDB & mdb,
		       DunbrackLib & db,
		       string outfile){
    
    //1 Create fake polychains
    vector<string> amino = aminoacids();
    
    map<string,Molecules> mtemplate;

    for(int i=0;i<amino.size();i++){
      if(amino[i]=="PRO") continue;
      for(int p = 0;p<molecules.protein_size();p++){
	stringstream str;
	str<<"T"<<p+1;
	Protein pr(str.str());
	string aa = amino[i];
	vector<string> seq;
	for(int a=0;a<molecules.protein(p).fragment_size();a++){
	  stringstream str;
	  str<<aa<<"_"<<a+1;
	  seq.push_back(str.str());
	}
	pr.build_template(seq,mdb,"DEFAULT","DEFAULT");
	assign_backbone(pr,molecules.protein(p));
	pr.generate_missing(mdb,"DEFAULT","DEFAULT");
	mtemplate[aa].add_protein(pr);
      }
    }

    ofstream out;
    out.open(outfile.c_str());
    //self
    {
      map<string,Molecules>::iterator i1 = mtemplate.begin(),
	e1 = mtemplate.end();
      while(i1!=e1){
	int size = i1->second.fragment_size();
	for(int aa1=0;aa1<size;aa1++)
	  self_energy(i1->second,aa1,options,db,mdb,out);
	++i1;
      }
    }
    out<<"EndSelf\n";
    //inter
    {
      map<string,Molecules>::iterator i1 = mtemplate.begin(),
	i2 =  mtemplate.begin(), e1 = mtemplate.end(), e2 = mtemplate.end();

      while(i1!=e1){
	i2 =  mtemplate.begin();
	while(i2!=e2){
	  int size = i1->second.fragment_size();
	  for(int aa1=0;aa1<size;aa1++){
	    for(int aa2 = aa1+2; aa2<size;aa2++){
	      //	    cout<<"DEBUG iter ene "<<aa1<<" "<<aa2<<endl;
	      inter_energy(i1->second,i2->second,aa1,aa2,options,db,mdb,out);
	    }
	  }
	  ++i2;
	}
	++i1;
      }
    }
  }


  class MFTSolve {
    //Data
    typedef map<string, double>  SiteDouble;
    typedef map<string, vector<double> > SiteVector;
    typedef map<string, vector<int> > SiteIVector;
    typedef map<string, map<string,vector<vector<double> > > > SiteMatrix;

    vector<SiteVector>           W;
    //    vector<SiteIVector>          Size;
    vector<SiteVector>           E;
    vector<SiteVector>           Eref;
    vector<vector< SiteMatrix> > Eint;

  public:
    MFTSolve(){}

    void read_tables(string file){ 
      ifstream in;
      in.open(file.c_str());
      istream_iterator<string> iter(in),end;
      while(iter!=end){
	string tok = *iter; ++iter;
	if(tok=="EndSelf") break;
	string aa = *iter; ++iter;
	int pos = atoi(iter->c_str()); ++iter;
	if(pos>=W.size()){
	  W.resize(pos+1); E.resize(pos+1);
	}
	int size = atoi(iter->c_str()); ++iter;
	W[pos][aa].resize(size);
	E[pos][aa].resize(size);
	for(int i=0;i<size;i++){
	  double e = 0;
	  double tmp;
	  tmp = atof(iter->c_str()); ++iter;
	  tmp = atof(iter->c_str()); ++iter;
	  e += tmp;
	  tmp = atof(iter->c_str()); ++iter;
	  e += tmp;
	  tmp = atof(iter->c_str()); ++iter;
	  e += tmp;
	  E[pos][aa][i] = e;
	}
      }
      cout<<"Site size "<<W.size()<<" "<<E.size()<<endl;
      Eint.resize(W.size());

      for(int i=0;i<W.size();i++){
	Eint[i].resize(W.size());
	for(int j=i+2;j<W.size();j++){
	  //	  Eint[i].resize(Eint[i].size()+1);
	  SiteVector::iterator i1,i2,e1,e2;
	  i1=W[i].begin();
	  e1=W[i].end();
	  i2=W[j].begin();
	  e2=W[j].end();
	  while(i1!=e1){
	    string aa1 = i1->first;
	    i2=W[j].begin();
	    while(i2!=e2){
	      string aa2 = i2->first;
// 	      cout<<"DEBUG res Eint "<<i<<" ("<<Eint.size()<<") "<< aa1<<" "
// 		  <<j<<" ("<<Eint[i].size()<<") "<< aa2<<endl;
	      Eint[i][j][aa1][aa2].resize(W[i][aa1].size());
	      for(int qq=0;qq<W[i][aa1].size();qq++){
		Eint[i][j][aa1][aa2][qq].resize(W[j][aa2].size());
	      }
	      ++i2;
	    }
	    ++i1;
	  }
	}
      }


      //Now iter
     while(iter!=end){
	string aa1  = *iter; ++iter;
	int    pos1 = atoi(iter->c_str()); ++iter;
	int    rot1 = atoi(iter->c_str()); ++iter;
	string aa2  = *iter; ++iter;
	int    pos2 = atoi(iter->c_str()); ++iter;
	int    rot2 = atoi(iter->c_str()); ++iter;
	double e = 0;
	double tmp;
	tmp = atof(iter->c_str()); ++iter;
	tmp = atof(iter->c_str()); ++iter;
	e += tmp;
	tmp = atof(iter->c_str()); ++iter;
	e += tmp;
	tmp = atof(iter->c_str()); ++iter;
	e += tmp;
// 	cout<<"DEB "<<pos1<<" "<<pos2 <<" "
// 	    <<aa1<<" "<<aa2<<" "<<rot1<<" "<<rot2<<endl;
// 	cout<<"DEB "<<	Eint[pos1][pos2][aa1][aa2].size()<<" "<<W[pos1][aa1].size()<<" "
// 	    <<Eint[pos1][pos2][aa1][aa2][rot1].size()<<" "<<W[pos2][aa2].size()<<" "<<endl;
	if(rot1>=Eint[pos1][pos2][aa1][aa2].size()){
	  cout<<pos1<<" "<<pos2 <<" "
	      <<aa1<<" "<<aa2<<" "<<rot1<<" "<<rot2<<" skipped "<<endl;
	  cout<<"Size " <<Eint[pos1][pos2][aa1][aa2].size()<<endl;
	}
	else if(rot2>=Eint[pos1][pos2][aa1][aa2][rot1].size()){
	  cout<<pos1<<" "<<pos2 <<" "
	      <<aa1<<" "<<aa2<<" "<<rot1<<" "<<rot2<<" skipped2 "<<endl;;
	  cout<<"Size " <<Eint[pos1][pos2][aa1][aa2][rot1].size()<<endl;
	}
	Eint[pos1][pos2][aa1][aa2][rot1][rot2] = e;
// 	  min(e 
// 	  - E[pos1][aa1][rot1]
// 	      - E[pos2][aa2][rot2],0.0) ;
// 	double q = e 
// 	  - E[pos1][aa1][rot1]
// 	  - E[pos2][aa2][rot2];
// 	if(q>0) cout<<q<<endl;
// 	  min(e 
// 	      - E[pos1][aa1][rot1]
// 	      - E[pos2][aa2][rot2],
// 	      0.0);
// 	;
     }
    }

    void print_amino_frq(){
      cout<<"\n\n";
      vector<map<double,string> > aa_prob;
      aa_prob.resize(W.size());
      for(int i=0;i<W.size();i++){
	SiteVector::iterator iter = W[i].begin(),
	  end = W[i].end();
	double WW = 0;
	cout<<i+1<<" ";
	while(iter!=end){
	  string aa = iter->first;
	  double w = 0;
	  for(int qq=0;qq<iter->second.size();qq++){
	    w += iter->second[qq];
	  }
	  

	  cout<<aa<<" ";//<<w<<" ";
	  printf("%4.4f ",w);
	  WW += w;
	  aa_prob[i][-log(w)] = aa;
	  ++iter;
	}
	cout<<"("<<WW<<")"<<endl;
      }

      cout<<"\n\n>\n";
      for(int i=0;i<aa_prob.size();i++){
	cout<<three2one(aa_prob[i].begin()->second);
      }
      cout<<"\n";

//       map<double, string> seqs;
//       for(int i=0;i<
    }

    void solve(double T){
      //init W
      init_W();
      double error;
      int step = 0;
      double beta = 1/T;
      //
      {
	Eref = E;
	for(int i=0;i<Eref.size();i++){
	  SiteVector::iterator iter = Eref[i].begin(),
	    end = Eref[i].end();
	  while(iter!=end){
	    double qref = (double)iter->second.size();
	    qref = T*log(qref);
	    for(int q=0;q<iter->second.size();q++)
	      iter->second[q] = qref;
	    ++iter;
	  }
	}
      }
      error = iterate_W(beta);
      ++step;
      while(error> 0.000001){
	error = iterate_W(beta);
	++step;
	aout<<"Step "<<step<<" "<<error<<"\n";
      }
      print_amino_frq();
    }

  private:
    void init_W(){
      int size = W.size();
      for(int i=0;i<size;i++){
	SiteVector::iterator iter = W[i].begin(),
	  end = W[i].end();
	while(iter!=end){
	  int s = iter->second.size();
	  double w = 1.0/(double)s;
	  for(int q=0;q<s;q++) iter->second[q] = w;
	  ++iter;
	}
      }
    }

    double iterate_W(double beta){
      int size = W.size();
      double error =0;
      for(int i=0;i<size;i++){
	error += iterate_Site(i,beta);
      }
      return error;
    }

    double iterate_Site(int site, double beta){
      SiteVector w;
      SiteVector & w_old = W[site];
      SiteVector ene;
      double logZ =0;;
      init_site_vector(w,site);
      compute_Ene(ene,site);

      //
      {
	SiteVector::iterator iter = ene.begin(), 
	  end = ene.end();
	while(iter!=end){
	  vector<double> & e = iter->second;
	  for(int i=0;i<e.size();i++){
	    //	    cout<<"DEBUG energy "<<e[i]<<" "<< logZ<<" "<<endl;
	    if(logZ == 0) logZ = -beta*e[i];
	    else
	      logZ = logZ +log1p(exp(-beta*e[i]-logZ));
	  }
	  ++iter;
	}
      }
      //      cout<<"DEBUG Z "<<site<<" "<<logZ<<endl;
      //
      {
	double error = 0;
 	SiteVector::iterator iter = w.begin(), 
	  end =w.end();
	while(iter!=end){ 
	  int s = iter->second.size();
	  for(int q=0;q<s;q++){
	    double logw = -beta*ene[iter->first][q]-logZ;
	    iter->second[q]= exp(logw);
	    error += (iter->second[q]-w_old[iter->first][q])*
	      (iter->second[q]-w_old[iter->first][q]);
// 	    cout<<"DEBUG err "<<error<<" "<<iter->second[q]<<" "
// 		<<w_old[iter->first][q]<<
//	      endl;
	  }
	  ++iter;
	}
	W[site] = w;
	return error;
      }
    }

    void compute_Ene( SiteVector & ene, int site){
      init_site_vector(ene,site);

      int size = W.size();
      SiteVector::iterator iter = ene.begin(),
	end = ene.end();
      while(iter!=end){
	string aa1 = iter->first;
	int rot_size1 = iter->second.size();
	for(int r1 = 0;r1<rot_size1;r1++){
	  double energy=0;
	  energy  = Eref[site][aa1][r1];//E[site][aa1][r1]; //self
	  for(int qq = 0; qq < size; qq++){
	    if((qq==site)||(qq+1==site)||(qq-1==site)) continue;
	    int idx1 = min(qq,site);
	    int idx2 = max(qq,site);

	    SiteVector::iterator i2 = W[qq].begin(),
	      e2 = W[qq].end();
	    while(i2!=e2){
	      string aa2 = i2->first;
	      int rot_size2 = i2->second.size();
	      for(int r2=0;r2<rot_size2;r2++){
		if(qq>site)
		  energy += W[qq][aa2][r2]*Eint[idx1][idx2][aa1][aa2][r1][r2];
		else
		  energy += W[qq][aa2][r2]*Eint[idx1][idx2][aa2][aa1][r2][r1];
	      }
	      ++i2;
	    }
	  }
	  ene[aa1][r1] = energy;
	}
	++iter;
      }
    }

    void init_site_vector(SiteVector & ene,int site){
      // ene = W[site];
      SiteVector::iterator iter = W[site].begin(),
	end = W[site].end();
      while(iter!=end){
	ene[iter->first].resize(iter->second.size());
	for(int i=0;i<iter->second.size();i++){
	 ene[iter->first][i] = 0;
	}
	++iter;
      }
    }
  };

#endif
  class MFTResSolve {
    typedef map<string,double> SiteVector;
    typedef map<string,map<string,double> > SiteMatrix;

    vector<SiteVector> W;
    vector<SiteVector> E;
    vector<vector<SiteMatrix> > Eint;
  public:
    MFTResSolve(string tbl){
      read_tables(tbl);
    }
    void print_amino_frq(){
      vector<map<double,string> > aa_prob;
      aa_prob.resize(W.size());
      for(int i=0;i<W.size();i++){
	SiteVector::iterator iter= W[i].begin(),
	  end = W[i].end();
	while(iter!=end){
	  aa_prob[i][-log(iter->second)] = iter->first;
	  ++iter;
	}
      }

      cout<<"\n\n";
      for(int i=0;i<W.size();i++){
	SiteVector::iterator iter= W[i].begin(),
	  end = W[i].end();
	double WW = 0;
	printf("%4i ",i+1);
	while(iter!=end){	  
	  cout<<iter->first<<" ";
	  printf("%4.4f ",iter->second);
	  WW += iter->second;
	  ++iter;
	}
	cout<<"("<<WW<<")"<<endl;
      }

      cout<<"\n\n>\n";
      for(int i=0;i<aa_prob.size();i++){
	cout<<three2one(aa_prob[i].begin()->second);
      }
      cout<<"\n";
    }
    void solve(double T){
      //init W
      init_W();
      double error;
      int step = 0;
      double beta = 1/T;
      error = iterate_W(beta);
      ++step;
      while(error> 0.000001){
	error = iterate_W(beta);
	++step;
	aout<<"Step "<<step<<" "<<error<<"\n";
      }
      print_amino_frq();
    }

    void mc_solve_sub(double T, int steps, vector<string> init,
		      vector<int> pos){
      set<int> mut;
      for(int i=0;i<pos.size();i++)
	if(pos[i]==1) mut.insert(i);
      
      vector<string> seq = init;
      vector<string> seq_buff = init;
      steps = steps/10;

      double e = mc_energy(seq);
      for(int qq=10;qq>=1;qq--){
	for(int i=0;i<steps;i++){
	  mc_move(seq,mut);
	  double ee = mc_energy(seq);
	  double r = drand48();
	  double q = ee - e;
	  bool acc = false;
	  if(r<exp(-q/(qq*T))) acc =1;
	  if(acc){
	    e = ee;
	    seq_buff = seq;
	  }
	  else seq = seq_buff;
	  
	  if(i%100==0){
	    cout<<"Step "<<qq<<" "<<i+1<<" "<<e<<"\n";
	    for(int z=0;z<seq.size();z++){
	      cout<<three2one(seq[z]);
	    }
	    cout<<"\n";
	  }
	}
      }
    }


    void mc_solve(double T, int steps, vector<string> nat){
      double enat = mc_energy(nat);
      vector<string> seq;
      vector<string> seq_buff;
      steps = steps/10;
      seq.resize(E.size());
      for(int i=0;i<seq.size();i++){
	seq[i] = "ALA";
      }
      seq_buff = seq;

      double e = mc_energy(seq);
      e = e + (e-enat)*(e-enat) -enat;
      for(int qq=10;qq>=1;qq--){
	for(int i=0;i<steps;i++){
	  mc_move(seq);
	  double ee = mc_energy(seq);
	  ee = ee + (ee-enat)*(ee-enat)-enat;
	  double r = drand48();
	  double q = ee - e;
	  bool acc = false;
	  if(r<exp(-q/(qq*T))) acc =1;
	  if(acc){
	    e = ee;
	    seq_buff = seq;
	  }
	  else seq = seq_buff;
	  
	  if(i%100==0){
	    cout<<"Step "<<qq<<" "<<i+1<<" "<<e<<"\n";
	    for(int z=0;z<seq.size();z++){
	      cout<<three2one(seq[z]);
	    }
	    cout<<"\n";
	  }
	}
      }
    }


    double mc_energy(vector<string> & seq){
      double e = 0;
      for(int i=0;i<seq.size();i++) e += E[i][seq[i]];
      for(int i=0;i<seq.size();i++){
	for(int j=i+1;j<seq.size();j++){
	  e += Eint[i][j][seq[i]][seq[j]];
	}
      }

//       for(int i=0;i<seq.size()-2;i++){
// 	if(seq[i]==seq[i+1]&&seq[i]==seq[i+2]) e = e+5;
//       }
      return e;
    }

    double random(int q, double &sigma){
      static vector<string> qq = aminoacids();
      vector<string> seq;
      double e = 0;
      double e2;
      seq.resize(E.size());
      for(int i=0;i<q;i++){
	for(int q=0;q<seq.size();q++)
	  seq[q] = qq[lrand48()%qq.size()];
	double ee = mc_energy(seq);
	e += ee;
	e2 += ee*ee;
      }
      sigma = e2/q-(e/q)*(e/q);
      return e/q;

    }

  private:
    void mc_move(vector<string> & seq){
      static vector<string> qq = aminoacids();
      for(int ww=0;ww<2;ww++){
	int pos = lrand48()%seq.size();
	int a = lrand48()%20;
	seq[pos] = qq[a];
      }
    }


    void mc_move(vector<string> & seq,set<int> & pos){
      static vector<string> qq = aminoacids();
      int p = lrand48()%pos.size();
      int q = lrand48()%20;
      set<int>::iterator iter = pos.begin();
      int count = 0;
      while(count<p){ ++iter;++count;}
      seq[*iter] = qq[q];
    }
    void read_tables(string file){ 
      ifstream in;
      in.open(file.c_str());
      istream_iterator<string> iter(in),end;
      //self 
      //Site aa nsite double 
      while(iter!=end){
	string tok = *iter; ++iter;
	if(tok=="EndSelf") break;
	string aa = *iter; ++iter;
	int pos = atoi(iter->c_str()); ++iter;
	if(pos>=W.size()){
	  W.resize(pos+1); E.resize(pos+1);
	}

	double ene = atof(iter->c_str()); ++iter;
	E[pos][aa] = ene;
      }

      Eint.resize(W.size());
      for(int i=0;i<Eint.size();i++){
	Eint[i].resize(W.size());
      }

      
      //Pair
      //pos1 aa1 pos2 aa2 double
      while(iter!=end){
	int pos1 = atoi(iter->c_str());++iter;
	string aa1 = *iter; ++iter;
	int pos2 = atoi(iter->c_str());++iter;
	string aa2 = *iter; ++iter;
	double e  = atof(iter->c_str());++iter;
	Eint[pos1][pos2][aa1][aa2] = e;
	Eint[pos2][pos1][aa2][aa1] = e;
      }
    }
    void init_W(){
      W = E;
      int size = W.size();
      for(int i=0;i<size;i++){
	double q = (double)W[i].size();
	q = 1/q;
	SiteVector::iterator iter = W[i].begin(),
	  end = W[i].end();
	while(iter!=end){
	  iter->second = q;
	  ++iter;
	}
      }
    }

    double iterate_W(double beta){
      int size = W.size();
      double error =0;
      for(int i=0;i<size;i++){
	error += iterate_Site(i,beta);
      }
      return error;
    }

    double iterate_Site(int site, double beta){
      SiteVector w;
      SiteVector & w_old = W[site];
      SiteVector ene;
      double logZ =0;;
      init_site_vector(w,site);
      compute_Ene(ene,site);

      //
      {
	SiteVector::iterator iter = ene.begin(), 
	  end = ene.end();
	while(iter!=end){
	  if(logZ == 0) logZ = -beta*iter->second;
	  else
	    logZ = logZ +log1p(exp(-beta*iter->second-logZ));

	  //	  cout<<"DEBUG energy "<<iter->second<<" "<< logZ<<" "<<endl;
	  ++iter;
	}
      }
      //      cout<<"DEBUG Z "<<site<<" "<<logZ<<endl;
      //
      {
	double error = 0;
 	SiteVector::iterator iter = w.begin(), 
	  end =w.end();
	while(iter!=end){ 
	  double logw = -beta*ene[iter->first]-logZ;
	  iter->second= exp(logw);
	  error += (iter->second-w_old[iter->first])*
	    (iter->second-w_old[iter->first]);
	  ++iter;
	}
	W[site] = w;
	return error;
      }
    }
    
    void compute_Ene( SiteVector & ene, int site){
      init_site_vector(ene,site);
      int size = W.size();
      SiteVector::iterator iter = ene.begin(),
	end = ene.end();
      
      while(iter!=end){
	string aa1 = iter->first;
	double energy = 0;
	double eneint = 0;
	energy = E[site][aa1]; //self energy
	for(int qq = 0; qq < size; qq++){
	  if(qq==site) continue;
	  SiteVector::iterator i2 = W[qq].begin(),
	    e2 = W[qq].end();
	  while(i2!=e2){
	    string aa2 = i2->first;
	    eneint += W[qq][aa2]*Eint[site][qq][aa1][aa2];
	    ++i2;
	  }
	}
	//	cout<<"DEBUG "<<site<<" "<<aa1<<" "<<energy<<" "<<eneint<<"\n";
	ene[aa1] = energy+10*eneint;
	++iter;
      }
    }
    
    void init_site_vector(SiteVector & ene,int site){
      // ene = W[site];
      SiteVector::iterator iter = E[site].begin(),
	end = E[site].end();
      while(iter!=end){
	ene[iter->first] = 0;
	++iter;
      }
    }    
  };


  void read_res_map(map<string, double> & aa,istream_iterator<string> & iter){

    while(1){
      string s1 = *iter; ++iter;
      if(s1=="BEGIN"){// cout<<"BEGIN"<<endl; 
	continue;
      }
      if(s1=="END"){ //cout<<"DONE"<<endl; 
	return;
      }
      //      string s2 = *iter; ++iter;
      double p = atof(iter->c_str()); ++iter;
      aa[s1] = p;
    }
  }

  void read_pair(map<string, map<string,double> > & aa, istream_iterator<string> & iter){
    //    cout<<"DEBUG READ_PAIR"<<endl;
    while(1){
      string s1 = *iter; ++iter;
      //      cout<<"DEBUG "<<s1<<endl;
      if(s1=="BEGIN"){ //cout<<"BEGIN"<<endl; 
	continue;
      }
      if(s1=="END"){ cout<<"DONE"<<endl; return;}
      string s2 = *iter; ++iter;
      double p = atof(iter->c_str()); ++iter;
      aa[s1][s2] = p;
    }
  }

  void to3(string &s){
    for(int i=0;i<s.size();i++){
      char c=s[i];
      if(c=='H'||c=='G') c='H';
      else if(c=='E'||c=='B') c='E';
      else c = 'C';
      s[i] = c;
    }

  }


  void read_cont_stat(map<int,map<string,double> > & cstat,
		      istream_iterator<string> & iter){

    int n = -1;
    while(1){
      string s1 = *iter;++iter;
      if(s1=="BEGIN"){
	n = atoi(iter->c_str());++iter;
	cout<<"DEBUG BEGIN "<<n<<endl;
      }
      else if(s1=="END"){
	n = atoi(iter->c_str());++iter;
	cout<<"DEBUG END "<<n<<endl;
	if(n==19) return;
      }
      else {
	double p = atof(iter->c_str());++iter;
	cstat[n][s1] = p;
      }
    }
  }
		      
  void count_n(const Molecules & m, vector<int> &n){
    vector<string> seq = m.sequence(2);
    vector<int> CA;
    for(int i=0;i<m.atom_size();i++){
      if(m[i].name()=="CA") CA.push_back(i);
    }
    if(CA.size()!=m.fragment_size()) return;
    n.resize(CA.size());
    
    for(int i=0;i<n.size();i++){
      for(int j=i+2;j<n.size();j++){
	if(m.distance(CA[i],CA[j])<7.5){
	  ++n[i] ;
	  ++n[j];
	}
      }
    }
    for(int i=0;i<n.size();i++){
      cout<<"CONT "<<seq[i]<<" "<<n[i]<<"\n";
      if(n[i]>=20) n[i]=19;
    }
  }

  void mkres_tables(const Molecules  & m, string tbl_file, string file){
    ifstream tbl; tbl.open(tbl_file.c_str());
    istream_iterator<string> iter(tbl),end;  

    map<int,map<string,double> > cstat;

    map<string,double> aa;
    map<string,double> aa_H;
    map<string,double> aa_E;
    map<string,double> aa_C;

    map<string,map<string,double> > AA;

    map<string,map<string,double> > d_1;
    map<string,map<string,double> > d_2;
    map<string,map<string,double> > d_3;
    map<string,map<string,double> > d_4;

    map<string,map<string,double> > d_1HH;
    map<string,map<string,double> > d_2HH;
    map<string,map<string,double> > d_3HH;
    map<string,map<string,double> > d_4HH;
    
    map<string,map<string,double> > d_1EE;
    map<string,map<string,double> > d_2EE;
    map<string,map<string,double> > d_3EE;
    map<string,map<string,double> > d_4EE;

    map<string,map<string,double> > d_1CC;
    map<string,map<string,double> > d_2CC;
    map<string,map<string,double> > d_3CC;
    map<string,map<string,double> > d_4CC;
    
    map<string,map<string,double> > d_1HC;
    map<string,map<string,double> > d_2HC;
    map<string,map<string,double> > d_3HC;
    map<string,map<string,double> > d_4HC;
    
    map<string,map<string,double> > d_1CH;
    map<string,map<string,double> > d_2CH;
    map<string,map<string,double> > d_3CH;
    map<string,map<string,double> > d_4CH;


    map<string,map<string,double> > d_1EC;
    map<string,map<string,double> > d_2EC;
    map<string,map<string,double> > d_3EC;
    map<string,map<string,double> > d_4EC;

    map<string,map<string,double> > d_1CE;
    map<string,map<string,double> > d_2CE;
    map<string,map<string,double> > d_3CE;
    map<string,map<string,double> > d_4CE;


    map<string,map<string,double> > d_1HE;
    map<string,map<string,double> > d_2HE;
    map<string,map<string,double> > d_3HE;
    map<string,map<string,double> > d_4HE;

    map<string,map<string,double> > d_1EH;
    map<string,map<string,double> > d_2EH;
    map<string,map<string,double> > d_3EH;
    map<string,map<string,double> > d_4EH;
    
    read_cont_stat(cstat,iter);

    cout<<"AA"<<" ";
    read_res_map(aa,iter);
    cout<<"AA_H"<<" ";
    read_res_map(aa_H,iter);
    cout<<"AA_E"<<" ";
    read_res_map(aa_E,iter);
    cout<<"AA_C"<<" ";
    read_res_map(aa_C,iter);

    cout<<"AA"<<" ";
    read_pair(AA,iter);

    cout<<"d1"<<" ";
    read_pair(d_1,iter);
    cout<<"d2"<<" ";
    read_pair(d_2,iter);
    cout<<"d3"<<" ";
    read_pair(d_3,iter);
    cout<<"d4"<<" ";
    read_pair(d_4,iter);

    cout<<"d1_HH"<<" ";
    read_pair(d_1HH,iter);
    cout<<"d1_EE"<<" ";
    read_pair(d_1EE,iter);
    cout<<"d1_CC"<<" ";
    read_pair(d_1CC,iter);

    cout<<"d1_CH"<<" ";
    read_pair(d_1CH,iter);
    cout<<"d1_HC"<<" ";
    read_pair(d_1HC,iter);

    cout<<"d1_CE"<<" ";
    read_pair(d_1CE,iter);
    cout<<"d1_EC"<<" ";
    read_pair(d_1EC,iter);

    cout<<"d1_HE"<<" ";
    read_pair(d_1HE,iter);
    cout<<"d1_EH"<<" ";
    read_pair(d_1EH,iter);

    cout<<"d2_HH"<<" ";
    read_pair(d_2HH,iter);
    cout<<"d2_EE"<<" ";
    read_pair(d_2EE,iter);
    cout<<"d2_CC"<<" ";
    read_pair(d_2CC,iter);

    cout<<"d2_CH"<<" ";
    read_pair(d_2CH,iter);
    cout<<"d2_HC"<<" ";
    read_pair(d_2HC,iter);

    cout<<"d2_CE"<<" ";
    read_pair(d_2CE,iter);
    cout<<"d2_EC"<<" ";
    read_pair(d_2EC,iter);

    cout<<"d2_HE"<<" ";
    read_pair(d_2HE,iter);
    cout<<"d2_EH"<<" ";
    read_pair(d_2EH,iter);

    cout<<"d3_HH"<<" ";
    read_pair(d_3HH,iter);
    cout<<"d3_EE"<<" ";
    read_pair(d_3EE,iter);
    cout<<"d3_CC"<<" ";
    read_pair(d_3CC,iter);

    cout<<"d3_CH"<<" ";
    read_pair(d_3CH,iter);
    cout<<"d3_HC"<<" ";
    read_pair(d_3HC,iter);

    cout<<"d3_CE"<<" ";
    read_pair(d_3CE,iter);
    cout<<"d3_EC"<<" ";
    read_pair(d_3EC,iter);

    cout<<"d3_HE"<<" ";
    read_pair(d_3HE,iter);
    cout<<"d3_EH"<<" ";
    read_pair(d_3EH,iter);

    cout<<"d4_HH"<<" ";
    read_pair(d_4HH,iter);
    cout<<"d4_EE"<<" ";
    read_pair(d_4EE,iter);
    cout<<"d4_CC"<<" ";
    read_pair(d_4CC,iter);

    cout<<"d4_CH"<<" ";
    read_pair(d_4CH,iter);
    cout<<"d4_HC"<<" ";
    read_pair(d_4HC,iter);

    cout<<"d4_CE"<<" ";
    read_pair(d_4CE,iter);
    cout<<"d4_EC"<<" ";
    read_pair(d_4EC,iter);

    cout<<"d4_HE"<<" ";
    read_pair(d_4HE,iter);
    cout<<"d4_EH"<<" ";
    read_pair(d_4EH,iter);


    //Comp sec
    AlmSecStruct ss(m);
    Coor<double> c(m);
    string s = ss.secstruct(c);
    to3(s);
    ofstream out;
    out.open(file.c_str());
    //Self
    vector<int> N;
    count_n(m,N);
    int n = s.size();
    vector<string> ax = aminoacids();
    for(int i=0;i<n;i++){
      for(int q=0;q<ax.size();q++){
	string a1 = ax[q];
	double e = 0;
	double Pa = aa[a1];
	//e = -log(Pa);

	//SS bias
	double Pss;
	if(s[i]=='H')
	  Pss = aa_H[a1];
	else if(s[i]=='E')
	  Pss = aa_E[a1];
	else if(s[i]=='C')
	  Pss = aa_C[a1];
	
	e = e - log(Pss/Pa)-log(cstat[N[i]][a1]);
	//e = -log(cstat[N[i]][a1]);
	//	e= 0;
	//	cout<<"Site "<<a1<<" "<<i<<" "<<e<<" "<<s[i]<<"\n";
	out<<"Site "<<a1<<" "<<i<<" "<<e<<"\n";
      }
    }
    out<<"EndSelf\n";


    for(int i=0;i<n;i++){
      int ca1;
      vector<int> fra = m.fragment_atoms(i);
      for(int qq1=0;qq1<fra.size();qq1++)
	if(m[fra[qq1]].name()=="CA")
	  ca1 = fra[qq1];
      for(int q1=0;q1<ax.size();q1++){
	string a1 = ax[q1];
	for(int j=i+1;j<n;j++){
	  int ca2;
	  vector<int> fra2 = m.fragment_atoms(j);
	  for(int qq1=0;qq1<fra2.size();qq1++)
	    if(m[fra2[qq1]].name()=="CA")
	      ca2 = fra2[qq1];

	  double d = m.distance(ca1,ca2);
	  for(int q2=0;q2<ax.size();q2++){
	    string a2 = ax[q2];
	    if(abs(i-j)<5){
	      double e = 0;
	      if(j==i+1)
		e = -log(AA[a1][a2]/(aa[a1]*aa[a2]));
	      //	      cout<<"DEBUG 1"<<endl;
	      out<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<e<<"\n";
	    }
	    else {
	      char ss1 = s[i];
	      char ss2 = s[j];

	      double P;
	      if(ss1=='H'&&ss2=='H'){
		if(d<5){
		  P = d_1HH[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2HH[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2HH[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2HH[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      else if(ss1=='E'&&ss2=='E'){
		if(d<5){
		  P = d_1EE[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2EE[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2EE[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2EE[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      else if(ss1=='C'&&ss2=='C'){
		if(d<5){
		  P = d_1CC[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2CC[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2CC[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2CC[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      else if(ss1=='H'&&ss2=='C'){
		if(d<5){
		  P = d_1HC[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2HC[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2HC[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2HC[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      if(ss1=='C'&&ss2=='H'){
		if(d<5){
		  P = d_1CH[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2CH[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2CH[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2CH[a1][a2];
		}
		else {
		  P=1;
		}
	      }

	      else if(ss1=='E'&&ss2=='C'){
		if(d<5){
		  P = d_1EC[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2EC[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2EC[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2EC[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      if(ss1=='C'&&ss2=='E'){
		if(d<5){
		  P = d_1CE[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2CE[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2CE[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2CE[a1][a2];
		}
		else {
		  P=1;
		}
	      }

	      else if(ss1=='E'&&ss2=='H'){
		if(d<5){
		  P = d_1EH[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2EH[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2EH[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2EH[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      if(ss1=='H'&&ss2=='E'){
		if(d<5){
		  P = d_1HE[a1][a2];
		}
		else if(d<7.5&&d>=5){
		  P = d_2HE[a1][a2];
		}
		else if(d<10&&d>=7.5){
		  P = d_2HE[a1][a2];
		}
		else if(d<12&&d>=10){
		  P = d_2HE[a1][a2];
		}
		else {
		  P=1;
		}
	      }
	      double Pss1,Pss2;
	      if(ss1=='H')
		Pss1 = aa_H[a1];
	      else if(ss1=='E')
		Pss1 = aa_E[a1];
	      else if(ss1=='C')
		Pss1 = aa_C[a1];
	      if(ss2=='H')
		Pss2 = aa_H[a2];
	      else if(ss2=='E')
		Pss2 = aa_E[a2];
	      else if(ss2=='C')
		Pss2 = aa_C[a2];
// 	      if(d<12)
// 		cout<<"DEBUG 2"<<" "<<d<<" "<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<P
// 		    <<" "<<ss1<<" "<<ss2<<endl;
	      if(P==0){
		double e = 0;
		if(j==i+1)
		  e = -log(AA[a1][a2]/(aa[a1]*aa[a2]));
		cout<<"WARN "<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<100<<"\n";
		out<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<e<<"\n";
	      }
	      else {
		//		out<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<-log(P/(Pss1*Pss2))<<"\n";
		double e = 0;
		if(j==i+1)
		  e = -log(AA[a1][a2]/(aa[a1]*aa[a2]));
		//if(d<12&&P!=1)
		if(P==1)
		  out<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<-log(P)+e<<"\n";
		else
		  out<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<-log(P)+e<<"\n";
		  //		else 
		  //		  out<<i<<" "<<a1<<" "<<j<<" "<<a2<<" "<<0<<"\n";
	      }
	    }
	  }	
	}
      }
    }
    //Write enes
    
//     EnvConst<double> env(m,1);
//     EnvPairConst<double> envpair(m,1);
//     Coor<double> coor(m);
//     ofstream out;
//     out.open(file.c_str());
//     //Self
//     vector<string> aa = aminoacids();
//     int size = m.fragment_size();
//     for(int i=0;i<size;i++){
//       for(int q = 0;q<aa.size();q++){
// 	if(aa[q]=="CYS"||aa[q]=="PRO") continue;
// 	out<<"Site "<<" "<<aa[q]<<" "<<i<<" "<<10*0.23*env.site_energy(i,aa[q],coor)<<"\n";
//       }
//     }
//     out<<"EndSelf\n";

//     for(int i=0;i<size;i++){
//        for(int j=i+1;j<size;j++){
// 	 for(int q1= 0;q1<aa.size();q1++){
// 	   if(aa[q1]=="CYS"||aa[q1]=="PRO") continue;
// 	   for(int q2= 0;q2<aa.size();q2++){
// 	     if(aa[q2]=="CYS"||aa[q2]=="PRO") continue;
// 	     out<<i<<" "<<aa[q1]<<" "<<j<<" "<<aa[q2]<<" "
// 		 <<4*0.23*envpair.site_energy(i,j,aa[q1],aa[q2],coor)<<"\n";
// 	   } 
// 	 }
//        }
//     }
    
//   }
  }
}


#ifdef __MFT__SC__
template<>
inline string to_string(const mft::MFTSolve & pdb){
  return "<MFTSolve>";
}
#endif
template<>
inline string to_string(const mft::MFTResSolve & pdb){
  return "<MFTResSolve>";
}

extern "C" {
  void init_mft(){
    //declarations here

    Module mod = Module("mft");
#ifdef __MFT__SC__
    mod.def_function("mkenergy_tables","Creates mft energy tables",mft::mkenergy_tables);


    Class<mft::MFTSolve >(mod.self(),"solver")
      .def_method("solve",&mft::MFTSolve::solve)
      .def_method("read_tables",&mft::MFTSolve::read_tables);
#endif
    mod.def_function("mkres_tables","",mft::mkres_tables);
    Class<mft::MFTResSolve,ArgList1<string>   >(mod.self(),"ressolver")
      .def_method("solve",&mft::MFTResSolve::solve)
      .def_method("mc_solve",&mft::MFTResSolve::mc_solve)
      .def_method("mc_solve_sub",&mft::MFTResSolve::mc_solve_sub)
      .def_method("energy",&mft::MFTResSolve::mc_energy)
      .def_method("random",&mft::MFTResSolve::random)
      ;

  }
}