/***************************************************************************
                          roseff.cpp  -  description
                             -------------------
    begin                : Tue Nov 18 01:45:35 2003
    copyright            : (C) 1999-2006 by Cavalli Andrea
    author               : : cavalli $
    date                 : : 2003/12/09 11:10:33 $
    id                   : : coor.h,v 1.1.2.2 2003/12/09 11:10:33 cavalli Exp $
    email                : cavalli@bioc.unizh.ch amc82@cam.ac.uk
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   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 <almost.h>
#include <rose/roseff.h>
#include <molecules/molecules.h>
#include <pdb.h>
#include <algorithm/sec_packing_data.h>
namespace Almost {

  RoseFF::RoseFF(const Molecules & molecules){
    molecules_ = & molecules;
    for(int i=0;i<molecules.atom_size();i++){
      if(molecules[i].name()=="N") N.push_back(i);
      if(molecules[i].name()=="C") C.push_back(i);
      if(molecules[i].name()=="CA") CA.push_back(i);
      if((molecules[i].name()=="O")||
	 (molecules[i].name()=="OT1")) O.push_back(i);
      if(molecules[i].name()=="CB") CB.push_back(i);
    }
    W_ss = 0.1;
    w_ss= 0.07;
    w_hs= 0.16;
    w_hh =0.0;;

    apmf = NULL;
    bhb = NULL;
    hb = NULL;
    ff = NULL;
    rgyr = NULL;

    w_solv = 0.1;
    w_elec = w_vdw = 0.232566*0.2;

    w_env =0.1; w_pair = 2;
    idx.resize(molecules.fragment_size());
    frag_size = idx.size();
    for(int i=0;i<molecules.fragment_size();i++){
      string red = molecules.fragment_name(i);
      //string(seq[i].begin(),seq[i].begin()+3);
      idx[i] = one2index__(three2one(red));
      if(idx[i]>20){
	cout<<"Undefined amino-acid "<<red<<" "<<idx[i]<<endl;
	exit(1);
      }
    }
  }

  RoseFF::~RoseFF(){
    if(apmf!=NULL) delete apmf;
    if(bhb!=NULL) delete bhb; 
    if(hb!=NULL) delete hb; 
    if(ff!=NULL) delete ff;
    if(rgyr!=NULL) delete rgyr;
  }

  void RoseFF::init_apmf(string p1, int kind, double w){
    if(apmf) delete apmf;
    apmf = new APMF<double>(*molecules_,p1,kind,w);
  }

  void RoseFF::init_force_field(const EnergyOptions<double> & options,
				const MDB & mdb){
    if(ff) delete ff;
    ff = new ForceField<double>(*molecules_,options,mdb);

  }

  void RoseFF::init_bhb(string sec, double w, int mode){
    if(bhb) delete bhb;
    bhb = new BHB<double>(*molecules_,sec,w,mode);
  }

  void RoseFF::init_hb(double l){
    if(hb) delete hb;
    hb = new HBSecStruct<double>(*molecules_);
    hb->set_lambda(l);
  }

  void RoseFF::disable_hb(){
    if(hb) delete hb;
    hb = 0;
  }

  void RoseFF::init_rgyr(double l){
    if(rgyr) delete rgyr;
    rgyr = new RGYRConst(*molecules_,
			 NameSelection(*molecules_,"///"),
			 l,
			 2.83*exp(0.34*log((double)molecules_->fragment_size()))
			 );
  }


  void RoseFF::setW_ss(double d){ W_ss = d;}
  double RoseFF::getW_ss() const { return W_ss;}

  double RoseFF::ss_energy(const Coor<double> & coor, string ss){
    SecStruct secstruct(ss);
    double SS_energy_term = SS_energy(coor,secstruct);
    double HH_energy_term = 0;//HH_energy(coor,secstruct);
    double HS_energy_term = HS_energy(coor,secstruct);
    return W_ss*( w_ss*SS_energy_term + 
		  w_hh*HH_energy_term + 
		  w_hs*HS_energy_term );
  }

  double RoseFF::apmf_energy(const Coor<double> & coor){
    if(apmf)
      return apmf->energy(coor);
    else
      return 0;
  }

  double RoseFF::ff_energy(const Coor<double> & coor){
    if(ff){
      EnergyData<double> data = ff->energy(coor);
      return w_vdw*data["vdw"]+
	+w_elec*data["elec"] 
	+ w_solv*data["solv"]
	;
    }
    else
      return 0;
  }

  //   double RoseFF::cc_energy(const Coor<double> & coor){
  //     double cc = 0;
  //     for(int i=0;i<cc.size();i++){
  //       cc += cc.constraints()[i]->energy(coor);
  //     }

  //   }

  double RoseFF::bhb_energy(const Coor<double> & coor, string & sec){
    if(bhb)
      return bhb->energy2(coor,sec);
    else 
      return 0;
  }

  double RoseFF::hb_energy(const Coor<double> & coor, string & sec){
    if(hb)
      return hb->energy(coor,sec);
    else 
      return 0;
  }

  double RoseFF::rgyr_energy(const Coor<double> & coor){
    if(rgyr)
      return rgyr->energy2(coor);
    else 
      return 0;
  }

  double RoseFF::hh_energy(const Coor<double> & coor, string ss){
    SecStruct secstruct(ss);
    double HH_energy_term = HH_energy(coor,secstruct);
    return  W_ss*w_hh*HH_energy_term;
  }

  double RoseFF::he_energy(const Coor<double> & coor, string ss){
    SecStruct secstruct(ss);
    double HS_energy_term = HS_energy(coor,secstruct);
    return  W_ss*w_hs*HS_energy_term;
  }


  double RoseFF::ee_energy(const Coor<double> & coor, string ss){
    SecStruct secstruct(ss);
    double SS_energy_term = SS_energy(coor,secstruct);
    return  W_ss*w_ss*SS_energy_term;
  }

  

  double RoseFF::SS_energy(const Coor<double> & coor, const SecStruct & secstruct){
    //Get Beta >= 2
    vector<const SecStructElement * > beta;
    vector<int> dimer_strand;
    vector<int> dimer_aa;
    int ndimers = 0;
    for(int i=0;i<secstruct.ss_element.size();i++){
      if((secstruct.is_beta(i)) &&
	 (secstruct.ss_element[i].len>1)
	 ) {
	beta.push_back(&secstruct.ss_element[i]);
	int l = secstruct.ss_element[i].len-1;
	dimer_strand.resize(ndimers+l);
	dimer_aa.resize(ndimers+l);
	for(int k=ndimers;k<ndimers+l;k++){
	  dimer_strand[k] = i;
	  dimer_aa[k] = secstruct.ss_element[i].beg + (k-ndimers);
	}
	ndimers += l;
      }
    }

    vector<vector<double> > dimer_pair_score;

    dimer_pair_score.resize(ndimers);
    for(int i=0;i<ndimers;i++){
      dimer_pair_score[i].resize(ndimers);
    }

    
    for(int d1=0;d1<ndimers;d1++){
      //      cout<<"DEBUG "<<dimer_aa[d1]<<" "<<dimer_aa[d1]+1<<endl;
      double v1[3],co1[3],co2[3];
      int N1 = N[dimer_aa[d1]]*Coor<double>::DIM;
      int C1 = C[dimer_aa[d1]]*Coor<double>::DIM;
      int C2 = C[dimer_aa[d1]+1]*Coor<double>::DIM;
      int O1 = O[dimer_aa[d1]]*Coor<double>::DIM;
      int O2 = O[dimer_aa[d1]+1]*Coor<double>::DIM;
      //       cout<<"DEBUG "
      // 	  <<N[dimer_aa[d1]]+1<<" "
      // 	  <<C[dimer_aa[d1]]+1<<" "
      // 	  <<C[dimer_aa[d1]+1]+1<<" "
      // 	  <<O[dimer_aa[d1]]+1<<" "
      // 	  <<O[dimer_aa[d1]+1]+1<<endl;
      //       if(O2>=coor.size()*4){
      // 	cout<<O.size()<<" "<<C.size()<<" "<<N.size()<<endl;
      // 	exit(1);
      //       }

      for(int k=0;k<3;k++){
	v1[k] = coor.coor[C2+k]-coor.coor[N1+k];
	co1[k] = coor.coor[O1+k]-coor.coor[C1+k];
	co2[k] = coor.coor[O2+k]-coor.coor[C2+k];
      }
      norm(co1);norm(co2);
      for(int d2=d1+1;d2<ndimers;d2++){
	int loop;
	if(dimer_strand[d1]==dimer_strand[d2]) continue;
	else loop = secstruct.ss_element[dimer_strand[d2]].beg-
	  secstruct.ss_element[dimer_strand[d1]].end;
	double v2[3],co3[3],co4[3];
	int N2 = N[dimer_aa[d2]]*Coor<double>::DIM;
	C1 = C[dimer_aa[d2]]*Coor<double>::DIM;
	C2 = C[dimer_aa[d2]+1]*Coor<double>::DIM;
	O1 = O[dimer_aa[d2]]*Coor<double>::DIM;
	O2 = O[dimer_aa[d2]+1]*Coor<double>::DIM;
	for(int k=0;k<3;k++){
	  v2[k] = coor.coor[C2+k]-coor.coor[N2+k];
	  co3[k] = coor.coor[O1+k]-coor.coor[C1+k];
	  co4[k] = coor.coor[O2+k]-coor.coor[C2+k];
	}
	norm(co3);norm(co4);
	double r[3];
	for(int k=0;k<3;k++){
	  r[k] = coor.coor[N2+k]+v2[k]/2-
	    ( coor.coor[N1+k]+ v1[k]/2);
	}
	double dist = norm(r);
#warning SS_senergy modified (dist<7.5) Fri Oct 27 10:16:44 CEST 2006
	if(dist<7.5){
	  double dp1 = dotp(r,co1);
	  double dp2 = dotp(r,co2);
	  double dp3 = dotp(r,co3);
	  double dp4 = dotp(r,co4);
	  norm(v1);norm(v2);
#warning SS_senergy modified todeg 180.0/M_PI Tue Dec  5 14:22:54 CET 2006 
	  double theta = 180.0/M_PI*acos(dotp(v1,v2));
	  
	  //phi
	  double t1[3],t2[3];
	  
	  vectorProduct(v1,r,t1);
	  vectorProduct(v1,t1,t2);
	  double q1 = dotp(v2,t1);
	  double q2 = dotp(v2,t2);
	  double phi = 180.0/M_PI*atan2(q1,q2);
	  
	  double dp = 0.5*(abs(dp1)+abs(dp2)+abs(dp3)+abs(dp4));
	  double hbs = hb_score(dp);
	  double tps = theta_phi_score(theta,phi,loop);
	  //distance 
#warning SS_senergy modified (dscore) Fri Oct 27 10:16:44 CEST 2006
	  double dscore = 0;
	  double dd_ = max(0.0,dist-4.5)-3; //best 4.5 2 because dist
	  if(loop<10){
	    dscore = dd_;
	  }
	  else { dscore = dd_;}
	  //tps score not USED!!!!	  
	  //	  cout<<"DEBUG dscore "<<dscore<<"\n";
	  dimer_pair_score[d1][d2] =  dimer_pair_score[d2][d1] = hbs+tps+dscore;
	}
      }
    }
    
    //Trimm
    double e=0;
    for(int d1 = 0; d1<ndimers; d1++){
      map<int,double> str_int;
      set<double> best;
      for(int d2 = 0;d2<ndimers;d2++){
	if(dimer_strand[d1]!=dimer_strand[d2]){
	  double c = str_int[dimer_strand[d2]];
	  if(dimer_pair_score[d1][d2]<c){
	    str_int[dimer_strand[d2]] = dimer_pair_score[d1][d2];
	    //  if(best.find(dimer_pair_score[d1][d2])!=best.end())
	    //  cout<<"WARNING shadow\n";
	    best.insert(dimer_pair_score[d1][d2]);
	  }
	}
      }
      //Take the two best
      set<double>::iterator iter = best.begin(), end = best.end();
      int count = 0;
      while(iter!=end){
	e+=*iter;
	++iter;
	++count;
	if(count==2) break;
      }
    }
    //    cout<<"ENERGY "<<e<<endl;
    return e;
  }


  double RoseFF::HH_energy(const Coor<double> & coor, const SecStruct & secstruct){
    vector<const SecStructElement * > alpha;
    for(int i=0;i<secstruct.ss_element.size();i++){
      if((secstruct.is_alpha(i)) &&
	 (secstruct.ss_element[i].len>=4)
	 ) {
	alpha.push_back(&secstruct.ss_element[i]);
      }
    }
    double e = 0;
    for(int hh1 = 0; hh1<alpha.size();hh1++){
      for(int hh2 = hh1+1; hh2<alpha.size();hh2++){
	int sep = alpha[hh2]->beg - alpha[hh1]->end;
	for(int d1 = alpha[hh1]->beg;d1<alpha[hh1]->end-3;d1++){
	  double v1[3],m1[3];
	  helix_axe(d1,v1,m1,coor);//begin 4 HHHH
	  for(int d2 = alpha[hh2]->beg;d2<alpha[hh2]->end-3;d2++){
	    double v2[3],m2[3],r[3];
	    helix_axe(d2,v2,m2,coor);
	    for(int i=0;i<3;i++) r[i] = m2[i]-m1[i];
	    double dist = norm(r);
	    if(dist<14){
	      double theta = acos(dotp(v1,v2));
	      double t1[3],t2[3];
	      vectorProduct(v1,r,t1);
	      vectorProduct(v1,t1,t2);
	      double q1 = dotp(v2,t1);
	      double q2 = dotp(v2,t2);
	      double phi = atan2(q1,q2);
	      e += theta_phi_score_hh(theta,phi,sep);
	    }	      
	  }
	}
      }
    }
    return e; 
  }


  double RoseFF::HS_energy(const Coor<double> & coor, const SecStruct & secstruct){
    vector<const SecStructElement * > alpha;
    vector<const SecStructElement * > beta;
    for(int i=0;i<secstruct.ss_element.size();i++){
      if((secstruct.is_alpha(i)) &&
	 (secstruct.ss_element[i].len>=4)
	 ) {
	alpha.push_back(&secstruct.ss_element[i]);
      }
      else if((secstruct.is_beta(i)) &&
	      (secstruct.ss_element[i].len>=2)
	      ) {
	beta.push_back(&secstruct.ss_element[i]);
      }
    }
    if((alpha.size()==0)||(beta.size()==0)) return 0;
    double e = 0;
    for(int hh1 = 0; hh1<alpha.size();hh1++){
      for(int ss2 = 0; ss2<beta.size();ss2++){
	int sep;
	if(alpha[hh1]->beg<beta[ss2]->beg){

	  sep = beta[ss2]->beg - alpha[hh1]->end;
	}
	else {
	  sep = alpha[hh1]->beg - beta[ss2]->end;
	}
	
	for(int d1 = alpha[hh1]->beg;d1<alpha[hh1]->end-3;d1++){
	  //find v1
	  double v1[3],m1[3];
	  helix_axe(d1,v1,m1,coor);//begin 4 HHHH
	  for(int d2 = beta[ss2]->beg;d2<beta[ss2]->end-1;d2++){
	    double v2[3],m2[3],r[3];
	    for(int i=0;i<3;i++){
	      v2[i] = coor.coor[4*C[d2+1]+i]
		-coor.coor[4*N[d2]+i];
	      m2[i] = coor.coor[4*N[d2]+i]+0.5*v2[i];
	    }
	    norm(v2);
	    for(int i=0;i<3;i++) r[i] = m2[i]-m1[i];
	    double dist = norm(r);
	    if(dist<14){
	      double theta = 180.0/M_PI*acos(dotp(v1,v2));
	      double t1[3],t2[3];
	      vectorProduct(v1,r,t1);
	      vectorProduct(v1,t1,t2);
	      double q1 = dotp(v2,t1);
	      double q2 = dotp(v2,t2);
	      double phi = 180.0/M_PI*atan2(q1,q2);
	      e+=theta_phi_score_hs(theta,phi,sep);
	    }
	  }
	}	
      }
    }
    return e;
  }


  double RoseFF::env_energy(const Coor<double> & coor){
    static double data[20][32] = {
      { 10.00, 10.00, 7.26, 6.57, 5.63, 4.87, 4.43, 4.00, 3.65, 3.36, 3.16, 3.12, 3.13, 3.13, 3.14, 3.14, 3.15, 3.14, 3.12, 3.10, 2.99, 3.00, 2.95, 2.96, 2.97, 3.06, 3.18, 3.35, 3.66, 3.96, 4.37, 3.73 },
      { 9.21, 9.21, 8.52, 6.91, 6.81, 6.03, 5.60, 5.20, 4.93, 4.69, 4.19, 4.05, 3.83, 3.58, 3.50, 3.30, 3.08, 3.02, 2.87, 2.71, 2.66, 2.61, 2.52, 2.58, 2.65, 2.80, 2.90, 3.17, 3.54, 3.96, 4.57, 4.12 },
      { 10.00, 10.00, 7.42, 6.57, 5.24, 4.19, 3.68, 3.25, 2.95, 2.74, 2.59, 2.52, 2.52, 2.57, 2.63, 2.74, 2.85, 3.01, 3.14, 3.29, 3.50, 3.66, 3.83, 4.02, 4.17, 4.40, 4.63, 4.85, 5.18, 5.52, 5.88, 5.13 },
      { 10.00, 10.00, 7.13, 6.38, 5.17, 4.29, 3.82, 3.28, 2.90, 2.68, 2.49, 2.43, 2.43, 2.46, 2.55, 2.67, 2.80, 2.99, 3.19, 3.41, 3.65, 3.87, 4.10, 4.33, 4.56, 4.80, 5.12, 5.36, 5.99, 6.12, 6.73, 5.30 },
      { 10.00, 10.00, 7.82, 7.13, 6.07, 5.38, 5.15, 4.77, 4.43, 4.09, 3.90, 3.75, 3.64, 3.48, 3.32, 3.13, 2.95, 2.81, 2.75, 2.65, 2.63, 2.59, 2.61, 2.63, 2.75, 2.99, 3.19, 3.55, 3.97, 4.46, 5.02, 4.85 },
      { 10.00, 10.00, 7.26, 6.50, 5.32, 4.18, 3.49, 3.14, 2.99, 2.90, 2.86, 2.89, 2.94, 3.00, 3.02, 3.06, 3.14, 3.17, 3.22, 3.25, 3.30, 3.29, 3.34, 3.32, 3.38, 3.49, 3.58, 3.76, 3.94, 4.27, 4.62, 3.96 },
      { 10.00, 9.21, 6.91, 6.12, 5.43, 4.64, 4.37, 3.88, 3.52, 3.38, 3.14, 3.03, 2.90, 2.89, 2.81, 2.71, 2.71, 2.70, 2.82, 2.82, 2.97, 3.08, 3.17, 3.26, 3.39, 3.65, 3.84, 4.18, 4.51, 4.87, 5.45, 5.04 },
      { 10.00, 10.00, 8.11, 7.13, 6.21, 5.45, 5.24, 4.96, 4.55, 4.30, 3.90, 3.82, 3.63, 3.49, 3.28, 3.18, 3.04, 2.89, 2.86, 2.80, 2.70, 2.61, 2.59, 2.61, 2.68, 2.85, 3.05, 3.37, 3.74, 4.19, 4.78, 4.41 },
      { 10.00, 10.00, 7.13, 6.12, 5.08, 4.37, 3.90, 3.44, 3.08, 2.82, 2.60, 2.48, 2.44, 2.45, 2.49, 2.54, 2.67, 2.78, 3.05, 3.29, 3.51, 3.88, 4.09, 4.41, 4.73, 5.02, 5.40, 5.81, 6.32, 6.65, 7.13, 5.63 },
      { 10.00, 10.00, 7.82, 7.01, 6.17, 5.52, 5.24, 4.82, 4.46, 4.19, 3.85, 3.71, 3.54, 3.39, 3.18, 3.06, 2.94, 2.83, 2.77, 2.70, 2.65, 2.62, 2.60, 2.66, 2.77, 2.98, 3.24, 3.58, 4.00, 4.47, 5.13, 4.72 },
      { 10.00, 10.00, 6.07, 5.74, 5.52, 5.02, 4.61, 4.41, 4.10, 3.78, 3.55, 3.45, 3.27, 3.16, 3.12, 3.01, 2.90, 2.84, 2.89, 2.83, 2.85, 2.84, 2.80, 2.82, 2.87, 3.07, 3.25, 3.51, 3.91, 4.41, 4.98, 4.37 },
      { 10.00, 10.00, 7.42, 6.57, 5.24, 4.26, 3.77, 3.38, 3.07, 2.92, 2.71, 2.66, 2.58, 2.66, 2.71, 2.70, 2.83, 2.92, 3.06, 3.26, 3.29, 3.51, 3.57, 3.73, 3.83, 4.05, 4.25, 4.54, 4.83, 5.36, 5.71, 4.96 },
      { 10.00, 10.00, 6.65, 5.95, 4.76, 4.01, 3.55, 3.15, 2.90, 2.75, 2.71, 2.69, 2.70, 2.80, 2.87, 2.91, 2.98, 3.06, 3.17, 3.26, 3.31, 3.45, 3.48, 3.64, 3.76, 3.89, 4.10, 4.25, 4.62, 5.08, 5.57, 4.85 },
      { 10.00, 10.00, 7.26, 6.65, 5.45, 4.64, 4.19, 3.67, 3.25, 2.96, 2.67, 2.61, 2.54, 2.55, 2.57, 2.63, 2.68, 2.79, 2.94, 3.17, 3.30, 3.48, 3.62, 3.84, 4.10, 4.33, 4.55, 4.88, 5.17, 5.68, 6.17, 5.36 },
      { 10.00, 10.00, 7.82, 6.50, 5.52, 4.78, 4.35, 3.93, 3.51, 3.17, 2.86, 2.75, 2.62, 2.56, 2.56, 2.50, 2.58, 2.64, 2.76, 2.96, 3.14, 3.33, 3.55, 3.81, 4.02, 4.37, 4.69, 5.13, 5.63, 6.03, 6.65, 5.43 },
      { 10.00, 10.00, 7.01, 6.21, 5.30, 4.40, 3.99, 3.55, 3.18, 2.95, 2.77, 2.75, 2.71, 2.77, 2.81, 2.87, 2.93, 3.00, 3.08, 3.16, 3.25, 3.32, 3.39, 3.43, 3.57, 3.70, 3.86, 4.05, 4.35, 4.68, 5.08, 4.35 },
      { 10.00, 10.00, 7.42, 6.44, 5.45, 4.79, 4.27, 3.85, 3.44, 3.17, 2.95, 2.80, 2.74, 2.72, 2.75, 2.74, 2.80, 2.88, 3.00, 3.08, 3.16, 3.25, 3.38, 3.37, 3.51, 3.64, 3.76, 3.94, 4.27, 4.58, 5.05, 4.34 },
      { 10.00, 10.00, 7.60, 7.01, 6.03, 5.50, 5.17, 4.79, 4.37, 3.98, 3.77, 3.62, 3.51, 3.39, 3.25, 3.15, 3.06, 2.99, 2.92, 2.88, 2.78, 2.71, 2.69, 2.72, 2.73, 2.84, 3.00, 3.29, 3.61, 4.02, 4.47, 4.15 },
      { 9.21, 9.21, 7.60, 7.26, 6.50, 5.47, 5.01, 4.88, 4.37, 4.12, 3.90, 3.60, 3.46, 3.34, 3.18, 3.01, 2.83, 2.73, 2.70, 2.65, 2.55, 2.63, 2.66, 2.75, 2.89, 3.04, 3.38, 3.79, 4.22, 4.98, 5.26, 4.46 },
      { 10.00, 10.00, 7.82, 6.81, 5.99, 5.60, 4.99, 4.73, 4.34, 4.00, 3.73, 3.56, 3.44, 3.29, 3.11, 2.94, 2.84, 2.72, 2.68, 2.62, 2.62, 2.62, 2.65, 2.78, 2.93, 3.19, 3.45, 3.81, 4.33, 4.93, 5.45, 4.73 }
    };

    double e = 0;
    vector<int> nb;
    nb.resize(frag_size);

    for(int a = 0;a<frag_size;a++){
      for(int b = 0;b<frag_size;b++){
	double dx,dy,dz,d2;
	int ipos = Coor<double>::DIM*CB[a];
	int jpos = Coor<double>::DIM*CB[b];
	dx = coor.coor[ipos] - coor.coor[jpos];
	dy = coor.coor[ipos+1] - coor.coor[jpos+1];
	dz = coor.coor[ipos+2] - coor.coor[jpos+2];
	d2 = dx*dx+dy*dy+dz*dz;
	if(d2<100){
	  ++nb[a];
	  ++nb[b];
	}	
      }
    }

    for(int a=0;a<frag_size;a++){
      int q = nb[a];
      if(q>31) q=31;
      e += data[idx[a]][q];
    }

    return w_env*e;
  }


  double RoseFF::pair_energy(const Coor<double> & coor){
    static double pair0[20][20] = {
      { -0.2968, 0.3676, 0.2450, 0.2169, 0.0411, 0.3456, 0.1857, -0.2287, 0.1832, -0.1172, -0.0273, 0.3457, -0.0087, 0.1109, 0.2204, 0.2165, 0.1084, -0.3303, 0.1300, 0.0381 },
      { 0.3358, -2.3997, 0.6197, 0.5811, 0.1896, 0.2426, 0.5343, 0.4260, 0.3839, 0.2925, 0.0810, 0.4516, 0.2152, 0.4083, 0.4956, 0.1840, 0.4080, 0.2374, 0.2808, 0.2453 },
      { 0.2634, 0.5565, -0.5933, 0.2495, 0.3852, -0.1557, -0.2091, 0.4182, -0.5871, 0.4961, 0.4083, -0.3496, -0.0934, 0.1964, -0.6039, -0.3097, -0.2844, 0.3238, 0.3541, 0.3885 },
      { 0.2688, 0.5560, 0.1855, -0.8756, 0.3191, 0.1077, 0.0412, 0.2173, -0.8232, 0.3896, 0.3041, -0.1073, -0.1281, -0.0843, -0.6607, -0.1180, -0.0900, 0.2602, 0.2086, 0.4023 },
      { -0.0126, 0.0923, 0.4628, 0.2536, -0.4862, 0.0200, 0.1141, -0.0562, 0.2663, -0.1556, -0.1265, 0.1902, -0.0297, 0.1374, 0.1070, 0.1069, 0.3144, -0.0954, -0.1925, -0.1774 },
      { 0.2816, 0.1862, -0.0988, 0.1370, 0.0164, -0.1942, 0.1038, -0.0756, -0.0036, -0.0853, -0.0524, -0.2629, -0.0243, -0.0838, 0.0263, 0.0525, -0.1396, 0.0548, 0.1335, 0.0840 },
      { 0.2371, 0.0539, -0.1815, -0.0282, 0.1169, 0.0518, -0.9004, 0.2530, 0.0612, 0.0705, 0.1240, -0.0521, -0.0843, -0.1668, -0.1381, -0.2014, -0.0471, 0.1352, 0.0835, -0.0212 },
      { -0.2218, 0.3750, 0.4421, 0.1995, -0.1708, -0.0516, 0.1549, -0.5392, 0.3077, -0.1400, -0.0033, 0.4643, 0.3892, 0.2230, 0.3442, 0.2925, 0.3987, -0.1507, 0.0288, -0.1059 },
      { 0.2425, 0.6657, -0.7520, -0.8286, 0.2615, 0.0166, 0.0298, 0.3246, -0.5538, 0.4062, 0.3412, -0.1706, 0.0451, -0.3227, 0.5854, -0.0114, -0.0706, 0.2628, 0.0725, -0.0636 },
      { -0.1567, 0.1537, 0.4104, 0.4585, -0.1823, 0.0078, 0.0551, -0.1702, 0.4149, -0.4134, 0.0385, 0.4063, 0.0548, 0.2532, 0.2576, 0.3351, 0.2688, -0.1810, -0.0749, -0.0497 },
      { 0.0001, 0.1786, 0.1746, 0.3925, -0.0893, -0.0169, 0.1250, -0.0972, 0.3192, -0.0475, -1.0073, 0.1830, -0.1063, 0.1210, 0.3139, 0.0839, 0.0957, -0.0084, 0.0232, -0.1152 },
      { 0.3636, 0.4107, -0.3731, -0.0587, 0.1373, -0.2605, -0.1087, 0.3712, -0.1321, 0.4309, 0.1837, -0.8960, 0.1134, -0.1598, -0.0548, -0.2588, -0.2715, 0.4393, 0.1536, 0.1562 },
      { 0.0025, 0.3188, -0.1732, -0.0494, -0.0975, -0.0795, -0.1524, 0.4265, 0.0600, 0.2440, -0.0137, -0.0352, -0.5577, -0.0632, -0.0754, 0.0129, 0.1603, 0.2573, -0.4923, -0.0430 },
      { 0.1655, 0.4917, 0.0115, -0.0657, 0.1483, -0.0638, -0.1238, 0.1729, -0.2784, 0.1999, 0.1734, -0.0464, -0.1359, -0.7277, -0.1990, -0.1218, -0.1881, 0.1716, 0.1028, 0.0843 },
      { 0.2141, 0.3402, -0.4703, -0.6353, 0.1554, -0.0158, -0.0628, 0.3311, 0.5053, 0.1360, 0.5053, -0.0415, 0.1068, -0.2048, -0.6554, -0.0903, -0.0570, 0.2625, -0.1981, 0.0471 },
      { 0.2243, 0.3110, -0.3147, -0.0858, 0.2279, 0.0094, -0.1174, 0.2180, -0.0871, 0.3195, 0.0638, -0.2797, -0.0334, -0.0988, -0.1069, -0.5562, -0.2233, 0.2055, 0.1892, 0.2037 },
      { 0.0895, 0.4083, -0.1563, -0.1373, 0.3202, -0.1482, -0.1196, 0.4745, -0.0607, 0.3349, 0.0042, -0.2473, 0.0975, -0.1482, 0.0017, -0.2670, -0.6919, 0.3202, 0.4024, 0.2912 },
      { -0.3300, 0.3530, 0.3552, 0.1495, -0.1174, 0.0404, 0.1598, -0.0747, 0.2110, -0.1672, -0.0372, 0.4848, 0.2322, 0.1759, 0.2483, 0.1761, 0.3908, -0.2359, 0.1519, -0.0427 },
      { 0.2586, 0.0494, 0.5899, 0.5670, -0.2674, 0.1458, -0.0431, -0.0602, -0.0833, -0.2222, -0.2957, 0.0194, -0.2453, 0.0283, 0.0032, 0.3295, 0.3352, 0.0892, -1.3140, -0.2844 },
      { 0.0850, 0.0862, 0.4142, 0.3235, -0.1405, 0.0463, -0.1155, -0.1409, 0.0044, -0.0952, -0.0957, 0.1775, 0.0188, 0.1212, -0.0113, 0.1688, 0.2953, 0.0110, -0.2432, -0.6852 }
    };

    static double pair1[20][20] = {
      { -0.0752, 0.2037, 0.1188, 0.1208, -0.0584, -0.0992, 0.0395, -0.0650, 0.1953, -0.1507, -0.0841, 0.1079, 0.1440, 0.1387, 0.0840, 0.0970, 0.0698, -0.0090, 0.0392, 0.0392 },
      { 0.0938, -1.0733, 0.2138, 0.2760, -0.0789, 0.0420, -0.2157, 0.0686, 0.0487, 0.0326, -0.0108, 0.0449, -0.0763, 0.0058, 0.0971, -0.0372, 0.0917, 0.0304, -0.1203, 0.0516 },
      { 0.0801, 0.1846, -0.2026, 0.1461, 0.2898, -0.3111, -0.3983, 0.4348, -0.6587, 0.4630, 0.2615, -0.3749, -0.1712, -0.2363, -0.5387, -0.2516, -0.0990, 0.3852, 0.1626, 0.1642 },
      { 0.1324, 0.2195, 0.1754, 0.1346, 0.2107, -0.0121, -0.2557, 0.2602, -0.5808, 0.1346, 0.0274, -0.1096, -0.1425, -0.1096, -0.5415, -0.2021, -0.1634, 0.2338, 0.1808, 0.0420 },
      { -0.0364, -0.0947, 0.2734, 0.2308, -0.2456, 0.0781, 0.0043, -0.0803, 0.1055, -0.0642, -0.1432, 0.1302, -0.0149, 0.0581, 0.1017, 0.0690, 0.1311, -0.0251, -0.0922, -0.0408 },
      { -0.0985, 0.0510, -0.2859, -0.0562, 0.1020, -0.3586, -0.1014, 0.3422, -0.0989, 0.3302, 0.0768, -0.2261, -0.1269, -0.1346, -0.1157, -0.1961, -0.0706, 0.2573, -0.0348, 0.0514 },
      { 0.0401, -0.1981, -0.4631, -0.1876, 0.0589, -0.1421, -0.4257, 0.3003, 0.1792, 0.2221, -0.0746, -0.0620, -0.0865, -0.0299, 0.0068, -0.1795, -0.0765, 0.1821, -0.0056, 0.0223 },
      { -0.0554, 0.0710, 0.3747, 0.2848, -0.0784, 0.3653, 0.2351, -0.3479, 0.1668, -0.1817, -0.0688, 0.3487, 0.2653, 0.2228, 0.2777, 0.2383, 0.0482, -0.2513, 0.1294, 0.0075 },
      { 0.1486, 0.1223, -0.6354, -0.6082, 0.0992, -0.1470, 0.1305, 0.1669, 0.0356, 0.2175, 0.1971, -0.1924, 0.0677, -0.0545, 0.3415, -0.0708, -0.0640, 0.2229, 0.0320, -0.1539 },
      { -0.1322, 0.0310, 0.4299, 0.1466, -0.0692, 0.3365, 0.2234, -0.1886, 0.1880, -0.2458, -0.0444, 0.3954, 0.1781, 0.1469, 0.1318, 0.2214, 0.1539, -0.1963, 0.0374, -0.0010 },
      { -0.0615, 0.0163, 0.1831, -0.0352, -0.1320, 0.1385, 0.0214, -0.0655, 0.1709, -0.0334, -0.3079, 0.1831, 0.0609, 0.1137, 0.0850, 0.1371, 0.0593, -0.0757, 0.0106, -0.0675 },
      { 0.1210, -0.0303, -0.3083, -0.0954, 0.1525, -0.2222, -0.1085, 0.3760, -0.1908, 0.3966, 0.1910, -0.5818, -0.2494, -0.2325, -0.0771, -0.2977, -0.1914, 0.3631, 0.0064, -0.0087 },
      { 0.1322, -0.0453, -0.1219, -0.0707, 0.0531, -0.1686, -0.0936, 0.2057, 0.0878, 0.1683, 0.0621, -0.2054, -0.2716, -0.1429, -0.1230, -0.1119, -0.0418, 0.2010, -0.2525, -0.1544 },
      { 0.1101, -0.0318, -0.1527, -0.1156, 0.0960, -0.1283, -0.1166, 0.2620, -0.0208, 0.1396, 0.0984, -0.2426, -0.1872, -0.2600, -0.0346, -0.2407, -0.1079, 0.1741, -0.0987, 0.0754 },
      { 0.0911, 0.0476, -0.5269, -0.5731, 0.0436, -0.1134, 0.0024, 0.2564, 0.3374, 0.1498, 0.1435, -0.0274, -0.0986, -0.0479, 0.0520, -0.0756, 0.0612, 0.1681, -0.2406, -0.0737 },
      { 0.0785, -0.0293, -0.2341, -0.2246, 0.1147, -0.1952, -0.1708, 0.2576, -0.0544, 0.2345, 0.0582, -0.2797, -0.0903, -0.2266, -0.0837, -0.1970, -0.1100, 0.2508, -0.0887, 0.0358 },
      { 0.0626, -0.0034, -0.0971, -0.1406, 0.1154, -0.0622, -0.0173, 0.0800, -0.0756, 0.1388, 0.0999, -0.2040, -0.0618, -0.1495, 0.0339, -0.0800, -0.1786, 0.0536, 0.0467, 0.0766 },
      { -0.0180, 0.0396, 0.3510, 0.2503, -0.0402, 0.2846, 0.2033, -0.2593, 0.1487, -0.2044, -0.0113, 0.3706, 0.2230, 0.2221, 0.1749, 0.2841, 0.0316, -0.3000, 0.2100, 0.0604 },
      { 0.1236, -0.1529, 0.1420, 0.1919, -0.1264, -0.0097, 0.0586, 0.1182, -0.0227, 0.0401, -0.0316, -0.0418, -0.2967, -0.0752, -0.2745, -0.0659, 0.1745, 0.1810, -0.3958, -0.2124 },
      { 0.0542, 0.0141, 0.0933, 0.0566, -0.0797, 0.0406, 0.0178, 0.0361, -0.1596, 0.0056, -0.1160, 0.0204, -0.1238, 0.0463, -0.0763, 0.0103, 0.0768, 0.0706, -0.1565, -0.1235 }
    };

    static double pair2[20][20] = {
      { -0.1423, 0.0802, 0.0194, 0.0308, 0.0795, -0.0054, 0.0345, -0.0164, 0.0606, 0.0113, 0.0546, 0.0624, 0.0196, 0.0315, 0.0014, 0.0437, -0.0157, -0.0835, 0.0724, 0.0873 },
      { 0.0768, -0.6483, 0.0839, 0.0882, -0.0391, -0.0199, -0.1425, 0.0459, 0.0119, 0.0736, -0.0433, -0.0051, 0.0300, -0.0525, 0.0316, -0.0843, -0.0224, 0.0532, -0.0436, 0.0153 },
      { 0.0361, 0.0502, -0.1325, 0.0241, 0.1662, -0.0972, -0.1094, 0.1658, -0.2637, 0.1913, 0.1544, -0.1123, -0.0929, -0.0486, -0.2483, -0.0205, -0.0276, 0.1105, 0.0120, -0.0680 },
      { 0.0392, 0.0730, 0.0107, -0.0679, 0.1073, 0.0207, -0.0736, 0.0843, -0.3370, 0.1331, 0.1520, -0.0182, -0.1147, 0.0034, -0.2340, 0.0243, 0.0245, 0.0535, -0.0397, -0.0695 },
      { 0.0898, -0.0211, 0.1400, 0.1326, -0.1975, 0.1230, 0.0394, -0.1427, 0.1026, -0.1486, -0.1096, 0.1052, 0.1461, 0.0533, 0.0896, 0.0888, 0.0706, -0.0726, -0.1189, -0.0771 },
      { -0.0390, -0.0028, -0.0448, -0.0050, 0.1646, -0.1681, -0.0138, 0.0696, 0.0178, 0.1326, 0.0134, -0.0632, -0.0576, -0.0524, -0.0051, -0.0806, -0.0780, 0.0243, 0.0992, 0.0924 },
      { 0.0680, -0.1119, -0.1195, -0.1167, 0.0322, -0.0398, -0.2860, 0.0882, 0.0795, 0.1030, 0.0077, -0.0215, -0.0282, -0.0439, 0.0393, -0.0716, -0.0169, 0.0760, -0.1265, -0.0267 },
      { -0.0072, 0.0739, 0.1589, 0.1070, -0.1412, 0.1108, 0.1226, -0.2060, 0.0531, -0.1419, -0.1115, 0.1364, 0.1760, 0.1087, 0.0931, 0.0726, 0.0800, -0.0679, 0.0386, -0.0207 },
      { 0.0526, 0.0338, -0.3024, -0.3214, 0.0982, -0.0165, 0.1084, 0.0602, -0.1411, 0.1114, 0.2046, -0.1358, -0.0726, 0.0010, 0.1547, -0.0062, 0.0173, 0.0549, 0.1530, 0.0621 },
      { 0.0216, 0.0517, 0.1818, 0.1050, -0.1519, 0.1292, 0.1099, -0.1292, 0.1242, -0.2718, -0.1089, 0.2026, 0.1127, 0.0797, 0.0595, 0.1043, 0.1101, -0.0412, 0.0172, 0.0201 },
      { 0.0401, -0.0136, 0.1669, 0.0917, -0.1277, 0.0561, -0.0075, -0.1234, 0.1427, -0.0968, -0.2344, 0.0626, 0.0794, 0.0367, 0.0262, 0.1070, 0.0459, -0.0430, -0.0102, -0.0216 },
      { 0.0585, 0.0114, -0.1009, 0.0077, 0.1000, -0.0633, -0.0425, 0.1413, -0.1185, 0.2006, 0.0673, -0.2532, -0.1387, -0.0735, -0.0180, -0.1405, -0.0827, 0.1105, -0.0313, -0.0564 },
      { 0.0457, 0.0232, -0.0907, -0.1120, 0.1302, -0.0586, -0.0759, 0.2112, -0.0879, 0.1354, 0.0705, -0.1098, -0.1825, -0.1267, -0.1384, -0.0394, -0.0401, 0.1405, -0.0998, -0.0711 },
      { 0.0011, -0.0429, -0.0490, 0.0159, 0.1045, -0.0612, -0.0789, 0.1272, -0.0039, 0.0672, 0.0547, -0.1065, -0.1242, -0.1951, -0.0227, -0.0269, -0.0050, 0.0600, -0.0244, 0.0076 },
      { -0.0044, 0.0168, -0.2339, -0.2519, 0.1254, -0.0234, -0.0009, 0.1224, 0.1442, 0.0458, 0.0591, 0.0086, -0.1221, 0.0006, -0.0286, -0.0167, 0.0164, 0.0653, 0.0179, 0.0274 },
      { 0.0414, -0.0961, -0.0300, 0.0477, 0.0705, -0.0987, -0.0723, 0.0772, -0.0270, 0.1148, 0.0607, -0.1278, -0.0300, -0.0407, -0.0008, -0.1385, -0.0941, 0.0544, 0.0249, 0.0240 },
      { -0.0108, -0.0628, 0.0005, 0.0418, 0.0356, -0.0625, 0.0141, 0.0575, 0.0147, 0.1033, 0.0229, -0.1001, 0.0010, -0.0377, 0.0341, -0.0790, -0.1023, 0.0041, 0.0050, 0.0058 },
      { -0.0838, 0.0629, 0.1058, 0.0382, -0.0791, 0.0527, 0.0578, -0.0675, 0.0495, -0.0399, -0.0528, 0.1047, 0.1393, 0.0517, 0.0500, 0.0399, 0.0327, -0.0949, 0.0381, 0.0222 },
      { 0.0800, -0.0789, 0.0288, 0.0232, -0.0864, 0.0837, -0.0561, 0.0278, 0.1162, 0.0081, -0.0667, 0.0207, -0.0706, -0.0404, 0.0903, -0.0200, -0.0402, 0.0171, -0.2900, -0.1274 },
      { 0.0748, 0.0236, -0.0755, -0.0541, -0.0788, 0.0671, -0.0497, 0.0281, 0.0336, 0.0150, 0.0111, -0.0427, -0.0531, 0.0242, 0.0150, 0.0279, -0.0081, 0.0250, -0.1309, -0.0841 }
    };

    static double pair3[20][20] = {
      { -0.1004, 0.0696, 0.0571, 0.0340, 0.0366, -0.0496, 0.0168, 0.0038, 0.0682, -0.0209, -0.0151, 0.0394, -0.0125, 0.0402, 0.0098, 0.0270, -0.0060, -0.0223, 0.0385, 0.0549 },
      { 0.0545, -0.5672, 0.0521, 0.0840, -0.0165, -0.0227, -0.0725, 0.0389, 0.0393, 0.0458, -0.0082, -0.0034, -0.0492, -0.0322, 0.0211, -0.0649, -0.0351, 0.0349, 0.0430, -0.0027 },
      { 0.0373, 0.0713, -0.0611, 0.0099, 0.0199, -0.0401, 0.0014, -0.0030, -0.0282, 0.0280, 0.0249, -0.0158, -0.0513, -0.0031, -0.0351, -0.0011, 0.0168, 0.0305, -0.0447, -0.0238 },
      { 0.0379, 0.0428, -0.0135, -0.1065, 0.0251, 0.0726, 0.0210, -0.0123, -0.0826, -0.0412, 0.0016, 0.0263, -0.0258, 0.0355, -0.0633, 0.0753, 0.0604, 0.0002, 0.0679, -0.0121 },
      { 0.0246, -0.0083, 0.0160, 0.0245, -0.1080, 0.0380, -0.0315, 0.0002, -0.0228, -0.0090, -0.0343, 0.0045, 0.0419, -0.0267, 0.0185, 0.0107, 0.0189, 0.0299, -0.0791, -0.0401 },
      { -0.0516, -0.0041, -0.0411, 0.0586, 0.0228, -0.1226, -0.0100, 0.0556, 0.1084, 0.0713, 0.0277, -0.0406, -0.0330, 0.0143, 0.0195, -0.0496, -0.0262, -0.0055, -0.0146, 0.0122 },
      { 0.0277, -0.1384, -0.0028, 0.0482, -0.0183, -0.0406, -0.2029, 0.0509, 0.0901, 0.0008, -0.0559, -0.0052, -0.0594, -0.0054, 0.0214, 0.0237, 0.0231, 0.0428, -0.1159, -0.0485 },
      { 0.0010, 0.0814, 0.0203, -0.0189, -0.0020, 0.0400, 0.0530, -0.0884, -0.0899, -0.0005, -0.0123, -0.0107, 0.1305, 0.0258, 0.0341, 0.0081, 0.0009, -0.0315, 0.1051, 0.0161 },
      { 0.0668, -0.0082, -0.0500, -0.0948, -0.0016, 0.0683, 0.1048, -0.0771, -0.1935, -0.0078, 0.0700, -0.0185, -0.0007, 0.0416, 0.0862, 0.0486, 0.0426, -0.0014, 0.1031, -0.0077 },
      { -0.0272, 0.0309, 0.0280, -0.0442, 0.0050, 0.0789, -0.0112, 0.0060, -0.0231, -0.0743, 0.0033, 0.0646, 0.0685, -0.0282, -0.0601, 0.0247, 0.0292, 0.0050, 0.0177, 0.0614 },
      { -0.0005, -0.0232, 0.0152, -0.0251, -0.0474, 0.0381, -0.0484, -0.0152, 0.0118, -0.0079, -0.1434, 0.0727, 0.0232, -0.0167, -0.0155, 0.0476, 0.0279, 0.0310, -0.0076, -0.0324 },
      { 0.0623, -0.0339, -0.0387, 0.0517, -0.0223, -0.0278, 0.0240, -0.0091, 0.0038, 0.0714, -0.0009, -0.1482, -0.0148, -0.0237, 0.0625, -0.0684, -0.0307, 0.0401, -0.0423, -0.0750 },
      { 0.0216, -0.0548, -0.0595, -0.0445, 0.0591, -0.0393, -0.0464, 0.1159, -0.0002, 0.0630, 0.0395, -0.0351, -0.1484, -0.0621, -0.0638, -0.0182, 0.0093, 0.0444, -0.0198, -0.0442 },
      { 0.0401, -0.0190, 0.0045, 0.0440, -0.0081, 0.0022, -0.0238, 0.0268, 0.0381, -0.0424, 0.0076, -0.0135, -0.0358, -0.1051, 0.0227, -0.0196, 0.0175, 0.0284, -0.0248, -0.0262 },
      { 0.0071, 0.0189, -0.0089, -0.0794, 0.0123, 0.0045, -0.0074, 0.0189, 0.1034, -0.0380, -0.0233, 0.0576, -0.0492, 0.0225, -0.1107, 0.0346, 0.0380, 0.0137, 0.0383, 0.0254 },
      { 0.0309, -0.0739, 0.0133, 0.0793, -0.0163, -0.0345, 0.0085, 0.0181, 0.0603, 0.0187, 0.0297, -0.0348, -0.0229, -0.0173, 0.0327, -0.0811, -0.0537, 0.0144, -0.0190, -0.0445 },
      { -0.0034, -0.0249, 0.0006, 0.0532, -0.0059, -0.0230, 0.0386, -0.0053, 0.0409, 0.0305, 0.0282, -0.0273, -0.0100, 0.0235, 0.0532, -0.0236, -0.0476, -0.0354, 0.0083, -0.0264 },
      { -0.0228, 0.0444, 0.0143, 0.0127, 0.0389, 0.0061, 0.0285, -0.0365, -0.0051, -0.0039, 0.0191, 0.0457, 0.0432, 0.0282, 0.0215, 0.0179, -0.0353, -0.0812, 0.0598, 0.0674 },
      { 0.0594, 0.0041, 0.0257, 0.0428, -0.0873, -0.0061, -0.0820, 0.0770, 0.0941, 0.0094, -0.0167, -0.0486, 0.0198, -0.0329, 0.0250, -0.0273, -0.0359, 0.0538, -0.3227, -0.1280 },
      { 0.0481, 0.0553, -0.0228, -0.0069, -0.0539, 0.0247, -0.0284, 0.0231, -0.0458, 0.0641, -0.0588, -0.0699, -0.0527, -0.0083, 0.0065, -0.0393, 0.0044, 0.0662, -0.1812, -0.0679 }
    };

    double e = 0;

    
    for(int a=0;a<frag_size;a++){
      for(int b=a+9;b<frag_size;b++){
	double dx,dy,dz,d2;
	int ipos = Coor<double>::DIM*CB[a];
	int jpos = Coor<double>::DIM*CB[b];
	dx = coor.coor[ipos] - coor.coor[jpos];
	dy = coor.coor[ipos+1] - coor.coor[jpos+1];
	dz = coor.coor[ipos+2] - coor.coor[jpos+2];
	d2 = dx*dx+dy*dy+dz*dz;
	if(d2<25){
	  e += pair0[idx[a]][idx[b]];
	}
	else if((25<=d2)&&(d2<56.25)){
	  e += pair1[idx[a]][idx[b]];	  
	}
	else if((56.25<=d2)&&(d2<100)){
	  e += pair2[idx[a]][idx[b]];	  
	}
	else if((100<=d2)&&(d2<144)){
	  e += pair3[idx[a]][idx[b]];	  
	}
      }
    }
    return w_pair*e;
  }

  double RoseFF::theta_phi_score(double theta, double phi, int sep){
    phi = phi+180;
    int ix = (int)( theta/18);
    int iy = (int) (phi/36);
    if(sep>10) return theta_phi_SS_hi[ix][iy];
    else return theta_phi_SS_low[ix][iy];
  }


  double RoseFF::theta_phi_score_hh(double theta, double phi, int sep){
    phi = phi+180;
    int ix = (int)( theta/18);
    int iy = (int) (phi/36);
    if(sep>10) return theta_phi_HH_hi[ix][iy];
    else return theta_phi_HH_low[ix][iy];
  }


  double RoseFF::theta_phi_score_hs(double theta, double phi, int sep){
    phi = phi+180;
    int ix = (int)( theta/18);
    int iy = (int) (phi/36);
    if(sep>10) return theta_phi_HS_hi[ix][iy];
    else return theta_phi_HS_low[ix][iy];    
  }

}


// Local Variables:
// mode:C++
// End: