/***************************************************************************
                          rosemc.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/rosemc.h>
#include <algorithm/sec_packing_data.h>
#include <secstruct.h>
#include <kernel/rmsd.h>
#include <algorithm/molalgo.h>
extern double score_silent(const Almost::Molecules & mols, string file);
extern int sys_bflock(string f, int s);
extern int sys_uflock(string f, int fd);

namespace Almost {

  void RoseMC::thread(const Molecules & molecules, string mrg_file){
    vector<double> phi;
    vector<double> psi;
    string SS;
    {
      //Read MRG file
      ifstream in;
      in.open(mrg_file.c_str());
      istream_iterator<string> iter(in),end;
      while(iter!=end){
	++iter;
	++iter;
	phi.push_back(atof(iter->c_str())*M_PI/180.0);++iter;
	psi.push_back(atof(iter->c_str())*M_PI/180.0);++iter;
	++iter;
	for(int i=0;i<6;i++) ++iter;
	SS += *iter;++iter;
      }
    }
    Coor<double> coor;
    Protein p = molecules.protein(0);
    int size = p.fragment_size();
    coor.resize(p.atom_size());
    BondRotation<double> br(p);
    vector<int> N;
    vector<int> C;
    vector<int> CA;
    int b = p.atom_offset();
    int e = p.atom_size()+p.atom_offset();
    for(int i = b; i<e;i++){
      if(p[i].name()=="N")  N.push_back(i);
      if(p[i].name()=="C")  C.push_back(i);
      if(p[i].name()=="CA") CA.push_back(i); 
    }

    int Size = phi.size();
    if(size>Size){
      cout<<"Source protein to small ..."<<size<<" "<<Size<<" "<<mrg_file<<endl;
      return;
    }

    for(int i=0;i<= Size - size; i++){
      for(int j=0;j<size;j++){
	if(phi[i+j]!=999){
	  set_torsion2(p,C[j-1],N[j],CA[j],C[j],phi[i+j],br);
	}
	if(psi[i+j]!=999){
	  set_torsion2(p,N[j],CA[j],C[j],N[j+1],psi[j],br);
	}
      }

      //Eval energy:
      for(int a=0;a<p.atom_size();a++){
	int ipos =Coor<double>::DIM*a;
	coor.coor[ipos] = p[a][0];
	coor.coor[ipos+1] = p[a][1];
	coor.coor[ipos+2] = p[a][2];
      }

      //Ene
      double env = rose_ff->env_energy(coor);
      double pair= rose_ff->pair_energy(coor);
      //      double ff  = rose_ff->ff_energy(coor);
      //SS
      string ss_sub = string(SS.begin()+i,SS.begin()+i+size);
      double ss     = rose_ff->ss_energy(coor,ss_sub);
      double bhb    = rose_ff->bhb_energy(coor,ss_sub);
      cout<<"Energy thread "<<i+1<<" ";
      cout<<env<<" "
	  <<pair<<" "
	//	  <<ff<<" "
	  <<ss<<" "
	  <<bhb<<" "
	  <<env+pair
	//+ff
	+ss+bhb<<" "<<ss_sub
	  <<"\n";
    }

    
  }

  Molecules RoseMC::run(const Molecules & molecules, const MDB & mdb, 
			const RoseOpts options, int sc)
  {
    sc_ = sc;
    if(options.dcd != "none"){
      dcd = new DCDOStream(options.dcd);
    }
    else {
      dcd = 0;
    }

    srand48(options.seed);

    steps = options.steps;
    nprint = options.nprint;
    ndump = options.ndump;
    p_plain = options.p_plain;
    p_3 = options.p_3;
    begin_temp = options.begin_temp;
    end_temp = options.end_temp;
    sa_steps = options.sa_steps;
    sa_cycles = options.sa_cycles;
    pdb = options.pdb;
    best_pdb = options.best_pdb;
    setup(molecules);
    setup_rosemoves(molecules,options);
    prot.resize(molecules.fragment_size());
    prot_buffer.resize(molecules.fragment_size());
    ss_match = options.ss_match;

    bool end = false;
    Molecules mols(molecules);
    do {
      mols = molecules;
      setup(molecules);
      end  = monte_carlo(mols);
    }while(!end);

    cout<<"Run done"<<endl;
    cleanup();
    assign(mols,coor_best);
    {
      cout<<"BEST  energy "<<best_ene.energy<<"\n";
      cout<<best_ene<<"\n";
      for(int i=0;i<prot_buffer.size();i++){
	cout<<prot_buffer[i]<<" ";
	if((i+1)%5==0) cout<<"\n";
      }
      cout<<"\n";
      dump(mols,secstruct_best,best_ene.energy,best_pdb,true);
    }
    if(dcd)
      delete dcd;
    dcd = 0;

    //    if(0)
    {
      delete f3;
      delete f9;
      //      delete lf9;
      if(f15)
	delete f15;
    }
    return mols;
  }

  //   Molecules RoseMC::cont(const Molecules & molecules, const MDB & mdb, 
  // 			 const RoseOpts options, int sc)
  //   {
  //     Molecules mols(molecules);
  //     if(move_set_.size()==0)
  //     {
  //       add_backbone(0.01,4,123);
  //       add_backbone(0.02,4,12223);
  //       if(sc!=0){
  // 	add_sidechain(0.01,4,219375);
  // 	add_sidechain(0.02,4,2123245);
  //       }
  //       if(sc>1)
  // 	add_scramble(0.1,1,1213334);
  //     }
  //     srand48(options.seed);
  //     bond_rot = new BondRotation<double>(molecules);
  //     steps = options.steps;
  //     nprint = options.nprint;
  //     ndump = options.ndump;
  //     p_plain = options.p_plain;
  //     p_3 = options.p_3;
  //     begin_temp = options.begin_temp;
  //     end_temp = options.end_temp;
  //     sa_steps = options.sa_steps;
  //     pdb = options.pdb;
  //     best_pdb = options.best_pdb;
  //     setup(molecules,true);
  //     setup_rosemoves(molecules,options);
  //     prot.resize(molecules.fragment_size());
  //     prot_buffer.resize(molecules.fragment_size());


  //     monte_carlo(mols);

  //     cout<<"Run done"<<endl;
  //     cleanup();
  //     assign(mols,coor_best);
  //     {
  //       cout<<"BEST  energy "<<best_energy<<"\n";
  //       for(int i=0;i<prot_buffer.size();i++){
  // 	cout<<prot_buffer[i]<<" ";
  // 	if((i+1)%5==0) cout<<"\n";
  //       }
  //       cout<<"\n";
  //       dump(mols,secstruct_best,best_energy,best_pdb);
  //     }
  //     {
  //       delete f3;
  //       delete f9;
  //       //      delete lf9;
  //     }
  //     return mols;
  //   }

  //   Molecules RoseMC::run_rmsd(const Molecules & molecules, const MDB & mdb, 
  // 			     const RoseOpts options, int sc,
  // 			     const Molecules & ref_mol
  // 			     )
  //   {
  //     Molecules mols(molecules);
  //     {
  //       add_backbone(0.01,4,123);
  //       add_backbone(0.02,4,12223);
  //       if(sc!=0){
  // 	add_sidechain(0.01,4,219375);
  // 	add_sidechain(0.02,4,2123245);
  //       }
  //       if(sc>1)
  // 	add_scramble(0.1,1,1213334);
  //     }
  //     srand48(options.seed);
  //     bond_rot = new BondRotation<double>(molecules);
  //     steps = options.steps;
  //     nprint = options.nprint;
  //     ndump = options.ndump;
  //     p_plain = options.p_plain;
  //     p_3 = options.p_3;
  //     begin_temp = options.begin_temp;
  //     end_temp = options.end_temp;
  //     sa_steps = options.sa_steps;
  //     pdb = options.pdb;
  //     setup(molecules);
  //     setup_rosemoves(molecules,options);
  //     prot.resize(molecules.fragment_size());
  //     prot_buffer.resize(molecules.fragment_size());

  //     monte_carlo(mols,ref_mol);
  //     cout<<"Run done"<<endl;
  //     cleanup();
  //     assign(mols,coor_best);
  //     cout<<"BEST  rmsd "<<best_energy<<"\n";
  //     for(int i=0;i<prot_buffer.size();i++){
  //       cout<<prot_buffer[i]<<" ";
  //       if((i+1)%5==0) cout<<"\n";
  //     }
  //     cout<<"\n";
  //     {
  //       delete f3;
  //       delete f9;
  //       //      delete lf9;
  //     }
  //     return mols;
  //   }


  //   Molecules RoseMC::run_ref(const Molecules & molecules, const MDB & mdb, 
  // 			    const RoseOpts options, int sc, string sh_file,
  // 			    string ss, double t)
  //   {
  //     Molecules mols(molecules);
  //     {
  //       add_backbone(0.1,4,123);
  //       add_backbone(0.2,4,12223);
  //       if(sc!=0){
  // 	add_sidechain(0.1,4,219375);
  // 	add_sidechain(0.1,4,2123245);
  //       }
  //       if(sc>1)
  // 	add_scramble(0.1,1,1213334);
  //     }
  //     srand48(options.seed);
  //     bond_rot = new BondRotation<double>(molecules);
  //     steps = options.steps;
  //     nprint = options.nprint;
  //     ndump = options.ndump;
  //     p_plain = options.p_plain;
  //     p_3 = options.p_3;
  //     begin_temp = options.begin_temp;
  //     end_temp = options.end_temp;
  //     sa_steps = options.sa_steps;
  //     pdb = options.pdb;
  //     vdw_full = options.vdw_full;
  //     setup(molecules);
  //     setup_rosemoves(molecules,options);
  //     prot.resize(molecules.fragment_size());
  //     prot_buffer.resize(molecules.fragment_size());

  //     monte_carlo_shx(mols,sh_file,ss,t);
  //     cout<<"Run done"<<endl;
  //     cleanup();
  //     assign(mols,coor_best);
  //     cout<<"BEST  energy "<<best_energy<<"\n";
  //     {
  //       delete f3;
  //       delete f9;
  //       //      delete lf9;
  //     }
  //     return mols;
  //   }

  bool RoseMC::monte_carlo(Molecules & molecules){
    //init
    cout<<"DEBUG init\n";
    //    int init = 0;
    //    double ssfrac=0;
    while(1)
      {
	//      int acc_counter = 0;
	int rej = 0;
	int acc = 0;
	int rose = 0;
	int acc_rose = 0;
	bool was_rose = false;
	//      double dT = begin_temp;

	init_energy(current_ene);
	ene = current_ene;
	//       double energy = current_energy;
	//       {
	// 	ene.vdw_energy = current_ene.vdw_energy;
	// 	ene.ss_energy = current_ss_energy;
	// 	ene.cc_energy = current_cc_energy;
	// 	ene.hb_energy = current_hb_energy;
	// 	ene.bhb_energy = current_bhb_energy;
	// 	ene.apmf_energy = current_apmf_energy;
	// 	ene.dpmf_energy = current_dpmf_energy;
	// 	ene.env_energy = current_env_energy;
	// 	ene.rgyr_energy = current_rgyr_energy;
	//       }

	int best_counter = 0;
	best_ene = ene;
	coor_best.assign(coor_buffer);
	double temp = 2*begin_temp;
	for(int i=0;i<2*steps;i++){	
	  //RoseMove
	  was_rose = true;
	  int rose_move = 0;
	  while(rose_move<=100){
	    ++rose_move;
	    ++rose;
	    if(i<steps){
	      double q = drand48();
	      // 	    if(q<0.99)
	      if(!sc_){
		f9->apply(coor,secstruct,prot);
	      }
	      else {
		if(q<0.25)
		  apply(coor);
		else{
		  double qq = drand48();
		  if(qq<0.5)
		    f9->apply(coor,secstruct,prot);
		  else{
		    if(f15){
		      f15->apply(coor,secstruct,prot);
		    }
		    else f9->apply(coor,secstruct,prot);
		  }
		}
	      }
	    }
	    else{
	      double q = drand48();
	      if(q<0.5)
		f3->apply_turn(coor,coor_buffer,secstruct,prot);
	      else if(q<0.75)
		f9->apply_turn(coor,coor_buffer,secstruct,prot);
	      else{
		if(!sc_){
		  double qq = drand48();
		  if(qq<0.5)
		    f9->apply(coor,secstruct,prot);
		  else{
		    if(f15)
		      f15->apply(coor,secstruct,prot);
		    else f9->apply(coor,secstruct,prot);
		  }
		}
		else 
		  apply(coor);
	      }
	    }
	    if(rose_move==100) break;
	    if(!check_ss()){
	      secstruct = secstruct_buffer;
	      coor.assign(coor_buffer);
	    }
	    else break;
	  }
	
	  //Energy
	  init_energy(current_ene);
	  if(rose_metropolis(current_ene.energy,ene.energy,temp)){
	    //Accept
	    if(was_rose) ++acc_rose;
	    coor_buffer.assign(coor);
	    if(dcd)
	      dcd->write(coor_buffer);
	    secstruct_buffer = secstruct; 
	    prot_buffer = prot;
	    ene = current_ene;
	    rej = 0;
	    ++acc;
	    if(ene.energy<best_ene.energy){
	      best_ene = ene;
	      cout<<"BEST  energy "<<best_ene.energy<<" "<<i+1<<"\n";
	      cout<<best_ene<<"\n";
	      coor_best.assign(coor_buffer);
	      secstruct_best = secstruct;	    
	      prot_best = prot_buffer;
	      best_counter = 0;
	    }
	    else ++best_counter;
	  }
	  else {
	    coor.assign(coor_buffer);
	    if(dcd)
	      dcd->write(coor_buffer);
	    secstruct = secstruct_buffer; 
	    prot = prot_buffer;
	    ++rej;
	    acc = 0;
	  }

	  if((i%nprint)==0){
	    assign(molecules,coor_buffer);
	    cout<<"Status  "<<i+1<<" "<<temp<<"\n";
	    cout<<(double)acc_rose/(double)rose<<" "<<ene.energy<<"\n";
	    cout<<ene<<"\n";;
	    print();
	  }
	  if((i%ndump)==0) dump(molecules,secstruct_buffer,ene.energy,pdb);
	}




	assign(molecules,coor_best);
	coor_buffer.assign(coor_best);
	coor.assign(coor_best);
	secstruct_buffer = secstruct_best;
	secstruct = secstruct_best; 
	prot_buffer = prot_best;
	prot = prot_buffer;
	//Check virtual sec and sec
	{
	  string sec = secstruc_hbonds (molecules.protein(0));
	  int c=0;
	  for(int i=0;i<sec.size();i++){
	    char c1 = sec[i];
	    char c2 = secstruct_buffer[i];
	    if(c1=='E'||c1=='B') c1 = 'E';
	    else if(c1=='H'||c1=='G') c1 = 'H';
	    else c1='C';
	    if(c2=='E'||c2=='B') c2 = 'E';
	    else if(c2=='H'||c2=='G') c2 = 'H';
	    else c2='C';
	    if(c1==c2) ++c;
	  }
	  if(((double)c)/sec.size()>ss_match){
	    cout<<"Equal "<<((double)c)/sec.size()<<" "<<ss_match<<endl;
	    break;
	  }
	  else {
	    cout<<"SEC NOT EQ "<<((double)c)/sec.size()<<" "<<ss_match<<endl;
	    return false;
	    //break;
	  }

	}


      
      }


    {
      rose_ff->disable_hb();
      ene.hb_energy = current_ene.hb_energy = best_ene.hb_energy = 0;
    }
    cout<<"DEBUG init done\n";
    for(int KK=0;KK<sa_cycles;KK++)
      {
	{
	  assign(molecules,coor_best);
	  char b[16];
	  sprintf(b,"%i",KK);
	  dump(molecules,secstruct_best,best_ene.energy,best_pdb+string(b));
	  cout<<"DUMP "<<best_pdb+string(b)<<endl;
	}
	int rej = 0;
	int acc = 0;
	int rose = 0;
	int acc_rose = 0;
	bool was_rose = false;
	rose_energy(ene,4);//From almost-0.99.2
	current_ene = ene;
	int best_counter = 0;
	best_ene = ene;

	coor_best.assign(coor_buffer);
	double dT = (begin_temp-end_temp)/sa_steps;
	for(int s=0;s<sa_steps;s++){
	  rej = 0;
	  acc =0;
	  rose = 0;
	  acc_rose = 0;
	  double temp = begin_temp-dT*s;
	  for(int i=0;i<steps;i++){
	    double move_switch = drand48();
	    //	if((s>(sa_steps/2))&&(move_switch<p_plain)){
	    if((move_switch<p_plain)){
	      //Plain move
	      apply(coor);
	      was_rose = false;
	    }
	    else {
	      //RoseMove
	      was_rose = true;
	      int rose_move = 0;
	      while(rose_move<=100){
		++rose_move;
		++rose;
		double sw_39 =   drand48();
		if(sw_39<p_3){
		  //	      f3->apply_turn(coor,secstruct,prot);
		  f3->apply_turn(coor,coor_buffer,secstruct,prot);
		}
		else {
		  f3->apply(coor,secstruct,prot);
		}
		if(rose_move==100) break;
		if(!check_ss()){
		  secstruct = secstruct_buffer;
		  coor.assign(coor_buffer);
		}
		else break;
	      }
	    }
	    //Energy
	    rose_energy(current_ene,4);//From almost-0.99.2

	
	    if(rose_metropolis(current_ene.energy,ene.energy,temp)){
	      //Accept
	      if(was_rose) ++acc_rose;
	      coor_buffer.assign(coor);
	      if(dcd)
		dcd->write(coor_buffer);
	      secstruct_buffer = secstruct; 
	      prot_buffer = prot;
	      ene = current_ene;
	      rej = 0;
	      ++acc;
	      if(ene.energy<best_ene.energy){	    
		best_ene = ene;
		cout<<"BEST  energy "<<best_ene.energy<<" "<<i+1<<"\n";
		cout<<best_ene<<"\n";
		coor_best.assign(coor_buffer);
		secstruct_best = secstruct;	    
		prot_best = prot_buffer;
		best_counter = 0;
	      }
	      else ++best_counter;
	    }
	    else {
	      coor.assign(coor_buffer);
	      if(dcd)
		dcd->write(coor_buffer);
	      secstruct = secstruct_buffer; 
	      prot = prot_buffer;
	      ++rej;
	    }
	    if((i%nprint)==0){
	      assign(molecules,coor_buffer);
	      cout<<"almost@rose>  "<<i+1<<" "<<temp<<"\n";
	      // 	      <<(double)acc_rose/(double)rose<<" "
	      // 	      <<(double)acc/(i+1)<<"\n";
	      cout<<ene<<"\n";
	      cout<<secstruct_buffer<<"\n"
		  <<secstruc_hbonds (molecules.protein(0))<<"\n\n";
	      print();
	    }
	    if((i%ndump)==0) dump(molecules,secstruct_buffer,ene.energy,pdb);
	  }
	}
      }
    return true;
  }


  void RoseMC::monte_carlo(Molecules & molecules, const Molecules & ref_mols){
    //     NameSelection s1(molecules,"///CA");
    //     NameSelection s2(ref_mols,"///CA");
    //     Coor<double> ref_coor(ref_mols);
    //     double r[3][3];
    //     int rej = 0;
    //     int acc = 0;
    //     int rose = 0;
    //     int acc_rose = 0;
    //     bool was_rose = false;
    //     double energy = rose_energy()+20*rmsd2_mat(coor,ref_coor,s1.result(),s2.result(),r);
    //     double current_energy = energy;
    //     int best_counter = 0;
    //     best_energy = 10000;
    //     coor_best.assign(coor_buffer);
    //     double dT = (begin_temp-end_temp)/sa_steps;
    //     for(int s=0;s<sa_steps;s++){
 
    //       double temp = begin_temp-dT*s;
    //       for(int i=0;i<steps;i++){
    // 	double move_switch = drand48();
    // 	if((move_switch<p_plain)){
    // 	  //Plain move
    // 	  apply(coor);
    // 	  was_rose = false;
    // 	}
    // 	else {
    // 	  //RoseMove
    // 	  was_rose = true;
    // 	  int rose_move = 0;
    // 	  while(rose_move<=100){
    // 	    ++rose_move;
    // 	    ++rose;
    // 	    double sw_39 =   drand48();
    // 	    if(sw_39<p_3){
    // 	      f3->apply(coor,secstruct,prot);
    // 	    }
    // 	    else {
    // 	      f9->apply(coor,secstruct,prot);
    // 	    }
    // 	    if(rose_move==100) break;
    // 	    if(!check_ss()){
    // 	      secstruct = secstruct_buffer;
    // 	      coor.assign(coor_buffer);
    // 	    }
    // 	    else break;
    // 	  }
    // 	}
    // 	//Energy
    // 	double rmsd = rmsd2_mat(coor,ref_coor,s1.result(),s2.result(),r);
    // 	current_energy = rose_energy()+20*rmsd;

    // 	if(rose_metropolis(current_energy,energy,temp)//&&
    // 	   //(rmsd<best_energy)
    // 	   ){
    // 	  //Accept
    // 	  if(was_rose) ++acc_rose;
    // 	  coor_buffer.assign(coor);
    // 	  secstruct_buffer = secstruct; 
    // 	  prot_buffer = prot;
    // 	  energy = current_energy;
    // 	  rej = 0;
    // 	  ++acc;
    // 	  //cout<<"DEB HB "<<hb_energy<<"\n";
	  
    // 	  if(rmsd<best_energy){
    // 	    best_energy = rmsd;
    // 	    cout<<"BEST  rmsd "<<best_energy<<"\n";
    // 	    coor_best.assign(coor_buffer);
    // 	    secstruct_best = secstruct;	    
    // 	    prot_best = prot_buffer;
    // 	    best_counter = 0;
    // 	  }
    // 	  else ++best_counter;
    // 	  // 	  cout<<"Debug Acc "<<i+1<<" "<<temp<<" "<<energy<<"\n"
    // 	  // 	      <<secstruct_buffer<<"\n\n";
    // 	}
    // 	else {
    // 	  coor.assign(coor_buffer);
    // 	  secstruct = secstruct_buffer; 
    // 	  prot = prot_buffer;
    // 	  // 	  cout<<"Debug Rej "<<i+1<<" "<<temp<<" "<<energy<<"\n"
    // 	  // 	      <<secstruct_buffer<<"\n\n";	
    // 	  ++rej;
    // 	  acc = 0;
    // 	  // 	  if(best_counter>=1000){
    // 	  // 	    best_counter= 0;
    // 	  // 	    coor_buffer.assign(coor_best);
    // 	  // 	    secstruct_buffer = secstruct_best; 
    // 	  // 	    energy = best_energy;
    // 	  // 	    coor.assign(coor_best);
    // 	  // 	    secstruct = secstruct_best;
    // 	  // 	    current_energy = best_energy;
    // 	  // 	    cout<<"Resetting best"<<"\n";
    // 	  // 	  }
    // 	}

    // 	// 	if((p_plain!=1)&&((double)acc_rose/(double)rose<0.01)){
    // 	// 	  p_plain =1;
    // 	// 	  cout<<"p_plain  set to 1"<<endl;
    // 	// 	}
    // 	if((i%nprint)==0){
    // 	  cout<<"Debug  "<<i+1<<" "<<temp<<" "
    // 	      <<(double)acc_rose/(double)rose<<" "<<energy<<"\n"
    // 	      <<vdw_energy<<" "<<ss_energy<<" "<<" "<<cc_energy<<" "<<"\n"
    // 	      <<secstruct_buffer<<"\n"
    // 	      <<secstruc_hbonds (molecules.protein(0))<<"\n\n";
    // 	  print();
    // 	}
    // 	if((i%ndump)==0) dump(molecules,secstruct_buffer,energy,pdb);
    // 	// 	if(acc==1){ 
    // 	// 	  temp = 0.99*temp; acc=0;
    // 	// 	}
    // 	// 	if(rej==10){ 
    // 	// 	  temp = 1.01*temp; rej=0;
    // 	// 	}
    //       }
    //     }
  }


  void RoseMC::monte_carlo_shx(Molecules & molecules,string file,string ss,
			       double t){
    //     Molecules buff_mol(molecules);
    //     int rej = 0;
    //     int acc = 0;
    //     int rose = 0;
    //     int acc_rose = 0;
    //     bool was_rose = false;
    //     secstruct =ss;
    //     secstruct_buffer =ss;
    //     double energy = rose_energy_ref();
    //     double s = score_silent(molecules,file);
    //     s = max(s,t);
    //     energy = energy/log(1+0.01*s);
    //     double current_energy = energy;
    //     int best_counter = 0;
    //     best_energy = energy;
    //     coor_best.assign(coor_buffer);
    //     double dT = (begin_temp-end_temp)/sa_steps;
    //     for(int s=0;s<sa_steps;s++){
    //       double temp = begin_temp-dT*s;
    //       for(int i=0;i<steps;i++){
    // 	double move_switch = drand48();
    // 	//	if((s>(sa_steps/2))&&(move_switch<p_plain)){
    // 	if((move_switch<p_plain)){
    // 	  //Plain move
    // 	  apply(coor);
    // 	  was_rose = false;
    // 	}
    // 	else {
    // 	  //RoseMove
    // 	  was_rose = true;
    // 	  int rose_move = 0;
    // 	  while(rose_move<=100){
    // 	    ++rose_move;
    // 	    ++rose;
    // 	    double sw_39 =   drand48();
    // 	    if(sw_39<p_3){
    // 	      f3->apply(coor,secstruct,prot);
    // 	    }
    // 	    else {
    // 	      f9->apply(coor,secstruct,prot);
    // 	    }
    // 	    if(rose_move==100) break;
    // 	    if(!check_ss()){
    // 	      secstruct = secstruct_buffer;
    // 	      coor.assign(coor_buffer);
    // 	    }
    // 	    else break;
    // 	  }
    // 	}
    // 	//Energy
    // 	current_energy = rose_energy_ref();
    // 	assign(buff_mol,coor);
    // 	double sc = score_silent(buff_mol,file);
    // 	double sc_old = sc;
    // 	sc = max(sc,t);
    // 	current_energy = current_energy/log(1+0.01*sc);
    // 	if(rose_metropolis(current_energy,energy,temp)){
    // 	  //Accept
    // 	  if(was_rose) ++acc_rose;
    // 	  coor_buffer.assign(coor);
    // 	  secstruct_buffer = secstruct; 
    // 	  energy = current_energy;
    // 	  rej = 0;
    // 	  ++acc;
    // 	  //cout<<"DEB HB "<<hb_energy<<"\n";
    // 	  if(energy<best_energy){
    // 	    best_energy = energy;
    // 	    cout<<"BEST  energy "<<best_energy<<" "<<sc_old<<" "
    // 		<<calc_rgyr(coor)<<"\n";
    // 	    coor_best.assign(coor_buffer);
    // 	    secstruct_best = secstruct;
    // 	    best_counter = 0;
    // 	  }
    // 	  else ++best_counter;
    // 	  // 	  cout<<"Debug Acc "<<i+1<<" "<<temp<<" "<<energy<<"\n"
    // 	  // 	      <<secstruct_buffer<<"\n\n";
    // 	}
    // 	else {
    // 	  coor.assign(coor_buffer);
    // 	  secstruct = secstruct_buffer; 
    // 	  // 	  cout<<"Debug Rej "<<i+1<<" "<<temp<<" "<<energy<<"\n"
    // 	  // 	      <<secstruct_buffer<<"\n\n";	
    // 	  ++rej;
    // 	  acc = 0;
    // 	  // 	  if(best_counter>=1000){
    // 	  // 	    best_counter= 0;
    // 	  // 	    coor_buffer.assign(coor_best);
    // 	  // 	    secstruct_buffer = secstruct_best; 
    // 	  // 	    energy = best_energy;
    // 	  // 	    coor.assign(coor_best);
    // 	  // 	    secstruct = secstruct_best;
    // 	  // 	    current_energy = best_energy;
    // 	  // 	    cout<<"Resetting best"<<"\n";
    // 	  // 	  }
    // 	}
    // 	if((i%nprint)==0){
    // 	  cout<<"Debug  "<<i+1<<" "<<temp<<" "
    // 	      <<(double)acc_rose/(double)rose<<" "<<energy<<"\n"
    // 	      <<secstruct_buffer<<"\n\n";
    // 	  print();
    // 	}
    // 	if((i%ndump)==0) dump(molecules,secstruct_buffer,energy,pdb);
    // 	// 	if(acc==1){ 
    // 	// 	  temp = 0.99*temp; acc=0;
    // 	// 	}
    // 	// 	if(rej==10){ 
    // 	// 	  temp = 1.01*temp; rej=0;
    // 	// 	}
    //       }
    //     }
  }



  void RoseMC::setup(const Molecules & molecules,bool cont){
    if(!cont){
      coor.copy(molecules);
      coor_buffer.copy(molecules);
      coor_best.copy(molecules);
      secstruct.resize(molecules.fragment_size(),'C');
      secstruct = secstruc_hbonds (molecules.protein(0));
      for(int i=0;i<secstruct.size();i++){
	if(secstruct[i]==' ') secstruct[i] = 'C';
      }
      secstruct_buffer = secstruct;//.resize(molecules.fragment_size(),'C');
    }
    else{
      coor.copy(molecules);
      coor_buffer.copy(molecules);
      coor_best.copy(molecules);
      secstruct_buffer = secstruct = secstruct_best;
    }
  }


  void RoseMC::setup_rosemoves(const Molecules & molecules,
			       const RoseOpts &options){
    f3 = new Fragment<double>(molecules,*bond_rot,
			      options.seed3,
			      options.frg3);

    f9 = new Fragment9<double>(molecules,*bond_rot,
			       options.seed9,
			       options.frg9);

    if(options.frg15!="none")
      f15 = new Fragment15<double>(molecules,*bond_rot,
				   options.seed9,
				   options.frg15);
    //     lf9 = new LocalFragment<double,9>(molecules,*bond_rot,
    // 				      options.seed9,
    // 				      options.frg9);

  }

  double RoseMC::init_energy(RoseMCEnergy & mc_energy){
    mc_energy.clear();
    //
    mc_energy.vdw_energy = max(0.0,rose_ff->ff_energy(coor));

    mc_energy.cc_energy = 0;
    for(int i=0;i<cc.size();i++){
      mc_energy.cc_energy += cc.constraints()[i]->energy(coor);
    }
    mc_energy.hb_energy = rose_ff->hb_energy(coor,secstruct);

    mc_energy.ss_energy = 4*rose_ff->ss_energy(coor,secstruct);
    mc_energy.summ();
    return mc_energy.energy;
  }


  double RoseMC::rose_energy(RoseMCEnergy & mc_energy,double k){
    mc_energy.clear();
    mc_energy.vdw_energy = max(0.0,rose_ff->ff_energy(coor));

    //APMF
    mc_energy.apmf_energy = rose_ff->apmf_energy(coor);

    mc_energy.cc_energy = 0;
    for(int i=0;i<cc.size();i++){
      mc_energy.cc_energy += cc.constraints()[i]->energy(coor);
    }

    mc_energy.hb_energy = rose_ff->hb_energy(coor,secstruct);

    mc_energy.bhb_energy = rose_ff->bhb_energy(coor,secstruct);
    mc_energy.rgyr_energy = rose_ff->rgyr_energy(coor);
    mc_energy.ss_energy = k*rose_ff->ss_energy(coor,secstruct);
    mc_energy.summ();
    return mc_energy.energy;
  }

  double RoseMC::rose_energy_ref(RoseMCEnergy & mc_energy){
    mc_energy.clear();
    mc_energy.vdw_energy = rose_ff->ff_energy(coor);

    //APMF
    mc_energy.apmf_energy = rose_ff->apmf_energy(coor);

    mc_energy.cc_energy = 0;
    for(int i=0;i<cc.size();i++){
      mc_energy.cc_energy += cc.constraints()[i]->energy(coor);
    }

    mc_energy.hb_energy = rose_ff->hb_energy(coor,secstruct);

    mc_energy.bhb_energy = 0;
    mc_energy.bhb_energy = rose_ff->bhb_energy(coor,secstruct);
    mc_energy.rgyr_energy = rose_ff->rgyr_energy(coor);

    mc_energy.summ();


    return mc_energy.energy;
  }

}

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