/***************************************************************************
                          molecule.cpp  -  description
                             -------------------
    begin                : Tue Nov 18 01:45:35 2003
    copyright            : (C) 1999-2008 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 <almostexception.h>
#include <almosttypes.h>
#include <almosttemplate.h>
#include <pdb.h>
#include <mdb/mdb.h>
#include <moveset/bondrotation.h>
#include <algorithm>

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


void Almost::Molecule::build(PDBSegment & segment,MDB & mdb,
			    string fpatch, string lpatch) 
{
  try{
    //Step 1: Find sequence from PDB
    vector<FragmentData> data;
    {
      find_molecule_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;
	throw "FATAL ERROR";
      }
    } 
    //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;
	throw "FATAL ERROR";
      }
    }
    //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;
	throw "FATAL ERROR";
      }
    }
    //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 11
    //Add groups
    {
      add_groups(mdb,data);
    }

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



void Almost::Molecule::build_fast(PDBSegment & segment,MDB & mdb,
				 string fpatch, string lpatch) 
{
  try{
    //Step 1: Find sequence from PDB
    vector<FragmentData> data;
    {
      find_molecule_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;
	throw "FATAL ERROR";
      }
    } 
    //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;
	throw "FATAL ERROR";
      }
    }
    //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;
	throw "FATAL ERROR";
      }
    }
    //Step 7: Done

    //Step 8
    //Add angles
    {
      //      add_angles(mdb,data);
      logMsg(Almost::MoleculePr,"No angles build");
    }
    //Step 8: Done

    //Step 9
    //Add dihedrals
    {
      //      add_dihedrals(mdb,data);
      logMsg(Almost::MoleculePr,"No dihedrals build");
    }
    //Step 9: Done

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

    //Step 11
    //Add groups
    {
      add_groups(mdb,data);
    }
    //Step 12
    //Print info
    {
      build_info();
    }
  }
  catch(AlmostException e){
    throw;
    //AEX(e.error());
    //    throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,e.error());
  }
}


void Almost::Molecule::build_missing(PDBSegment & segment,MDB & mdb,
				    string fpatch, string lpatch){
  try{
    vector<FragmentData> data;
    vector<string> sequ = segment.sequence(PDBSegment::ICODE);
    {
      find_molecule_sequence(sequ,mdb,data,fpatch,lpatch);
    }

    {
      add_segment_fragments(data);
    }

    {
      add_segment_atoms(sequ,mdb,data);
      define_known_atoms(segment,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;
	throw "FATAL ERROR";
      }
    }

    {
      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;
  }
}



int Almost::Molecule::build_sequ(const vector<string> &sequ, MDB & mdb,
				string fpatch, string lpatch){
  int error = 0;
  try{
    vector<FragmentData> data;
    {
      find_molecule_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;
	throw "FATAL ERROR";
      }
    }

    {
      add_angles(mdb,data);
    }

    {
      add_dihedrals(mdb,data);
    }

    {
      add_impropers(mdb,data);
    }

    {
      add_groups(mdb,data);
    }

    {
      try {
	generate_coor(mdb,data);	
      }
      catch(...){
	error = 1;
      }
    }
    {
      setup_parameters(mdb);
    }
    {
      build_info();
    }

  }
  catch(AlmostException e){
    throw;
  }
  return error;
}

int Almost::Molecule::build_template(const vector<string> &sequ, MDB & mdb,
				    string fpatch, string lpatch){
  int error = 0;
  try{
    vector<FragmentData> data;
    {
      find_molecule_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;
	throw "FATAL ERROR";
      }
    }

    {
      add_angles(mdb,data);
    }

    {
      add_dihedrals(mdb,data);
    }

    {
      add_impropers(mdb,data);
    }

    {
      add_groups(mdb,data);
    }

    //     {
    //       try {
    // 	generate_coor(mdb,data);	
    //       }
    //       catch(...){
    // 	error = 1;
    //       }
    //     }
    {
      setup_parameters(mdb);
    }
    {
      build_info();
    }

  }
  catch(AlmostException e){
    throw;
  }
  return error;
}



void Almost::Molecule::hbuild(MDB & mdb,
			     string fpatch,
			     string lpatch){

  vector<FragmentData> data;
  //data.resize(fragment_size());
  for(int i = 0;i<fragment_size();i++){
    FragmentData d(fragment_[i],"NONE","NONE");
    if(i==0){ 
      if(fpatch=="DEFAULT"){
	d.fpatch = mdb.fragment(fragment_[i]).def_first_patch();
      }
      else d.fpatch = fpatch;
    }
    if(i==fragment_size()-1){
      if(lpatch=="DEFAULT"){
	d.lpatch = mdb.fragment(fragment_[i]).def_last_patch();
      }
      else d.lpatch = lpatch;
    }
    data.push_back(d);
  }

  for(int i=0;i<atom_.size();i++){
    if(atom_[i].name()[0]=='H')
      atom_[i][0] = atom_[i][1] = atom_[i][2] = 9999;
  }
  generate_missing(mdb,fpatch,lpatch);
}

void Almost::Molecule::generate_missing(MDB & mdb,
				       string fpatch,
				       string lpatch){

  vector<FragmentData> data;
  //data.resize(fragment_size());
  for(int i = 0;i<fragment_size();i++){
    FragmentData d(fragment_[i],"NONE","NONE");
    if(i==0){ 
      if(fpatch=="DEFAULT"){
	d.fpatch = mdb.fragment(fragment_[i]).def_first_patch();
      }
      else d.fpatch = fpatch;
    }
    if(i==fragment_size()-1){
      if(lpatch=="DEFAULT"){
	d.lpatch = mdb.fragment(fragment_[i]).def_last_patch();
      }
      else d.lpatch = lpatch;
    }
    data.push_back(d);
  }

  generate_coor(mdb,data);
  //  aout<<"DEBUG Almost::Molecule::generate_missing done "<<endl;
}


void Almost::Molecule::patch(string frag, string p, MDB & mdb){
  aout<<"Not implemented sorry ....."<<endl;
  //   int fd = find_fragment(frag);
  //   if(fd<0){
  //     aout<<"Couldn't find fragment "<<frag<<endl;
  //     return;
  //   }

  //   Template tpatch = mdb.patch(p);


  //   //delete

  //   //atoms
  //   {
  //     vector<string> atom = tpatch.del_atoms();
  //     for(int i=0;i<atom.size();i++){
  //       string a = frag+"/"+atom[i];
  //       int na = find_atom(a);
  //       if(na<0){
  // 	aout<<"Couldn't find atom "+a<<endl;
  // 	return;
  //       }

  //       remove_atom(na);
  //     }
  //   }

  //   //bonds
  //   {
  //     vector<TBond> bond = tpatch.del_bond();
  //     for(int i=0;i<bond.size();i++){
  //       string a1 = frag+"/"+bond[i].first();
  //       string a2 = frag+"/"+bond[i].second();
  //       int na1 = find_atom(a1);
  //       int na2 = find_atom(a2);
  //       if(na1<0){
  // 	aout<<"Couldn't find bond "+a1<<endl;
  // 	return;
  //       }

  //       if(na2<0){
  // 	aout<<"Couldn't find bond "+a2<<endl;
  // 	return;
  //       }

  //       if(bonded(na1,na2)){
  // 	remove_bond(na1,na2);
  //       }
  //       else {
  // 	aout<<"Atoms : "<<a1<<" "<<a2<<" not bonded"<<endl;
  // 	return;
  //       }
  //     }
  //   }

  //angles
  //   {
  //     vector<TAngle> angle = tpatch.del_angle();
  //     for(int i=0;i<angle.size();i++){
  //       string a1 = frag+"/"+angle[i].first();
  //       string a2 = frag+"/"+angle[i].second();
  //       string a2 = frag+"/"+angle[i].second();
  //       int na1 = find_atom(a1);
  //       int na2 = find_atom(a2);
  //       if(na1<0){
  // 	aout<<"Couldn't find angle "+a1<<endl;
  // 	return;
  //       }

  //       if(na2<0){
  // 	aout<<"Couldn't find angle "+a2<<endl;
  // 	return;
  //       }

  //       if(angleed(na1,na2)){
  // 	remove_angle(na1,na2);
  //       }
  //       else {
  // 	aout<<"Atoms : "<<a1<<" "<<a2<<" not angleed"<<endl;
  // 	return;
  //       }
  //     }
  //   }

  //   //dihedral
  //   {
  //     vector<TDihedral> dihedral = tpatch.del_dihedral();
  //     for(int i=0;i<dihedral.size();i++){
  //       string a[4];
  //       vector<int> na(4);
      
  //       for(int j=0;j<4;j++){
  // 	a[j] = frag+"/"+dihedral[i].atom(j);
  // 	na[j] = find_atom(a[j]);
  // 	if(na[j]<0){
  // 	  aout<<"Couldn't find dihedral "+a[j]<<endl;
  // 	  return;
  // 	}
  //       }
      
  //       vector<vector<int> >::iterator iter = dihedrals_.find(na);
  //       if(iter!=dihedrals_.end())
  // 	dihedrals_.remove(iter);
  //       else {
  // 	aout<<"Atoms: ";
  // 	for(int j=0;j<4;j++) aout<<a[j]<<" ";
  // 	aout<<"do not form a dihedral"<<endl;
  // 	return;
  //       }	
  //     }
  //   }

  //   //improper
  //   {
  //     vector<TImproper> improper = tpatch.del_improper();
  //     for(int i=0;i<improper.size();i++){
  //       string a[4];
  //       vector<int> na(4);
      
  //       for(int j=0;j<4;j++){
  // 	a[j] = frag+"/"+improper[i].atom(j);
  // 	na[j] = find_atom(a[j]);
  // 	if(na[j]<0){
  // 	  aout<<"Couldn't find improper "+a[j]<<endl;
  // 	  return;
  // 	}
  //       }
      
  //       vector<vector<int> >::iterator iter = impropers_.find(na);
  //       if(iter!=impropers_.end())
  // 	impropers_.remove(iter);
  //       else {
  // 	aout<<"Atoms: ";
  // 	for(int j=0;j<4;j++) aout<<a[j]<<" ";
  // 	aout<<"do not form an improper"<<endl;
  // 	return;
  //       }	
  //     }
  //   }
  
  //   //insert-modify
  //   //atoms
  

}



void Almost::Molecule::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;
      throw "FATAL ERROR";
    }
  }
}



void Almost::Molecule::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;
    if(segment[i].iCode()!=' ')
      resName += segment[i].iCode();
    int fragd = find_fragment_by_name(resName);
    if(fragd<0){
      cerr<<"No such fragment "<<resName<<endl;
      throw "FATAL ERROR";
    }
    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());
      atom.set_canonical_name(tatom.canonical_name());
      int atomd = add_atom(atom,fragd);
      if(atomd<0){
	aout<<"FATAL ERROR: couldn't add atom "<<resName+"/"+segment[i][j].name
	    <<" to segment"<<endl;
	throw "FATAL ERROR";
      }
    }
  }
}


void Almost::Molecule::add_segment_atoms( const vector<string> & sequ,
					 MDB & mdb,
					 vector<FragmentData> & data){
  for(int i=0;i<data.size();i++){
    string resName = data[i].fragment;
    //    sprintf(buff,"_%i",resSeq);
    //    resName +=buff;
    int fragd = find_fragment_by_name(resName);
    if(fragd<0){
      cerr<<"No such fragment "<<resName<<endl;
      throw "FATAL ERROR";
    }

    vector<string> atoms;
    string _fpatch;
    string _lpatch;
    _fpatch="NONE";
    _lpatch="NONE";
    if(i==0){ 
      if(data[i].fpatch=="DEFAULT"){
	_fpatch = mdb.fragment(resName).def_first_patch();
      }else{
	_fpatch=data[i].fpatch;
      }
    } 
    if(i==(data.size()-1)){
      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());
      atom.set_canonical_name(tatom.canonical_name());
      int atomd = add_atom(atom,fragd);
      if(atomd<0){
	aout<<"FATAL ERROR: couldn't add atom "<<resName+"/"+atoms[j]
	    <<" to segment"<<endl;
	throw "FATAL ERROR";
      }

    }
    //    for(int j=0;j<sequ.size();++j){
    //      aout<<"DEB building segment residue"<<sequ[i]<<endl;
    //       TAtom tatom = mdb.atom(resName+"/"+segment[i][j].name,
    // 			     data[i].fpatch,data[i].lpatch);
    //       Atom atom(segment[i][j].name,tatom.type_name(),
    // 		segment[i][j].x,segment[i][j].y,segment[i][j].z,
    // 		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;
    // 	throw "FATAL ERROR";
    //       }
    //    }
  }
}

void Almost::Molecule::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;
      throw "FATAL ERROR";
    }
    //    string fname = fragment_name(fd,Molecule::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_by_canonical_name(data[fd].fragment+string("/")+(*iter));
      if(indx<0){
	aout<<"FATAL ERROR: couldn't find atom "
	    <<data[fd].fragment+string("/")+(*iter)<<endl;
	throw "FATAL ERROR";
      }
      add_excl(i,indx);
      ++iter;
    }
  }
}

void Almost::Molecule::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::Molecule::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_by_canonical_name(fragment_[i]+"/"+h_atoms[k]);
	if(atomd<=-1){
	  // 	      aout <<atomd<<endl;
	  aout<<"WARNING: "<<"Hydrogen "+fragment_[i]+"/"
	    +h_atoms[k]+
	    " missing for atom "
	    + fragment_[i]+"/"+heavy_atoms[j]+
	    " in fragment "+fragment_[i]
	      <<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::Molecule::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){
	aout<< "WARNING 1: "<< "Atom " <<fragment_[i]<<"/"<< atoms[j]
	    << " not found in fragment "
	    << fragment_[i] <<endl;
	// 	throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,
	// 			      "Atom "+fragment_[i]+"/"+ atoms[j]+
	// 			      " not found in fragment "+
	// 			      fragment_[i]);
	continue;
      }
      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_by_canonical_name(search_name);
	if(atomd2<0){
	  aout<< "WARNING 2: "<<"Atom "<<search_name
	      << " not found in fragment "
	      << fragment_[i]<<endl;
	  // 	  throw AlmostException(__FILE__,__LINE__,__PRETTY_FUNCTION__,
	  // 				"Atom "+search_name+
	  // 				" not found in fragment "+
	  // 				fragment_[i]);
	  continue;
	}
	if(!add_bond(atomd1,atomd2,bond_kind[k])){
	  aout<<"Warning: bond beteewn "<<atomd1+1<<" "<<atomd2+1<<" not added!"<<endl;
	}
      }
    }
  } 
}



int Almost::Molecule::find_atom(string name) const{
  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);
  //   cout<<"DEBUG nfrag "<<nfrag<<endl;
  if(nfrag==-1) return -1;
  for(int i=0;i<fragment_atoms_[nfrag].size();++i){
    if(atom_[fragment_atoms_[nfrag][i]-atom_offset()].name()==atom) 
      return fragment_atoms_[nfrag][i]-atom_offset();
    //     else cout<<"DEBUG "<<i<<" "<<" "<<fragment_atoms_[nfrag][i]<<" "
    // 	     <<atom_[fragment_atoms_[nfrag][i]].name()<<endl;
  }


  return -1;
}

int Almost::Molecule::find_atom_by_canonical_name(string name) const{
  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);
  //   cout<<"DEBUG nfrag "<<nfrag<<endl;
  if(nfrag==-1) return -1;

  for(int i=0;i<fragment_atoms_[nfrag].size();++i){
    if(atom_[fragment_atoms_[nfrag][i]-atom_offset()].canonical_name()==atom) 
      return fragment_atoms_[nfrag][i]-atom_offset();
    //     else cout<<"DEBUG "<<i<<" "<<" "<<fragment_atoms_[nfrag][i]<<" "
    // 	     <<atom_[fragment_atoms_[nfrag][i]].name()<<endl;
  }

  return -1;
}

void Almost::Molecule::add_angles(MDB & mdb,vector<FragmentData> & data){
  angle_generate();
  //   aout<<"DEBUG ANGLES ...\n";
  //   for(int i=0;i<angle_.size();i++){
  //     for(int j=0;j<angle_[i].size();j++){
  //       int a = angle_[i][j];
  //       aout<< a+1<<" "<<atom_[a].name()<<" "<<atom_[a].type_name()<<" ";
  //     }
  //     aout<<"\n";
  //   }
}
    
void Almost::Molecule::add_dihedrals(MDB & mdb,vector<FragmentData> & data){

  //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];
      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_by_canonical_name(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;
	  throw "FATAL ERROR";
	}
	dihe_.push_back(index[k]);
      }

      if(!dihedral_.add_dihedral(dihe_)) 
	aout<<"Warning: cant add dihedral"<<endl;
    }
  }
  if(dihedral_.size()==0){  
    dihedral_generate();
  }
}

void Almost::Molecule::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
    //aout<<fragment_[i]<<" number if imph "<<imph.size()<<endl;
    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];
      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
	      <<" "<<name[k]
	      <<endl;
	  throw "FATAL ERROR";
	}
	imph_.push_back(index[k]);
      }
      if(!improper_.add_improper(imph_))
	aout<<"Warning: can't add improper"<<endl;
    }
  }
}


void Almost::Molecule::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;
	name  = fragment_[i]+"/"+name;
	index = find_atom_by_canonical_name(name);
	tmp.push_back(index);
      }
      sort(tmp.begin(),tmp.end());
      group_.push_back(tmp);
    }
  }
}

void Almost::Molecule::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_by_canonical_name(name[k]);
	if(atomd[k]<0){
	  aout<<"FATAL ERROR: "<<name[k]<<
	    " not found\n";
	  for(int q=0;q<atom_.size();q++){
	    cout<<atom_[q].name()<<" ";
	  }
	  cout<<endl;
	  aout<<"Fragment "<<i<<" "<< data[i].fpatch
	      <<" "<< data[i].lpatch<<" "<<ic.size()<<"\n";
	  throw "FATAL ERROR";
	}
      }

      //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;
  {
    for(int i=0;i<atom_.size();i++){
      if((atom_[i][0]!=9999)&&
	 (atom_[i][1]!=9999)&&
	 (atom_[i][2]!=9999)) known.insert(i);
    }
    //    aout<<"DEBUG known "<<known.size()<<endl;
  }
  if((ic_man.ic.size()==0)&&(known.size()!=atom_.size())){
    aout<<"\t\t WARNING: Undefined atoms\n";
    return;
  }

  if(known.size()<3){
    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;
    }
  }

  //2nd try ...
//   if(known.size()!=atom_.size()){
//     aout<<"\t>> Generating synthetic internal coordinates.\n";
//     set<int> unkn;
//     for(int i=0;i<atom_.size();i++){
//       if(known.find(i)==known.end()){
// 	unkn.insert(i);
//       }
//     }

//     //Search dihe with i as first atom
//     int count = 0;
//     const vector<vector<int> > & dd = dihedral_.dihedrals();
//     set<int>::iterator iter = unkn.begin(),
//       end = unkn.end();
//     while(iter!=end){
//       int aa = *iter;
//       //find a dihe for aa
//       int stop = 0;
//       aout<<">> Searching for "<<aa<<"\n";
//       for(int qq=0;qq<dd.size();qq++){
// 	if(aa==dd[qq][0]||aa==dd[qq][3]){
// 	  //found 
// 	  ic_man.ic.push_back(ICManager::IC());
// 	  for(int zz=0;zz<4;zz++){
// 	    ic_man.ic.back().atom[zz] = atom_[dd[qq][zz]].name();
// 	    ic_man.ic.back().atom_type[zz] = atom_[dd[qq][zz]].type_name();
// 	    ic_man.ic.back().atomd[zz] = dd[qq][zz];
// 	  }
// 	  for(int zz=0;zz<5;zz++){
// 	    ic_man.ic.back().data[zz] = 0;
// 	  }
// 	  ic_man.ic.back().improper = false;
// 	  aout<<"\t>>  synthetic internal coordinates generated for:\n\t   ";
// 	  for(int zz=0;zz<4;zz++){
// 	    aout<<atom_name(dd[qq][zz])<<" ";
// 	  } 
// 	  aout<<"\n";
// 	  ++count;
// 	  stop = 1;
// 	}
// 	else {
// 	  for(int zz=0;zz<4;zz++){
// 	    cout<<dd[qq][zz]<<" ";
// 	  }
// 	  cout<<endl;
// 	}
// 	if(stop)
// 	  break;
//       }
//       ++iter;
//     }
//     if(count){      
//       ic_man.sort();
//       ic_man.fill(mdb);
//       aout<<"\t>> Adding other atoms...\n";
//       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;
// 	}
//       }      
//     }
//   }  
  //end of 2nd try
  if(almost::verbose())
    if(known.size()==atom_.size()) aout<<"\t>> Generation of molecule 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)<<" "<<atom_name(i,Molecule::SHORT)<<endl;
	  if(atom_name(i,Molecule::SHORT)=="H"||atom_name(i,Molecule::SHORT)=="HN"){
	    //Try to recover   ....
	    int fd = find_fragment(i);
	    int ca = find_atom_by_canonical_name(fragment_name(fd)+"/CA");
	    int cb = find_atom_by_canonical_name(fragment_name(fd)+"/CB");
	    int N = find_atom_by_canonical_name(fragment_name(fd)+"/N");
	    if(fd==-1||ca==-1||cb==-1||N==-1)
	      throw "UNDEFINED ATOMS";

	    cout<<"DEBUG "<<fd<<" "<<ca<<" "<<cb<<endl;
	    double n[3];
	    n[0] = atom_[cb][0] - atom_[ca][0];
	    n[1] = atom_[cb][1] - atom_[ca][1];
	    n[2] = atom_[cb][2] - atom_[ca][2];
	    double nn = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
	    nn = 1.0/sqrt(nn);
	    atom_[i][0] = atom_[N][0]+n[0]*nn;
	    atom_[i][1] = atom_[N][1]+n[1]*nn;
	    atom_[i][2] = atom_[N][2]+n[2]*nn;
	  }
	  else 
	    throw "UNDEFINED ATOMS";
	}
    }
}


void Almost::Molecule::define_known_atoms(PDBSegment & segment,
					 MDB & mdb,
					 vector<FragmentData> & data){

  char buff[16];
  int def = 0;
  for(int i=0;i<data.size();i++){
    int resSeq = segment[i].fragSeq();
    string resName = segment[i].name();
    sprintf(buff,"_%i",resSeq);
    resName +=buff;
    if(segment[i].iCode()!=' ')
      resName += segment[i].iCode();
    int fragd = find_fragment_by_name(resName);
    if(fragd<0){
      cerr<<"No such fragment "<<resName<<endl;
      throw "FATAL ERROR";
    }
    for(int j=0;j<segment[i].size();++j){
      int atomd = find_atom_by_canonical_name(resName+"/"+segment[i][j].name);
      if(atomd<0){
	//	aout<<"FATAL ERROR: couldn't find atom "<<resName+"/"+segment[i][j].name
	//	    <<" in segment template"<<endl;
	//	throw "FATAL ERROR";
      }
      else {
	atom_[atomd][0] = segment[i][j].x;
	atom_[atomd][1] = segment[i][j].y;
	atom_[atomd][2] = segment[i][j].z;
	++def;
      }
    }
  }
  if(almost::verbose())
    aout<<"\t>> Segment contains "<<def<<" "<<atom_size()-def<<" ("<<atom_size()<<") defined and missing atoms\n";
}

void Almost::Molecule::find_molecule_sequence(PDBSegment & segment,
					    MDB & mdb,
					    vector<FragmentData> & data,
					    string fpatch,
					    string lpatch){
  vector<string> sequ = segment.sequence(PDBSegment::ICODE);  
  for(int i=0;i<sequ.size();++i){
    string fpatch_ = "NONE";
    string lpatch_ = "NONE";
    if(i==0){
      if(fpatch=="DEFAULT"){
	fpatch_ = mdb.fragment(sequ[i]).def_first_patch();
      }else{
	fpatch_ = fpatch;
      }
      //check if is AMINOACID PATCH
      if(fpatch_!="NONE"){
	if(mdb.patch_kind(fpatch_)!=TemplatePatch::AMINOACID){
	  cerr<<"FATAL ERROR: Patch "<<fpatch_<<" is not an aminoacid patch!"
	      <<endl;
	  throw "FATAL ERROR";
	}
      }
    }
    if(i==sequ.size()-1){
      if(lpatch=="DEFAULT"){
	lpatch_ = mdb.fragment(sequ[i]).def_last_patch();
      }else{
	lpatch_ = lpatch;
      }
      //check if is AMINOACID PATCH
      if(lpatch_!="NONE"){
	if(mdb.patch_kind(lpatch_)!=TemplatePatch::AMINOACID){
	  cerr<<"FATAL ERROR: Patch "<<lpatch_<<" is not an aminoacid patch!"
	      <<endl;
	  throw "FATAL ERROR";
	}
      }
    }
      
    //check if is AMINOACID ...
    int fcode = mdb.fragment_kind(sequ[i]);
    if(fcode!=TemplateFragment::AMINOACID){
      cerr<<"FATAL ERROR: Fragment "<<sequ[i]<<" is not an aminoacid! "
	  <<" code "<< fcode
	  <<endl;
      throw "FATAL ERROR";
    }

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


void Almost::Molecule::find_molecule_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(i==0){
      if(fpatch=="DEFAULT"){
	fpatch_ = mdb.fragment(sequ[i]).def_first_patch();
      }else{
	fpatch_ = fpatch;
      }
      //check if is AMINOACID PATCH
      if(fpatch_!="NONE"){
	if(mdb.patch_kind(fpatch_)!=TemplatePatch::AMINOACID){
	  cerr<<"FATAL ERROR: Patch "<<fpatch_<<" is not an aminoacid patch!"
	      <<endl;
	  throw "FATAL ERROR";
	}
      }
    }
    if(i==sequ.size()-1){
      if(lpatch=="DEFAULT"){
	lpatch_ = mdb.fragment(sequ[i]).def_last_patch();
      }else{
	lpatch_ = lpatch;
      }
      //check if is AMINOACID PATCH
      if(lpatch_!="NONE"){
	if(mdb.patch_kind(lpatch_)!=TemplatePatch::AMINOACID){
	  cerr<<"FATAL ERROR: Patch "<<lpatch_<<" is not an aminoacid patch!"
	      <<endl;
	  throw "FATAL ERROR";
	}
      }
    }
      
    //check if is AMINOACID ...
    int fcode = mdb.fragment_kind(sequ[i]);
    if(fcode!=TemplateFragment::AMINOACID){
      cerr<<"FATAL ERROR: Fragment "<<sequ[i]<<" is not an aminoacid! "
	  <<" code "<< fcode
	  <<endl;
      throw "FATAL ERROR";
    }

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

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

  aout<<"\n\t";
  int l=0;
  for(int i=0;i<fragment_.size();i++){
    if(l<fragment_[i].size()) l = fragment_[i].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::Molecule::is_forward(TDihedral & d){
  for(int i=0;i<4;i++) if(d.atom(i)[0]=='+') return true;
  return false;
}

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

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

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

 

void Almost::Molecule::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::Molecule::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::Molecule::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::Molecule::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::Molecule::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::Molecule::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);
}


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

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




void Almost::Molecule::setup_parameters(const MDB & mdb){
  //Atoms
  {
    int b = atom_offset();
    int e = b + atom_size();
    Molecule & 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();
    Molecule & 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();
    Molecule & 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();
    Molecule & 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+name4<<" ignored"<<endl;
      }
    }
  }  

  //Improper
  {
    int size = impropers().size();
    Molecule & 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+" "+name4<<" ignored"<<endl;
      }
    }
  }  


  //NB
  {
    set<string> reg_types;
    Molecule &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;
    }
  }
}


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