/***************************************************************************
                          solvent.cpp  -  description
                             -------------------
    begin                : Mon Jan 27 06:38:48 2003
    copyright            : (C) 2002 by Cavalli Andrea
    email                : cavalli@bioc.unizh.ch
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <almostexception.h>
#include <almosttypes.h>
#include <almosttemplate.h>
#include <pdb.h>
#include <mdb/mdb.h>
#include <buildutils.h>
#include <molecules/solvent.h>

void Almost::Solvent::build(PDBSegment & segment,MDB & mdb,
			    string fpatch, string lpatch) 
{
  try{
    //Step 1: Find sequence from PDB
    vector<FragmentData> data;
    {
      find_solvent_sequence(segment,mdb,data,fpatch,lpatch);
    }
    //Step 1: Done

    //Step 2
    //Add fragments
    {
      add_segment_fragments(data);
    }
    //Step 2: Done

    //Step 3
    //Add atoms
    {
      add_segment_atoms(segment,mdb,data);
    }
    //Step 3: Done


    //Step 4
    //Add fixed exclusions
    {
      add_fixed_exclusions(mdb,data);
    }
    //Step 4: Done

    //Step 5
    // Check heavy atoms
    {
      try{
	check_heavy_atoms(mdb,fpatch,lpatch);
      }
      catch(AlmostException exc){
	aout<<exc.error()<<endl;
	throw;
      }
      catch(...){
	aout<<"Exception caught. Check your PDB file"<<endl;
	exit(-1);
      }
    } 
    //Step 5: Done
    
    //Step 6
    //Check H-atoms
    {
      try{
	check_h_atoms(mdb,fpatch,lpatch);
      }
      catch(AlmostException exc){
	aout<<exc.error()<<endl;
	throw;
      }
      catch(...){
	aout<<"Exception caught. Check your PDB file"<<endl;
	exit(-1);
      }
    }
    //Step 6: Done
    
    //Step 7
    //Make bonds
    {
      try{
	make_bonds(mdb,data);
      }
      catch(AlmostException exc){
	aout<<exc.error()<<endl;
	throw;
      }
      catch(...){
	aout<<"Exception caught. Check your PDB file"<<endl;
	exit(-1);
      }
    }
    //Step 7: Done

    //Step 8
    //Add angles
    {
      add_angles(mdb,data);
    }
    //Step 8: Done

    //Step 9
    //Add dihedrals
    {
      add_dihedrals(mdb,data);
    }
    //Step 9: Done

    //Step 10
    //Add impropers
    {
      add_impropers(mdb,data);
    }
    //Step 10: Done

    //Step 11
    //Add groups
    {
      add_groups(mdb,data);
    }

    {
      setup_parameters(mdb);
    }

    //Step 12
    //Print info
    {
      build_info();
    }
    //Step 12: Done
  }
  catch(AlmostException e){
    throw;
    //AEX(e.error());
    //    throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,e.error());
  }
} 


void Almost::Solvent::build_sequ(const vector<string> &sequ, MDB & mdb,
  string fpatch, string lpatch){
  try{
    vector<FragmentData> data;
    {
      find_solvent_sequence(sequ,mdb,data,fpatch,lpatch);
    }

    {
      add_segment_fragments(data);
    }

    {
      add_segment_atoms(sequ,mdb,data);
    }
    {
      add_fixed_exclusions(mdb,data);
    }

    {
      try{
	make_bonds(mdb,data);
      }
      catch(AlmostException exc){
	aout<<exc.error()<<endl;
	throw;
      }
      catch(...){
	aout<<"Exception caught. Check your PDB file"<<endl;
	exit(-1);
      }
    }

    {
      add_angles(mdb,data);
    }

    {
      add_dihedrals(mdb,data);
    }

    {
      add_impropers(mdb,data);
    }

    {
      add_groups(mdb,data);
    }

    {
      generate_coor(mdb,data);
    }

    {
      setup_parameters(mdb);
    }


    {
      build_info();
    }

  }
  catch(AlmostException e){
    throw;
  }
}


void Almost::Solvent::add_segment_fragments(vector<FragmentData> & data){
  for(int i=0;i<data.size();i++){
    int fragd = add_fragment(data[i].fragment);
    if(fragd<0){
      aout<<"FATAL ERROR: can't add fragment "<<data[i].fragment<<" to segment"
	  <<endl;
      exit(-1);
    }
  }
}

void Almost::Solvent::add_segment_atoms( PDBSegment & segment,
					 MDB & mdb,
					 vector<FragmentData> & data){
  char buff[16];
  for(int i=0;i<data.size();i++){
    int resSeq = segment[i].fragSeq();
    string resName = segment[i].name();
    sprintf(buff,"_%i",resSeq);
    resName +=buff;
    int fragd = find_fragment_by_name(resName);
    if(fragd<0){
      cerr<<"No such fragment "<<resName<<endl;
      exit(-1);
    }
    for(int j=0;j<segment[i].size();++j){
      TAtom tatom = mdb.atom(resName+"/"+segment[i][j].name,
			     data[i].fpatch,data[i].lpatch);
      AtomT atomt = mdb.atomT(tatom.type_name());
      Atom atom(segment[i][j].name,tatom.type_name(),
		segment[i][j].x,segment[i][j].y,segment[i][j].z,
		atomt.mass(),
		tatom.charge());
      int atomd = add_atom(atom,fragd);
      if(atomd<0){
	aout<<"FATAL ERROR: couldn't add atom "<<resName+"/"+segment[i][j].name
	    <<" to segment"<<endl;
	exit(-1);
      }
    }
  }
}


void Almost::Solvent::add_segment_atoms( const vector<string> & sequ,
					 MDB & mdb,
					 vector<FragmentData> & data){
  //  char buff[16];
  for(int i=0;i<data.size();i++){
    //    int resSeq = i;
    string resName = data[i].fragment;
    int fragd = find_fragment_by_name(resName);
    if(fragd<0){
      cerr<<"No such fragment "<<resName<<endl;
      exit(-1);
    }

    vector<string> atoms;
    string _fpatch;
    string _lpatch;
    _fpatch="NONE";
    _lpatch="NONE";
    if(data[i].fpatch=="DEFAULT"){
      _fpatch = mdb.fragment(resName).def_first_patch();
    }else{
      _fpatch=data[i].fpatch;
    }
    if(data[i].lpatch=="DEFAULT"){
      _lpatch = mdb.fragment(resName).def_last_patch();
    }else{
      _lpatch=data[i].lpatch;
    }
    atoms = mdb.atoms(fragment_[i],
		      MDB::ALL_ATOMS,
		      _fpatch,_lpatch);  
    for(int j=0;j<atoms.size();++j){
      TAtom tatom = mdb.atom(resName+"/"+atoms[j],
			     data[i].fpatch,data[i].lpatch);
      AtomT atomt = mdb.atomT(tatom.type_name());
      Atom atom(atoms[j],tatom.type_name(),
 		9999,9999,9999,
		atomt.mass(),
 		tatom.charge());
      int atomd = add_atom(atom,fragd);
      if(atomd<0){
	aout<<"FATAL ERROR: couldn't add atom "<<resName+"/"+atoms[j]
	    <<" to segment"<<endl;
	exit(-1);
      }
    }
  }
}

void Almost::Solvent::add_fixed_exclusions(MDB & mdb,
					   vector<FragmentData> & data){
  
  for(int i=0;i<atom_.size();i++){
    string name = atom_name(i);
    int fd = find_fragment(i);
    if(fd<0){
      aout<<"FATAL ERROR: couldn't find fragment of atom "<<name<<endl;
      exit(-1);
    }
    //    string fname = fragment_name(fd,Solvent::SHORT);
    set<string> excl = mdb.atom(name,
				data[fd].fpatch,
				data[fd].lpatch).exclusion_list();
    set<string>::iterator iter = excl.begin();
    while(iter!=excl.end()){
      int indx = find_atom(data[fd].fragment+string("/")+(*iter));
      if(indx<0){
	aout<<"FATAL ERROR: couldn't find atom "
	    <<data[fd].fragment+string("/")+(*iter)<<endl;
	exit(-1);
      }
      add_excl(i,indx);
      ++iter;
    }
  }
}

void Almost::Solvent::check_heavy_atoms(MDB & mdb,
					string fpatch, 
					string lpatch){
  string _fpatch;
  string _lpatch;
  vector<string> heavy_atoms;
  for(int i=0;i<fragment_.size();i++){
    _fpatch="NONE";
    _lpatch="NONE";
    if(i==0){ 
      if(fpatch=="DEFAULT"){
	_fpatch = mdb.fragment(fragment_[i]).def_first_patch();
      }else{
	_fpatch=fpatch;
      }
    } 
    if(i==(fragment_.size()-1)){
      if(lpatch=="DEFAULT"){
	_lpatch = mdb.fragment(fragment_[i]).def_last_patch();
      }else{
	_lpatch=lpatch;
      }
    }       
      
    heavy_atoms = mdb.atoms(fragment_[i],
			    MDB::HEAVY_ATOMS,
			    _fpatch,_lpatch);
    // 	aout <<"Fragment "<<_fragment[i].name()<< " patch "<<_fpatch<<" "
    // 	     <<_lpatch<<endl;
    int atomcount =0; 
    for(int j = 0; j< heavy_atoms.size();j++){
      int atomd = find_atom(fragment_[i]+"/"+heavy_atoms[j]);
      if(atomd!=-1) 
	++atomcount;
      else aout <<fragment_[i]+"/"+heavy_atoms[j]<< " not found"<<endl;
    }
    if(atomcount!= heavy_atoms.size()){
      aout <<"error in frag "<<i<<" atoms "<<atomcount<<endl;
      aout <<"atomcount "<<atomcount<<" heavy_atoms.size() "<<heavy_atoms.size()<<endl;
      throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,
			    "missing heavy atom in fragment "+fragment_[i]);
    }
  }
}


void Almost::Solvent::check_h_atoms(MDB & mdb,
				    string fpatch, 
				    string lpatch){
  string _fpatch;
  string _lpatch;
  vector<string> heavy_atoms;
  for(int i=0;i<fragment_.size();i++){
    _fpatch="NONE";
    _lpatch="NONE";
    if(i==0){ 
      if(fpatch=="DEFAULT"){
	_fpatch = mdb.fragment(fragment_[i]).def_first_patch();
      }else{
	_fpatch=fpatch;
      }
    } 
    if(i==(fragment_.size()-1)){
      if(lpatch=="DEFAULT"){
	_lpatch = mdb.fragment(fragment_[i]).def_last_patch();
      }else{
	_lpatch=lpatch;
      }
    } 
    heavy_atoms = mdb.atoms(fragment_[i],
			    MDB::HEAVY_ATOMS,
			    _fpatch,_lpatch);
    for(int j = 0; j< heavy_atoms.size();j++){
      vector<string> h_atoms = 
	mdb.bonds(fragment_[i]+"/"+heavy_atoms[j],
		  MDB::H_BONDS,_fpatch,_lpatch);
      for(int k=0;k<h_atoms.size();k++){
	int atomd = find_atom(fragment_[i]+"/"+h_atoms[k]);
	if(atomd<=-1){
	  // 	      aout <<atomd<<endl;
	  throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,
				"Hydrogen "+fragment_[i]+"/"
				+h_atoms[k]+
				" missing for atom "
				+ fragment_[i]+"/"+heavy_atoms[j]+
				" in fragment "+fragment_[i]);
	}
      }
    }
  }
}



void Almost::Solvent::make_bonds(MDB & mdb,
				 vector<FragmentData> & data){
  for(int i=0;i<data.size();i++){
    vector<string> atoms = mdb.atoms(fragment_[i],
				     MDB::ALL_ATOMS,
				     data[i].fpatch,
				     data[i].lpatch);
    for(int j = 0; j< atoms.size();j++){
      vector<string> bond_atoms = 
	mdb.bonds(fragment_[i]+"/"+atoms[j],
		  MDB::OUT_BONDS,
		  data[i].fpatch,data[i].lpatch);

      vector<int> bond_kind = 
	mdb.bonds_kind(fragment_[i]+"/"+atoms[j],
		       MDB::OUT_BONDS,
		       data[i].fpatch,data[i].lpatch);

      int atomd1 = find_atom(fragment_[i]+"/"+ atoms[j]);
      if(atomd1<=-1){
	throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,
			      "Atom "+fragment_[i]+"/"+ atoms[j]+
			      " not found in fragment "+
			      fragment_[i]);
      }
      for(int k=0;k<bond_atoms.size();++k){
	string search_name;
	if(bond_atoms[k][0]=='+'){
	  if(i!=(fragment_.size()-1))
	    search_name =fragment_[i+1]+"/"+string(bond_atoms[k].begin()+1,
							  bond_atoms[k].end());
	  else continue;
	}
	else search_name =fragment_[i]+"/"+bond_atoms[k];
	int atomd2 = find_atom(search_name);
	if(atomd2<0){
	  throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,
				"Atom "+search_name+
				" not found in fragment "+
				fragment_[i]);
	}
	if(!add_bond(atomd1,atomd2,bond_kind[k])){
	  aout<<"Warning: bond beteewn "<<atomd1+1<<" "<<atomd2+1<<" not added!"<<endl;
	}
      }
    }
  } 
}



int Almost::Solvent::find_atom(string name){
  string frag;
  string atom;
  int pos=name.find("/");
  if(pos==0) return -1;
  if(pos==string::npos) return -1;
  frag = string(name.begin(),name.begin()+pos);
  atom = string(name.begin()+pos+1 ,name.end());
  int nfrag = find_fragment_by_name(frag);
  if(nfrag==-1) return -1;
  for(int i=0;i<fragment_atoms_[nfrag].size();++i)
    if(atom_[fragment_atoms_[nfrag][i]].name()==atom) return fragment_atoms_[nfrag][i];
  return -1;
}

void Almost::Solvent::add_angles(MDB & mdb,vector<FragmentData> & data){
  logMsg(Almost::SolventPr,"No angle generation for solvent");
  return;
  angle_generate();
}
    
void Almost::Solvent::add_dihedrals(MDB & mdb,vector<FragmentData> & data){
  logMsg(Almost::SolventPr,"No dihedral generation for solvent");
  return;
  //auto generate
  //WHAT ???
  for(int i=0;i<fragment_.size();i++){
    vector<TDihedral> dihe;
    
    dihe = mdb.dihedrals(fragment_[i],
			 data[i].fpatch,
			 data[i].lpatch);

    //register
    for(int j=0;j<dihe.size();j++){
      if((i==0)&&is_backward(dihe[j])) continue;
      if((i==fragment_.size()-1)&&is_forward(dihe[j])) continue;


      //names
      string name[4];
      int index[4];
      //      int frag[4];
      for(int k=0;k<4;k++){
	name[k] = dihe[j].atom(k);
	if(name[k][0]=='+'){
	  name[k] = 
	    fragment_[i+1]+"/"+string(name[k].begin()+1,name[k].end());
	}
	else if(name[k][0]=='-'){
	  name[k] = 
	    fragment_[i-1]+"/"+string(name[k].begin()+1,name[k].end());
	}  
	else {
	  name[k] = fragment_[i]+"/"+name[k];
	}
      }
   
      for(int k=0;k<4;k++) index[k] = find_atom(name[k]);

      //register
      vector<int> dihe_;
      for(int k=0;k<4;k++) {
	if(index[k]==-1){
	  aout<<"Couldn't find atom "<< k+1<< " of dihedral "<<j+1<<endl;
	  exit(-1);
	}
	dihe_.push_back(index[k]);
      }
      if(!dihedral_.add_dihedral(dihe_)) 
	aout<<"Warning: cant add dihedral"<<endl;
    }
  }
  if(dihedral_.size()==0){  
    //    aout<<"No dihedral information in MDB assuming (lazy) autogeneration ....\n"; 
    dihedral_generate();
  }
}

void Almost::Solvent::add_impropers(MDB & mdb,vector<FragmentData> & data){
  for(int i=0;i<fragment_.size();i++){
    vector<TImproper> imph;
    
    imph = mdb.impropers(fragment_[i],
			 data[i].fpatch,
			 data[i].lpatch);

    //register
    for(int j=0;j<imph.size();j++){
      if((i==0)&&is_backward(imph[j])) continue;
      if((i==fragment_.size()-1)&&is_forward(imph[j])) continue;


      //names
      string name[4];
      int index[4];
      //      int frag[4];
      for(int k=0;k<4;k++){
	name[k] = imph[j].atom(k);
	if(name[k][0]=='+'){
	  name[k] = 
	    fragment_[i+1]+"/"+string(name[k].begin()+1,name[k].end());
	}
	else if(name[k][0]=='-'){
	  name[k] = 
	    fragment_[i-1]+"/"+string(name[k].begin()+1,name[k].end());
	}  
	else {
	  name[k] = fragment_[i]+"/"+name[k];
	}
      }
   
      for(int k=0;k<4;k++) index[k] = find_atom(name[k]);

      //register
      vector<int> imph_;
      for(int k=0;k<4;k++) {
	if(index[k]==-1){
	  aout<<"Couldn't find atom "<< k+1<< " of improper "<<j+1<<endl;
	  exit(-1);
	}
	imph_.push_back(index[k]);
      }
      if(!improper_.add_improper(imph_))
	aout<<"Warning can't add improper"<<endl;
    }
  }
}

void Almost::Solvent::add_groups(MDB & mdb,vector<FragmentData> & data){
  for(int i=0;i<fragment_.size();i++){
    vector<TGroup> group;
    group = mdb.groups(fragment_[i],
		       data[i].fpatch,
		       data[i].lpatch);
    //register
    for(int j=0;j<group.size();j++){
      vector<int > tmp;
      set<string>::const_iterator iter = group[j].atoms().begin(),
	end = group[j].atoms().end();
      while(iter!=end){
	string  name = *iter;
	++iter;
	int index;
	//	int frag;
	name  = fragment_[i]+"/"+name;
	index = find_atom(name);
	tmp.push_back(index);
      }
      sort(tmp.begin(),tmp.end());
      group_.push_back(tmp);
    }
  }
}


void Almost::Solvent::generate_coor(MDB & mdb,
				    vector<FragmentData> & data){


  ICManager ic_man;
  
  for(int i=0;i<fragment_.size();i++){
    vector<TIC> ic = mdb.ic(data[i].fragment,
			    data[i].fpatch,
			    data[i].lpatch);
    for(int j= 0;j<ic.size();j++){
      if((i==0)&&(ic_man.is_backward(ic[j]))) continue;
      if((i==data.size()-1)&&(ic_man.is_forward(ic[j]))) continue;

      string name[4];
      bool is_improper = false;
      int atomd[4];
      for(int k=0;k<4;k++){
	string atomk = ic[j].atom(k);
	if(atomk[0] == '+'){ 
	  name[k] = string(atomk.begin()+1,
			   atomk.end());
	  name[k] = data[i+1].fragment+"/"+name[k];
	}
	else if(atomk[0] == '-'){ 
	  name[k] = string(atomk.begin()+1,
			   atomk.end());
	  name[k] = data[i-1].fragment+"/"+name[k];
	}
	else if(atomk[0] == '*'){ 
	  name[k] = string(atomk.begin()+1,
			   atomk.end());
	  name[k] = data[i].fragment+"/"+name[k];
	  is_improper = true;
	}
	else{
	  name[k] = data[i].fragment+"/"+atomk;
	}
      }
      for(int k=0;k<4;k++){
	atomd[k] = find_atom(name[k]);
	if(atomd[k]<0){
	  aout<<"FATAL ERROR: "<<name[k]<<
	    " not found\n";
	  aout<<"Fragment "<<i<<" "<< data[i].fpatch
	      <<" "<< data[i].lpatch<<" "<<ic.size()<<"\n";
	  exit(-1);
	}
      }

      //add to ICManager
      ic_man.ic.push_back(ICManager::IC());
      for(int k=0;k<4;k++){
	ic_man.ic.back().atom[k] = name[k];
	ic_man.ic.back().atom_type[k] = atom_[atomd[k]].type_name();
	ic_man.ic.back().atomd[k] = atomd[k];
      }
      for(int k=0;k<5;k++) ic_man.ic.back().data[k] = ic[j].data(k);
      ic_man.ic.back().improper = is_improper;
    }
  }

  ic_man.sort();
  ic_man.fill(mdb);

  //now build
  
  //Step 1
  set<int> known;
  if((ic_man.ic.size()==0)&&(known.size()!=atom_.size())){
    aout<<"\t\t WARNING: Undefined atoms\n";
    return;
  }
  ICManager::IC & ic = ic_man.ic[0];
  if(!ic.improper){
    //build first 3 atoms
    atom_[ic.atomd[0]][0] = 0;
    atom_[ic.atomd[0]][1] = 0;
    atom_[ic.atomd[0]][2] = 0;
    known.insert(ic.atomd[0]);

    atom_[ic.atomd[1]][0] = ic.data[0];
    atom_[ic.atomd[1]][1] = 0;
    atom_[ic.atomd[1]][2] = 0;
    known.insert(ic.atomd[1]);
  
    //
    double d = mdb.bondT(ic.atom_type[0],
			 ic.atom_type[1]).equ_length();
    double a = ic.data[1];
    atom_[ic.atomd[2]][0] = ic.data[0] + d*cos(M_PI-a);
    atom_[ic.atomd[2]][1] = d*sin(M_PI-a);
    atom_[ic.atomd[2]][2] = 0;
    known.insert(ic.atomd[2]);
  }
  else {
    //build first 3 atoms
    atom_[ic.atomd[0]][0] = 0;
    atom_[ic.atomd[0]][1] = 0;
    atom_[ic.atomd[0]][2] = 0;
    known.insert(ic.atomd[0]);

    atom_[ic.atomd[2]][0] = ic.data[0];
    atom_[ic.atomd[2]][1] = 0;
    atom_[ic.atomd[2]][2] = 0;
    known.insert(ic.atomd[2]);
  
    //
    double d = mdb.bondT(ic.atom_type[1],
			 ic.atom_type[2]).equ_length();
    double a = ic.data[1];
    atom_[ic.atomd[1]][0] = ic.data[0] + d*cos(M_PI-a);
    atom_[ic.atomd[1]][1] = d*sin(M_PI-a);
    atom_[ic.atomd[1]][2] = 0;
    known.insert(ic.atomd[1]);
  }
  
  while(known.size()!=atom_.size()){
    bool stop;
    bool left;
    ICManager::IC * ic = ic_man.next(known,left,stop);
    if(stop) break;

    //build
    if(left){
      double coor[3];
      build_left(coor,
		 atom_[ic->atomd[1]].coor(),
		 atom_[ic->atomd[2]].coor(),
		 atom_[ic->atomd[3]].coor(),
		 ic);
      atom_[ic->atomd[0]][0] = coor[0];
      atom_[ic->atomd[0]][1] = coor[1];
      atom_[ic->atomd[0]][2] = coor[2];
      known.insert(ic->atomd[0]);
      ic->used = true;
    }
    else {
      known.insert(ic->atomd[3]);
      double coor[3];
      build_right(coor,
		  atom_[ic->atomd[0]].coor(),
		  atom_[ic->atomd[1]].coor(),
		  atom_[ic->atomd[2]].coor(),
		  ic);
      atom_[ic->atomd[3]][0] = coor[0];
      atom_[ic->atomd[3]][1] = coor[1];
      atom_[ic->atomd[3]][2] = coor[2];
      ic->used = true;
    }
  }
  if(almost::verbose())
    if(known.size()==atom_.size()) aout<<"\t\t>> Generation of protein coor done \n";
  else {
    aout<<"\t\t>> Missing coordinates!!!"<<endl;
    for(int i=0;i<atom_.size();i++)
      if(known.find(i)==known.end()) 
	aout<<"WARNING: UNDEF ATOM "<<atom_name(i)<<endl;
  }
}


void Almost::Solvent::find_solvent_sequence(PDBSegment & segment,
					    MDB & mdb,
					    vector<FragmentData> & data,
					    string fpatch,
					    string lpatch){
  vector<string> sequ = segment.sequence();
  for(int i=0;i<sequ.size();++i){
    string fpatch_ = "NONE";
    string lpatch_ = "NONE";

    if(fpatch=="DEFAULT"){
      fpatch_ = mdb.fragment(sequ[i]).def_first_patch();
    }else{
      fpatch_ = fpatch;
    }
    //check if is SOLVENT PATCH
    if(fpatch_!="NONE"){
      if(mdb.patch_kind(fpatch_)!=TemplatePatch::SOLVENT){
	cerr<<"FATAL ERROR: Patch "<<fpatch_<<" is not an solvent patch!"
	    <<endl;
	exit(-1);
      }
    }
      
    //check if is SOLVENT ...
    int fcode = mdb.fragment_kind(sequ[i]);
    if(fcode!=TemplateFragment::SOLVENT){
      cerr<<"FATAL ERROR: Fragment "<<sequ[i]<<" is not a solvent! "
	  <<" code "<< fcode
	  <<endl;
      exit(-1);
    }

    FragmentData fd(sequ[i],fpatch_,lpatch_);
    data.push_back(fd);
  }
}


void Almost::Solvent::find_solvent_sequence(const vector<string> & sequ,
					    MDB & mdb,
					    vector<FragmentData> & data,
					    string fpatch,
					    string lpatch){
  for(int i=0;i<sequ.size();++i){
    string fpatch_ = "NONE";
    string lpatch_ = "NONE";

    if(fpatch=="DEFAULT"){
      fpatch_ = mdb.fragment(sequ[i]).def_first_patch();
    }else{
      fpatch_ = fpatch;
    }
    //check if is SOLVENT PATCH
    if(fpatch_!="NONE"){
      if(mdb.patch_kind(fpatch_)!=TemplatePatch::SOLVENT){
	cerr<<"FATAL ERROR: Patch "<<fpatch_<<" is not an solvent patch!"
	    <<endl;
	exit(-1);
      }
    }
      
    //check if is SOLVENT ...
    int fcode = mdb.fragment_kind(sequ[i]);
    if(fcode!=TemplateFragment::SOLVENT){
      cerr<<"FATAL ERROR: Fragment "<<sequ[i]<<" is not a solvent! "
	  <<" code "<< fcode
	  <<endl;
      exit(-1);
    }

    FragmentData fd(sequ[i],fpatch_,lpatch_);
    data.push_back(fd);
  }
}

void Almost::Solvent::build_info(){
  if(!almost::verbose()) return;
  char buff[256];
  sprintf(buff,"%lu",fragment_.size());
  logMsg(Almost::SolventPr,string("Solvent generation compleated")
	 +buff+
	 string(" fragments"));

  aout<<"\n\t";
  int l= fragment_.back().size();
  for(int i=0;i<fragment_.size();i++){
    aout<<fragment_[i]<<string(l-fragment_[i].size(),' ')<<" ";
    if(((i+1)%5)==0){
      aout<<"\n\t";
    }
  }
  aout<<"\n\n";
  aout<<"\t>> Build details:"<<"\n\n";
  aout<<"\t>> Number of fragments:   "<<fragment_.size()<<"\n";
  aout<<"\t>> Number of atoms:       "<<atom_.size()<<"\n";
  aout<<"\t>> Number of bonds:       ";
  int nbonds= 0;
  for(int i=0;i<bonds_.size();i++){
    nbonds += bonds_[i].size();
  }
  aout<<nbonds/2<<"\n";
  aout<<"\t>> Number of angles:      "<<angle_.size()<<"\n";
  aout<<"\t>> Number of dihedrals:   "<<dihedral_.size()<<"\n";
  aout<<"\t>> Number of impropers:   "<<improper_.size()<<"\n";
  aout<<"\t>> Number of groups:      "<<group_.size()<<"\n";
  int nfixed=0;
  for(int i=0;i<fixed_excl_.size();i++){
    nfixed += fixed_excl_[i].size();
  }
  aout<<"\t>> Number of fixed excl.: "<<nfixed/2<<"\n";
  aout<<"\n\n";
}

bool Almost::Solvent::is_forward(TDihedral & d){
  for(int i=0;i<4;i++) if(d.atom(i)[0]=='+') return true;
  return false;
}

bool Almost::Solvent::is_backward(TDihedral &d){
  for(int i=0;i<4;i++) if(d.atom(i)[0]=='-') return true;
  return false;
}

bool Almost::Solvent::is_forward(TImproper & d){
  for(int i=0;i<4;i++) if(d.atom(i)[0]=='+') return true;
  return false;
}

bool Almost::Solvent::is_backward(TImproper &d){
  for(int i=0;i<4;i++) if(d.atom(i)[0]=='-') return true;
  return false;
}


void Almost::Solvent::translate(const Almost::Coor3D & coor){
  for(int i=0;i<atom_.size();i++){
    for(int j=0;j<3;j++)
      atom_[i][j] += coor[j];

  }
}

void Almost::Solvent::rotate(double angle, const Coor3D & axis,
		const Coor3D & point){

  int nmove = atom_.size();
  if(nmove==0) return;
    
  vector<vector<double> > M;
  vector<vector<double> > U;
  vector<vector<double> > S;
  vector<vector<double> > I;
  double tmp[3];
 
  //Init matrices ...
  M.resize(3);
  U.resize(3);
  S.resize(3);
  I.resize(3);
  for(int i=0;i<3;i++){
    M[i].resize(3);
    U[i].resize(3);
    S[i].resize(3);
    I[i].resize(3);
  }
      
  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++){
      if(i!=j) I[i][j]=0;
      else I[i][j]=1;
    }

  Coor3D n = axis;
  double norm= 0;
  for(int i=0;i<3;i++) norm +=n[i]*n[i];

  norm = sqrt(norm);
  for(int i=0;i<3;i++) n[i] /=norm;

  for(int i =0; i<3;++i){
    for(int j = 0;j<3;++j){
      U[i][j]=n[i]*n[j];
      S[i][j]=0; 
    }
  }

  S[0][1] = -n[2];
  S[1][0] = -S[0][1];
  
  S[0][2] = n[1];
  S[2][0] =-S[0][2];
  
  S[1][2] = -n[0];
  S[2][1] = -S[1][2];
  
  double cosv = cos(angle);
  double sinv = sin(angle);
  
  for(int i =0; i<3;++i){
    for(int j = 0;j<3;++j){
      M[i][j] = U[i][j]+ cosv*(I[i][j]-U[i][j])+sinv*S[i][j];
    }
  }

  for(int i=0;i<nmove;++i){
    //translate
    for(int j=0;j<3;j++){
      atom_[i][j] -=point[j];
    }
    //rota
    for(int j=0;j<3;++j){
      tmp[j] = 0;
      for(int k=0;k<3;++k){
	tmp[j] += M[j][k]*atom_[i][k];
      }
    }

    for(int j=0;j<3;++j){      
      atom_[i][j] =tmp[j];
    }


    //TRANSLATE BACK
    for(int j=0;j<3;++j){
      atom_[i][j] += point[j];
    }
  }
}

Almost::Coor3D Almost::Solvent::cog() const {
  Coor3D c;
  if(atom_.size()==0) return c;

  for(int i=0;i<atom_.size();i++)
    for(int j=0;j<3;j++)
      c[j] += atom_[i][j];


  for(int j=0;j<3;j++)
    c[j] /= atom_.size();

  return c;
}

Almost::Coor3D Almost::Solvent::com(const MDB & mdb) const {
  Coor3D c;
  if(atom_.size()==0) return c;
  double M = 0;
  for(int i=0;i<atom_.size();i++){
    double m = mdb.atomT(atom_[i].type_name()).mass(); 
    M += m;
    for(int j=0;j<3;j++){
      c[j] += m*atom_[i][j];
    }
  }
  for(int j=0;j<3;j++)
    c[j] /= M;

  return c;
}

double Almost::Solvent::distanceP2(int a1, int a2) const {
  double d2 = 0;
  for(int i=0;i<3;i++)
    d2 += (atom_[a1][i]-atom_[a2][i])*(atom_[a1][i]-atom_[a2][i]);
  return d2;
}

double Almost::Solvent::distance(int a1, int a2) const {
  double d2 = 0;
  for(int i=0;i<3;i++)
    d2 += (atom_[a1][i]-atom_[a2][i])*(atom_[a1][i]-atom_[a2][i]);
  return sqrt(d2);
}

void Almost::Solvent::setup_parameters(const MDB & mdb){
  //Atoms
  {
    int b = atom_offset();
    int e = b + atom_size();
    Solvent & p = * this;
    for(int i=b;i<e;i++){
      try {
	SasaT sasat = mdb.sasaT(p[i].type_name());
	p[i].set_sasa_solv(sasat.solv());
	p[i].set_sasa_vdw(sasat.vdw());
	p[i].set_sasa_atom_par(sasat.atom_par());
	p[i].set_has_sasa(true);
      }
      catch(...){
	p[i].set_has_sasa(false);
// 	aout<<"\t>> WARNING: no SASA parameters for atom: "<<p[i].name()
// 	    <<" "<<p[i].type_name()<<endl;
      }

      try {
	EEF1T eef1t = mdb.eef1T(p[i].type_name());
	p[i].set_has_eef1(true);
	p[i].set_eef1_vol(eef1t.Vol());
	p[i].set_eef1_gref(eef1t.Gref());
	p[i].set_eef1_gfree(eef1t.Gfree());
	p[i].set_eef1_href(eef1t.Href());
	p[i].set_eef1_cpref(eef1t.CPref());
	p[i].set_eef1_sigw(eef1t.Sigw());
      }
      catch(...){
	p[i].set_has_eef1(false);
// 	aout<<"\t>> WARNING: no EEF1 parameters for atom: "<<p[i].name()
// 	    <<" "<<p[i].type_name()<<endl;
      }
    }
    
  }
  {
    //Bonds pars
    int b = atom_offset();
    int e = b + atom_size();
    Solvent & p = * this;
    for(int i=b;i<e;i++){
      const vector<int> & tmp_bonds = bonds(i);
      string name1 = p[i].type_name();
      for(int j=0;j<tmp_bonds.size();j++){
	//      if(bonds[j]<i) continue;
	string name2 = p[tmp_bonds[j]].type_name();
	BondT bt = mdb.bondT(name1,name2);
	pair<int,int> aa(i,tmp_bonds[j]);
	bond_data_[aa].set_bond_const(bt.bond_const());
	bond_data_[aa].set_equ_length(bt.equ_length());
      }
    }
  }
  //Angles
  {
    int size = angles().size();
    Solvent & p = * this;
    for(int i=0;i<size;i++){
      string name1 = p[angles()[i][0]].type_name();
      string name2 = p[angles()[i][1]].type_name();
      string name3 = p[angles()[i][2]].type_name();
      try {
	vector<AngleT> at = mdb.angleT(name1,name2,name3);
	for(int aa=0;aa<at.size();aa++){
	  Kernel::AngleData ad;
	  ad.set_angle_const(at[aa].angle_const());
	  ad.set_equ_angle(at[aa].equ_angle());

	  ad.set_ub_angle_const(at[aa].ub_angle_const());
	  ad.set_ub_equ_angle(at[aa].ub_equ_angle());

	  ad.set_has_ub(at[aa].has_ub());
	  angle_.add_angle_data(i,ad);
	  
	}
      }
      catch(...){
	aout<<"WARNING: angle data not in MDB for "<<name1+"-"+name2+"-"+name3<<" ignored"<<endl;
      }
    }
  }
  //Dihedral
  {
    int size = dihedrals().size();
    Solvent & p = * this;
    for(int i=0;i<size;i++){
      string name1 = p[dihedrals()[i][0]].type_name();
      string name2 = p[dihedrals()[i][1]].type_name();
      string name3 = p[dihedrals()[i][2]].type_name();
      string name4 = p[dihedrals()[i][3]].type_name();
      try {
	vector<DihedralT> at = mdb.dihedralT(name1,name2,name3,name4);
	for(int aa=0;aa<at.size();aa++){
	  Kernel::DihedralData dd;
	  dd.set_dihe_const(at[aa].dihe_const());
	  dd.set_n(at[aa].n());
	  dd.set_delta(at[aa].delta());	  
	  dihedral_.add_dihedral_data(i,dd);	  
	}
      }
      catch(...){
	aout<<"WARNING: dihedral data not in MDB for "<<name1+"-"+name2+"-"+name3<<" ignored"<<endl;
      }
    }
  }  

  //Improper
  {
    int size = impropers().size();
    Solvent & p = * this;
    for(int i=0;i<size;i++){
      string name1 = p[impropers()[i][0]].type_name();
      string name2 = p[impropers()[i][1]].type_name();
      string name3 = p[impropers()[i][2]].type_name();
      string name4 = p[impropers()[i][3]].type_name();
      try {
	vector<ImproperT> at = mdb.improperT(name1,name2,name3,name4);
	for(int aa=0;aa<at.size();aa++){
	  Kernel::ImproperData id;
	  id.set_imph_const(at[aa].imph_const());
	  id.set_n(at[aa].n());
	  id.set_equ_angle(at[aa].equ_angle());	  
	  improper_.add_improper_data(i,id);	  
	}
      }
      catch(...){
	aout<<"WARNING: improper data not in MDB for "<<name1+"-"+name2+"-"+name3<<" ignored"<<endl;
      }
    }
  }  


  //NB
  {
    set<string> reg_types;
    Solvent &p = *this;
    int b = atom_offset();
    int e = b + atom_size();
    for(int i=b;i<e;i++){
      string type = p[i].type_name();
      if(reg_types.find(type)==reg_types.end()){
	try {
	  NBondT bnt = mdb.nbondT(p[i].type_name());
	  nbdata_.set_emin(type,bnt.emin());
	  nbdata_.set_emin14(type,bnt.emin14());
	  nbdata_.set_sigma(type,bnt.sigma());
	  nbdata_.set_sigma14(type,bnt.sigma14());
	  reg_types.insert(type);
	}
	catch(...){
	  aout<<"\t>>WARNING atom type "<<type<<" not defined in MDB"<<endl;
	}
      }
    }

    set<string>::iterator iter1 = reg_types.begin(), iter2 = reg_types.begin(),
      end = reg_types.end();
    while(iter1!=end){
      iter2 = reg_types.begin();
      while(iter2!=end){
	if(mdb.is_nbfixed(*iter1,*iter2)){
	  NBFix nbfix = mdb.nbfix(*iter1,*iter2);
	  nbdata_.set_emin_fix(*iter1,*iter2,nbfix.e_min());
	  nbdata_.set_emin14_fix(*iter1,*iter2,nbfix.e_min14());
	  nbdata_.set_sigma_fix(*iter1,*iter2,nbfix.r_vdw());
	  nbdata_.set_sigma14_fix(*iter1,*iter2,nbfix.r_vdw14());	  
	}
	++iter2;
      }
      ++iter1;
    }
  }
}