/***************************************************************************
                          molalgo.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.                                   *
 *                                                                         *
 ***************************************************************************/
#ifndef __MOLALGO__
#define __MOLALGO__

#include <almost.h>
#include <molecules/molecules.h>
#include <io/formostream.h>
#include <almosttypes.h>
#include <mdb.h>
#include <select.h>
#include <algorithm/molalgo.h>
#include <kernel/rmsd.h>
#include <kernel/rgyr.h>
#include <moveset.h>

extern "C" {
  int  nsc_dclm(
		double * ,	/* atom coordinates xyz0, xyz1, ... */
		double * ,	/* atom radii r0, r1, r2, ... */
		int ,         /* number of atoms */
		int ,		/* number of dots per fully accessible sphere */
		int ,		/* flag : dots, volume and/or area per atom */
		double * ,	/* 1 output: overall area */
		double ** ,   /* 2 output: pointer to list of areas per atom */
		double * ,	/* 3 output: overall volume */
		double ** ,	/* 4 output: pointer to list of surface dots x0, y0, z0, ... */
		int *		/* 5 output: number of surface dots */
		);
}

#define FLAG_DOTS       01
#define FLAG_VOLUME     02
#define FLAG_ATOM_AREA  04
namespace Almost {
  //Solvate
  void boxsolvate(const Protein & p, const Solvent & wat,
		  double t, string file){
    //Compute BB of protein
    if(p.atom_size()==0){
      cout<<"\t>> Protein is empty: nothing to solvate!!"<<endl;
      return;
    }

    double max[3],min[3],wmax[3],wmin[3];
    max[0] = p[0][0];
    max[1] = p[0][1];
    max[2] = p[0][2];
    min[0] = p[0][0];
    min[1] = p[0][1];
    min[2] = p[0][2];

    for(int i=0;i<p.atom_size();i++){
      if(p[i][0]>max[0]) max[0] = p[i][0];
      if(p[i][1]>max[1]) max[1] = p[i][1];
      if(p[i][2]>max[2]) max[2] = p[i][2];
      if(p[i][0]<min[0]) min[0] = p[i][0];
      if(p[i][1]<min[1]) min[1] = p[i][1];
      if(p[i][2]<min[2]) min[2] = p[i][2];
    }

    for(int i=0;i<3;i++) max[i] += t;
    for(int i=0;i<3;i++) min[i] -= t;

    cout<<"\t>> Generating box of size:\n"
	<<"\t         "<<max[0]-min[0]<<"x"
	<<max[1]-min[1]<<"x"
	<<max[2]-min[2]<<"\n";

    wmax[0] = wat[0][0];
    wmax[1] = wat[0][1];
    wmax[2] = wat[0][2];
    wmin[0] = wat[0][0];
    wmin[1] = wat[0][1];
    wmin[2] = wat[0][2];

    for(int i=0;i<wat.atom_size();i++){
      if(wat[i][0]>wmax[0]) wmax[0] = wat[i][0];
      if(wat[i][1]>wmax[1]) wmax[1] = wat[i][1];
      if(wat[i][2]>wmax[2]) wmax[2] = wat[i][2];
      if(wat[i][0]<wmin[0]) wmin[0] = wat[i][0];
      if(wat[i][1]<wmin[1]) wmin[1] = wat[i][1];
      if(wat[i][2]<wmin[2]) wmin[2] = wat[i][2];
    }


    cout<<"\t>> Water box of size:\n"
	<<"\t         "<<wmax[0]-wmin[0]<<"x"
	<<wmax[1]-wmin[1]<<"x"
	<<wmax[2]-wmin[2]<<"\n";
  

    Solvent wtrans = wat;
    Coor3D tr;
    tr[0] = min[0] - wmin[0];
    tr[1] = min[1] - wmin[1];
    tr[2] = min[2] - wmin[2];
    wtrans.translate(tr);

    {
      ofstream f;
      f.open("deb1.pdb");
      PDBOStream pdb(f);
      pdb<<wtrans;
    }


    int Nx = (int)(ceil((max[0]-min[0])/(wmax[0]-wmin[0])));
    int Ny = (int)(ceil((max[1]-min[1])/(wmax[1]-wmin[1])));
    int Nz = (int)(ceil((max[2]-min[2])/(wmax[2]-wmin[2])));

    cout<<"\t>> Construncting  water grid of size: \n"
	<<"\t      "<<Nx<<"x"<<Ny<<"x"<<Nz<<endl;

    double dx = wmax[0]-wmin[0];
    double dy = wmax[1]-wmin[1];
    double dz = wmax[2]-wmin[2];

    Molecules tip;
    for(int i=0;i<Nx;i++){
      for(int j=0;j<Ny;j++){
	for(int k=0;k<Nz;k++){
	  Coor3D tplus;
	  tplus[0] = dx*i;
	  tplus[1] = dy*j;
	  tplus[2] = dz*k;

	  wtrans.translate(tplus);
	  tip.add_solvent(wtrans);
	  Coor3D tminus;
	  tminus[0] = -dx*i;
	  tminus[1] = -dy*j;
	  tminus[2] = -dz*k;
	
	  wtrans.translate(tminus);
	}
      }
    }

    //remove unnecessary waters
    int nwat = tip.fragment_size();
    if(tip.atom_size()!=3*nwat){
      cout<<"\t>> FATAL ERROR: WATER is not of type TIP3P!\n";
      return;
    }

    vector<int> remove;
    remove.resize(nwat);
    int atom_counter = 0;
    for(int i=0;i<nwat;i++){
      //check in is in box
      for(int a = atom_counter;a<atom_counter+3;a++){
	if((tip[a][0]<min[0]+1)||
	   (tip[a][1]<min[1]+1)||
	   (tip[a][2]<min[2]+1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
	//       else {
	// 	cout<<"KEEP "<<a<<" "<<tip[a][0]<<" "<<min[0]<<endl;
	//       }

	if((tip[a][0]>max[0]-1)||
	   (tip[a][1]>max[1]-1)||
	   (tip[a][2]>max[2]-1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
      }
      atom_counter += 3;    
    }

    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){ atom_counter += 3;continue;}
      for(int j = 0;j<p.atom_size();j++){
	double dx = p[j][0]-tip[atom_counter][0];
	double dy = p[j][1]-tip[atom_counter][1];
	double dz = p[j][2]-tip[atom_counter][2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2<(2.5*2.5)){
	  remove[i] = 1;
	  break;
	}
      }
      atom_counter += 3;
    }

    //writr PDB
    ofstream pdb;
    pdb.open(file.c_str());

    atom_counter = 0;
    int natom = 0;
    int nfrag = 0;
    int removed = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){removed++; atom_counter += 3;continue;}
      char buff[180];
      for(int l=0;l<3;l++){
	sprintf(buff,
		"%-6s%5i %-4s%5s%5i    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s\n",
		"ATOM",
		++natom,
		tip[atom_counter].name().c_str(),
		"TIP3",
		nfrag+1,
		tip[atom_counter].x(),
		tip[atom_counter].y(),
		tip[atom_counter].z(),
		(double)0,
		(double)0,
		"WAT"
		);
	pdb<<buff;
	++atom_counter;
      }
      ++nfrag;
    }
    pdb<<"END\n";
    cout<<"\t>> "<<removed<<" water molecules removed because of overlaps\n";
  }



  double spheresolvate(const Protein & p, const Solvent & wat,
		       double t, string file, const MDB& mdb ){
    //Compute BB of protein
    if(p.atom_size()==0){
      cout<<"\t>> Protein is empty: nothing to solvate!!"<<endl;
      return 0;
    }


    Coor3D com = p.com(mdb);
    //compute spherical radius
    double R2 = 0;
    double R = 0;
    for(int i=0;i<p.atom_size();i++){
      double dx,dy,dz;
      dx = p[i][0]-com[0];
      dy = p[i][1]-com[1];
      dz = p[i][2]-com[2];
      double d2 = dx*dx+dy*dy+dz*dz;
      if(d2>R2) R2=d2;
    }
    R = sqrt(R2)+t;
    R2 = R*R;



    double max[3],min[3],wmax[3],wmin[3];
    max[0] = com[0]+R;
    max[1] = com[1]+R;
    max[2] = com[2]+R;
    min[0] = com[0]-R;
    min[1] = com[1]-R;
    min[2] = com[2]-R;



  
    cout<<"\t>> Generating box on size:\n"
	<<"\t         "<<max[0]-min[0]<<"x"
	<<max[1]-min[1]<<"x"
	<<max[2]-min[2]<<"\n";

    wmax[0] = wat[0][0];
    wmax[1] = wat[0][1];
    wmax[2] = wat[0][2];
    wmin[0] = wat[0][0];
    wmin[1] = wat[0][1];
    wmin[2] = wat[0][2];

    for(int i=0;i<wat.atom_size();i++){
      if(wat[i][0]>wmax[0]) wmax[0] = wat[i][0];
      if(wat[i][1]>wmax[1]) wmax[1] = wat[i][1];
      if(wat[i][2]>wmax[2]) wmax[2] = wat[i][2];
      if(wat[i][0]<wmin[0]) wmin[0] = wat[i][0];
      if(wat[i][1]<wmin[1]) wmin[1] = wat[i][1];
      if(wat[i][2]<wmin[2]) wmin[2] = wat[i][2];
    }


    cout<<"\t>> Water box on size:\n"
	<<"\t         "<<wmax[0]-wmin[0]<<"x"
	<<wmax[1]-wmin[1]<<"x"
	<<wmax[2]-wmin[2]<<"\n";
  

    Solvent wtrans = wat;
    Coor3D tr;
    tr[0] = min[0] - wmin[0];
    tr[1] = min[1] - wmin[1];
    tr[2] = min[2] - wmin[2];
    wtrans.translate(tr);

    {
      ofstream f;
      f.open("deb1.pdb");
      PDBOStream pdb(f);
      pdb<<wtrans;
    }


    int Nx = (int)(ceil((max[0]-min[0])/(wmax[0]-wmin[0])));
    int Ny = (int)(ceil((max[1]-min[1])/(wmax[1]-wmin[1])));
    int Nz = (int)(ceil((max[2]-min[2])/(wmax[2]-wmin[2])));

    cout<<"\t>> Construncting  water grid of size: \n"
	<<"\t         "<<Nx<<"x"<<Ny<<"x"<<Nz<<endl;

    double dx = wmax[0]-wmin[0];
    double dy = wmax[1]-wmin[1];
    double dz = wmax[2]-wmin[2];

    Molecules tip;
    for(int i=0;i<Nx;i++){
      for(int j=0;j<Ny;j++){
	for(int k=0;k<Nz;k++){
	  Coor3D tplus;
	  tplus[0] = dx*i;
	  tplus[1] = dy*j;
	  tplus[2] = dz*k;

	  wtrans.translate(tplus);
	  tip.add_solvent(wtrans);
	  Coor3D tminus;
	  tminus[0] = -dx*i;
	  tminus[1] = -dy*j;
	  tminus[2] = -dz*k;
	
	  wtrans.translate(tminus);
	}
      }
    }

    //remove unnecessary waters
    int nwat = tip.fragment_size();
    if(tip.atom_size()!=3*nwat){
      cout<<"\t>> FATAL ERROR: WATER is not of type TIP3P!\n";
      return 0;
    }

    vector<int> remove;
    remove.resize(nwat);
    int atom_counter = 0;
    for(int i=0;i<nwat;i++){
      //check in is in box
      for(int a = atom_counter;a<atom_counter+3;a++){
	if((tip[a][0]<min[0]+1)||
	   (tip[a][1]<min[1]+1)||
	   (tip[a][2]<min[2]+1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
	//       else {
	// 	cout<<"KEEP "<<a<<" "<<tip[a][0]<<" "<<min[0]<<endl;
	//       }

	if((tip[a][0]>max[0]-1)||
	   (tip[a][1]>max[1]-1)||
	   (tip[a][2]>max[2]-1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
      }
      atom_counter += 3;    
    }


    //remove sphere
    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      for(int a = atom_counter;a<atom_counter+3;a++){
	double dx,dy,dz;
	dx = tip[a][0]-com[0];
	dy = tip[a][1]-com[1];
	dz = tip[a][2]-com[2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2>R2){ remove[i] = 1; break;}
      }
      atom_counter +=3;
    }


    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){ atom_counter += 3;continue;}
      for(int j = 0;j<p.atom_size();j++){
	double dx = p[j][0]-tip[atom_counter][0];
	double dy = p[j][1]-tip[atom_counter][1];
	double dz = p[j][2]-tip[atom_counter][2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2<(2.5*2.5)){
	  remove[i] = 1;
	  break;
	}
      }
      atom_counter += 3;
    }



    //writr PDB
    ofstream pdb;
    pdb.open(file.c_str());

    atom_counter = 0;
    int natom = 0;
    int nfrag = 0;
    int removed = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){removed++; atom_counter += 3;continue;}
      char buff[180];
      for(int l=0;l<3;l++){
	sprintf(buff,
		"%-6s%5i %-4s%5s%5i    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s\n",
		"ATOM",
		++natom,
		tip[atom_counter].name().c_str(),
		"TIP3",
		nfrag+1,
		tip[atom_counter].x(),
		tip[atom_counter].y(),
		tip[atom_counter].z(),
		(double)0,
		(double)0,
		"WAT"
		);
	pdb<<buff;
	++atom_counter;
      }
      ++nfrag;
    }
    pdb<<"END\n";
    return sqrt(R2);
  }


  void shellsolvate(const Protein & p, const Solvent & wat,
		    double t, string file, const MDB& mdb ){
    //Compute BB of protein
    if(p.atom_size()==0){
      cout<<"\t>> Protein is empty: nothing to solvate!!"<<endl;
      return;
    }


    Coor3D com = p.com(mdb);
    //compute spherical radius
    double R2 = 0;
    double R = 0;
    for(int i=0;i<p.atom_size();i++){
      double dx,dy,dz;
      dx = p[i][0]-com[0];
      dy = p[i][1]-com[1];
      dz = p[i][2]-com[2];
      double d2 = dx*dx+dy*dy+dz*dz;
      if(d2>R2) R2=d2;
    }
    R = sqrt(R2)+t;
    R2 = R*R;



    double max[3],min[3],wmax[3],wmin[3];
    max[0] = com[0]+R;
    max[1] = com[1]+R;
    max[2] = com[2]+R;
    min[0] = com[0]-R;
    min[1] = com[1]-R;
    min[2] = com[2]-R;



  
    cout<<"\t>> Generating box on size:\n"
	<<"\t         "<<max[0]-min[0]<<"x"
	<<max[1]-min[1]<<"x"
	<<max[2]-min[2]<<"\n";

    wmax[0] = wat[0][0];
    wmax[1] = wat[0][1];
    wmax[2] = wat[0][2];
    wmin[0] = wat[0][0];
    wmin[1] = wat[0][1];
    wmin[2] = wat[0][2];

    for(int i=0;i<wat.atom_size();i++){
      if(wat[i][0]>wmax[0]) wmax[0] = wat[i][0];
      if(wat[i][1]>wmax[1]) wmax[1] = wat[i][1];
      if(wat[i][2]>wmax[2]) wmax[2] = wat[i][2];
      if(wat[i][0]<wmin[0]) wmin[0] = wat[i][0];
      if(wat[i][1]<wmin[1]) wmin[1] = wat[i][1];
      if(wat[i][2]<wmin[2]) wmin[2] = wat[i][2];
    }


    cout<<"\t>> Water box on size:\n"
	<<"\t         "<<wmax[0]-wmin[0]<<"x"
	<<wmax[1]-wmin[1]<<"x"
	<<wmax[2]-wmin[2]<<"\n";
  

    Solvent wtrans = wat;
    Coor3D tr;
    tr[0] = min[0] - wmin[0];
    tr[1] = min[1] - wmin[1];
    tr[2] = min[2] - wmin[2];
    wtrans.translate(tr);

    {
      ofstream f;
      f.open("deb1.pdb");
      PDBOStream pdb(f);
      pdb<<wtrans;
    }


    int Nx = (int)(ceil((max[0]-min[0])/(wmax[0]-wmin[0])));
    int Ny = (int)(ceil((max[1]-min[1])/(wmax[1]-wmin[1])));
    int Nz = (int)(ceil((max[2]-min[2])/(wmax[2]-wmin[2])));

    cout<<"\t>> Construncting  water grid of size: \n"
	<<"\t         "<<Nx<<"x"<<Ny<<"x"<<Nz<<endl;

    double dx = wmax[0]-wmin[0];
    double dy = wmax[1]-wmin[1];
    double dz = wmax[2]-wmin[2];

    Molecules tip;
    for(int i=0;i<Nx;i++){
      for(int j=0;j<Ny;j++){
	for(int k=0;k<Nz;k++){
	  Coor3D tplus;
	  tplus[0] = dx*i;
	  tplus[1] = dy*j;
	  tplus[2] = dz*k;

	  wtrans.translate(tplus);
	  tip.add_solvent(wtrans);
	  Coor3D tminus;
	  tminus[0] = -dx*i;
	  tminus[1] = -dy*j;
	  tminus[2] = -dz*k;
	
	  wtrans.translate(tminus);
	}
      }
    }

    //remove unnecessary waters
    int nwat = tip.fragment_size();
    if(tip.atom_size()!=3*nwat){
      cout<<"\t>> FATAL ERROR: WATER is not of type TIP3P!\n";
      return;
    }

    vector<int> remove;
    remove.resize(nwat);
    int atom_counter = 0;
    for(int i=0;i<nwat;i++){
      //check in is in box
      for(int a = atom_counter;a<atom_counter+3;a++){
	if((tip[a][0]<min[0]+1)||
	   (tip[a][1]<min[1]+1)||
	   (tip[a][2]<min[2]+1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
	//       else {
	// 	cout<<"KEEP "<<a<<" "<<tip[a][0]<<" "<<min[0]<<endl;
	//       }

	if((tip[a][0]>max[0]-1)||
	   (tip[a][1]>max[1]-1)||
	   (tip[a][2]>max[2]-1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
      }
      atom_counter += 3;    
    }


    //remove sphere
    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      for(int a = atom_counter;a<atom_counter+3;a++){
	double dx,dy,dz;
	dx = tip[a][0]-com[0];
	dy = tip[a][1]-com[1];
	dz = tip[a][2]-com[2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2>R2){ remove[i] = 1; break;}
      }
      atom_counter +=3;
    }


    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){ atom_counter += 3;continue;}
      double m = 10000;
      for(int j = 0;j<p.atom_size();j++){
	double dx = p[j][0]-tip[atom_counter][0];
	double dy = p[j][1]-tip[atom_counter][1];
	double dz = p[j][2]-tip[atom_counter][2];
	double d2 = dx*dx+dy*dy+dz*dz;
	double d = sqrt(d2);
	if(d2<(2.5*2.5)){
	  remove[i] = 1;
	  break;
	}
	NBondT bnt = mdb.nbondT(p[j].type_name());
	double vdw = bnt.sigma();
	if(d-vdw-t < m){
	  m = d-vdw-t;
	}
      }
      if(m>0){ remove[i] = 1;}
      atom_counter += 3;
    }



    //writr PDB
    ofstream pdb;
    pdb.open(file.c_str());

    atom_counter = 0;
    int natom = 0;
    int nfrag = 0;
    int removed = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){removed++; atom_counter += 3;continue;}
      char buff[180];
      for(int l=0;l<3;l++){
	sprintf(buff,
		"%-6s%5i %-4s%5s%5i    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s\n",
		"ATOM",
		++natom,
		tip[atom_counter].name().c_str(),
		"TIP3",
		nfrag+1,
		tip[atom_counter].x(),
		tip[atom_counter].y(),
		tip[atom_counter].z(),
		(double)0,
		(double)0,
		"WAT"
		);
	pdb<<buff;
	++atom_counter;
      }
      ++nfrag;
    }
    pdb<<"END\n";
  }

  void molshellsolvate(const Molecules & p, const Solvent & wat,
		       double t, string file, const MDB& mdb ){
    //Compute BB of protein
    if(p.atom_size()==0){
      cout<<"\t>> Moecules are empty: nothing to solvate!!"<<endl;
      return;
    }


    Coor3D com = p.com(mdb);
    //compute spherical radius
    double R2 = 0;
    double R = 0;
    for(int i=0;i<p.atom_size();i++){
      double dx,dy,dz;
      dx = p[i][0]-com[0];
      dy = p[i][1]-com[1];
      dz = p[i][2]-com[2];
      double d2 = dx*dx+dy*dy+dz*dz;
      if(d2>R2) R2=d2;
    }
    R = sqrt(R2)+t;
    R2 = R*R;



    double max[3],min[3],wmax[3],wmin[3];
    max[0] = com[0]+R;
    max[1] = com[1]+R;
    max[2] = com[2]+R;
    min[0] = com[0]-R;
    min[1] = com[1]-R;
    min[2] = com[2]-R;



  
    cout<<"\t>> Generating box on size:\n"
	<<"\t         "<<max[0]-min[0]<<"x"
	<<max[1]-min[1]<<"x"
	<<max[2]-min[2]<<"\n";

    wmax[0] = wat[0][0];
    wmax[1] = wat[0][1];
    wmax[2] = wat[0][2];
    wmin[0] = wat[0][0];
    wmin[1] = wat[0][1];
    wmin[2] = wat[0][2];

    for(int i=0;i<wat.atom_size();i++){
      if(wat[i][0]>wmax[0]) wmax[0] = wat[i][0];
      if(wat[i][1]>wmax[1]) wmax[1] = wat[i][1];
      if(wat[i][2]>wmax[2]) wmax[2] = wat[i][2];
      if(wat[i][0]<wmin[0]) wmin[0] = wat[i][0];
      if(wat[i][1]<wmin[1]) wmin[1] = wat[i][1];
      if(wat[i][2]<wmin[2]) wmin[2] = wat[i][2];
    }


    cout<<"\t>> Water box on size:\n"
	<<"\t         "<<wmax[0]-wmin[0]<<"x"
	<<wmax[1]-wmin[1]<<"x"
	<<wmax[2]-wmin[2]<<"\n";
  

    Solvent wtrans = wat;
    Coor3D tr;
    tr[0] = min[0] - wmin[0];
    tr[1] = min[1] - wmin[1];
    tr[2] = min[2] - wmin[2];
    wtrans.translate(tr);

    {
      ofstream f;
      f.open("deb1.pdb");
      PDBOStream pdb(f);
      pdb<<wtrans;
    }


    int Nx = (int)(ceil((max[0]-min[0])/(wmax[0]-wmin[0])));
    int Ny = (int)(ceil((max[1]-min[1])/(wmax[1]-wmin[1])));
    int Nz = (int)(ceil((max[2]-min[2])/(wmax[2]-wmin[2])));

    cout<<"\t>> Construncting  water grid of size: \n"
	<<"\t         "<<Nx<<"x"<<Ny<<"x"<<Nz<<endl;

    double dx = wmax[0]-wmin[0];
    double dy = wmax[1]-wmin[1];
    double dz = wmax[2]-wmin[2];

    Molecules tip;
    for(int i=0;i<Nx;i++){
      for(int j=0;j<Ny;j++){
	for(int k=0;k<Nz;k++){
	  Coor3D tplus;
	  tplus[0] = dx*i;
	  tplus[1] = dy*j;
	  tplus[2] = dz*k;

	  wtrans.translate(tplus);
	  tip.add_solvent(wtrans);
	  Coor3D tminus;
	  tminus[0] = -dx*i;
	  tminus[1] = -dy*j;
	  tminus[2] = -dz*k;
	
	  wtrans.translate(tminus);
	}
      }
    }

    //remove unnecessary waters
    int nwat = tip.fragment_size();
    if(tip.atom_size()!=3*nwat){
      cout<<"\t>> FATAL ERROR: WATER is not of type TIP3P!\n";
      return;
    }

    vector<int> remove;
    remove.resize(nwat);
    int atom_counter = 0;
    for(int i=0;i<nwat;i++){
      //check in is in box
      for(int a = atom_counter;a<atom_counter+3;a++){
	if((tip[a][0]<min[0]+1)||
	   (tip[a][1]<min[1]+1)||
	   (tip[a][2]<min[2]+1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
	//       else {
	// 	cout<<"KEEP "<<a<<" "<<tip[a][0]<<" "<<min[0]<<endl;
	//       }

	if((tip[a][0]>max[0]-1)||
	   (tip[a][1]>max[1]-1)||
	   (tip[a][2]>max[2]-1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
      }
      atom_counter += 3;    
    }


    //remove sphere
    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      for(int a = atom_counter;a<atom_counter+3;a++){
	double dx,dy,dz;
	dx = tip[a][0]-com[0];
	dy = tip[a][1]-com[1];
	dz = tip[a][2]-com[2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2>R2){ remove[i] = 1; break;}
      }
      atom_counter +=3;
    }


    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){ atom_counter += 3;continue;}
      double m = 10000;
      for(int j = 0;j<p.atom_size();j++){
	double dx = p[j][0]-tip[atom_counter][0];
	double dy = p[j][1]-tip[atom_counter][1];
	double dz = p[j][2]-tip[atom_counter][2];
	double d2 = dx*dx+dy*dy+dz*dz;
	double d = sqrt(d2);
	if(d2<(2.5*2.5)){
	  remove[i] = 1;
	  break;
	}
	NBondT bnt = mdb.nbondT(p[j].type_name());
	double vdw = bnt.sigma();
	if(d-vdw-t < m){
	  m = d-vdw-t;
	}
      }
      if(m>0){ remove[i] = 1;}
      atom_counter += 3;
    }



    //writr PDB
    ofstream pdb;
    pdb.open(file.c_str());

    atom_counter = 0;
    int natom = 0;
    int nfrag = 0;
    int removed = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){removed++; atom_counter += 3;continue;}
      char buff[180];
      for(int l=0;l<3;l++){
	sprintf(buff,
		"%-6s%5i %-4s%5s%5i    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s\n",
		"ATOM",
		++natom,
		tip[atom_counter].name().c_str(),
		"TIP3",
		nfrag+1,
		tip[atom_counter].x(),
		tip[atom_counter].y(),
		tip[atom_counter].z(),
		(double)0,
		(double)0,
		"WAT"
		);
	pdb<<buff;
	++atom_counter;
      }
      ++nfrag;
    }
    pdb<<"END\n";
  }

  double molspheresolvate(const Molecules & mols, const Solvent & wat,
			  double t, string file, const MDB& mdb ){
    //Compute BB of protein
    if(mols.atom_size()==0){
      cout<<"\t>> Moleculs empty: nothing to solvate!!"<<endl;
      return 0;
    }


    Coor3D com = mols.com(mdb);
    //compute spherical radius
    double R2 = 0;
    double R = 0;
    for(int i=0;i<mols.atom_size();i++){
      double dx,dy,dz;
      dx = mols[i][0]-com[0];
      dy = mols[i][1]-com[1];
      dz = mols[i][2]-com[2];
      double d2 = dx*dx+dy*dy+dz*dz;
      if(d2>R2) R2=d2;
    }
    R = sqrt(R2)+t;
    R2 = R*R;



    double max[3],min[3],wmax[3],wmin[3];
    max[0] = com[0]+R;
    max[1] = com[1]+R;
    max[2] = com[2]+R;
    min[0] = com[0]-R;
    min[1] = com[1]-R;
    min[2] = com[2]-R;



  
    cout<<"\t>> Generating box on size:\n"
	<<"\t         "<<max[0]-min[0]<<"x"
	<<max[1]-min[1]<<"x"
	<<max[2]-min[2]<<"\n";

    wmax[0] = wat[0][0];
    wmax[1] = wat[0][1];
    wmax[2] = wat[0][2];
    wmin[0] = wat[0][0];
    wmin[1] = wat[0][1];
    wmin[2] = wat[0][2];

    for(int i=0;i<wat.atom_size();i++){
      if(wat[i][0]>wmax[0]) wmax[0] = wat[i][0];
      if(wat[i][1]>wmax[1]) wmax[1] = wat[i][1];
      if(wat[i][2]>wmax[2]) wmax[2] = wat[i][2];
      if(wat[i][0]<wmin[0]) wmin[0] = wat[i][0];
      if(wat[i][1]<wmin[1]) wmin[1] = wat[i][1];
      if(wat[i][2]<wmin[2]) wmin[2] = wat[i][2];
    }


    cout<<"\t>> Water box on size:\n"
	<<"\t         "<<wmax[0]-wmin[0]<<"x"
	<<wmax[1]-wmin[1]<<"x"
	<<wmax[2]-wmin[2]<<"\n";
  

    Solvent wtrans = wat;
    Coor3D tr;
    tr[0] = min[0] - wmin[0];
    tr[1] = min[1] - wmin[1];
    tr[2] = min[2] - wmin[2];
    wtrans.translate(tr);

    {
      ofstream f;
      f.open("deb1.pdb");
      PDBOStream pdb(f);
      pdb<<wtrans;
    }


    int Nx = (int)(ceil((max[0]-min[0])/(wmax[0]-wmin[0])));
    int Ny = (int)(ceil((max[1]-min[1])/(wmax[1]-wmin[1])));
    int Nz = (int)(ceil((max[2]-min[2])/(wmax[2]-wmin[2])));

    cout<<"\t>> Construncting  water grid of size: \n"
	<<"\t         "<<Nx<<"x"<<Ny<<"x"<<Nz<<endl;

    double dx = wmax[0]-wmin[0];
    double dy = wmax[1]-wmin[1];
    double dz = wmax[2]-wmin[2];

    Molecules tip;
    for(int i=0;i<Nx;i++){
      for(int j=0;j<Ny;j++){
	for(int k=0;k<Nz;k++){
	  Coor3D tplus;
	  tplus[0] = dx*i;
	  tplus[1] = dy*j;
	  tplus[2] = dz*k;

	  wtrans.translate(tplus);
	  tip.add_solvent(wtrans);
	  Coor3D tminus;
	  tminus[0] = -dx*i;
	  tminus[1] = -dy*j;
	  tminus[2] = -dz*k;
	
	  wtrans.translate(tminus);
	}
      }
    }

    //remove unnecessary waters
    int nwat = tip.fragment_size();
    if(tip.atom_size()!=3*nwat){
      cout<<"\t>> FATAL ERROR: WATER is not of type TIP3P!\n";
      return 0;
    }

    vector<int> remove;
    remove.resize(nwat);
    int atom_counter = 0;
    for(int i=0;i<nwat;i++){
      //check in is in box
      for(int a = atom_counter;a<atom_counter+3;a++){
	if((tip[a][0]<min[0]+1)||
	   (tip[a][1]<min[1]+1)||
	   (tip[a][2]<min[2]+1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
	//       else {
	// 	cout<<"KEEP "<<a<<" "<<tip[a][0]<<" "<<min[0]<<endl;
	//       }

	if((tip[a][0]>max[0]-1)||
	   (tip[a][1]>max[1]-1)||
	   (tip[a][2]>max[2]-1)){
	  remove[i] = 1;
	  //	cout<<"BOX size exceded "<<endl;
	}
      }
      atom_counter += 3;    
    }


    //remove sphere
    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      for(int a = atom_counter;a<atom_counter+3;a++){
	double dx,dy,dz;
	dx = tip[a][0]-com[0];
	dy = tip[a][1]-com[1];
	dz = tip[a][2]-com[2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2>R2){ remove[i] = 1; break;}
      }
      atom_counter +=3;
    }


    atom_counter = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){ atom_counter += 3;continue;}
      for(int j = 0;j<mols.atom_size();j++){
	double dx = mols[j][0]-tip[atom_counter][0];
	double dy = mols[j][1]-tip[atom_counter][1];
	double dz = mols[j][2]-tip[atom_counter][2];
	double d2 = dx*dx+dy*dy+dz*dz;
	if(d2<(2.5*2.5)){
	  remove[i] = 1;
	  break;
	}
      }
      atom_counter += 3;
    }



    //writr PDB
    ofstream pdb;
    pdb.open(file.c_str());

    atom_counter = 0;
    int natom = 0;
    int nfrag = 0;
    int removed = 0;
    for(int i=0;i<nwat;i++){
      if(remove[i]==1){removed++; atom_counter += 3;continue;}
      char buff[180];
      for(int l=0;l<3;l++){
	sprintf(buff,
		"%-6s%5i %-4s%5s%5i    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s\n",
		"ATOM",
		++natom,
		tip[atom_counter].name().c_str(),
		"TIP3",
		nfrag+1,
		tip[atom_counter].x(),
		tip[atom_counter].y(),
		tip[atom_counter].z(),
		(double)0,
		(double)0,
		"WAT"
		);
	pdb<<buff;
	++atom_counter;
      }
      ++nfrag;
    }
    pdb<<"END\n";
    return sqrt(R2);
  }

  double fit(const Molecules & m1, Molecules & m2, const NameSelection &s){

    int size;
    double r[3][3],cog[3],cog1[3];
    Coor<double> c1(m1);
    Coor<double> c2(m2);
    double rmsd = rmsd_mat(c1,c2,s.result(),r);
    //check selection !!!!
    //now rotate
    for(int i=0;i<3;i++) cog[i] = cog1[i] = 0; 
    set<int>::iterator iter = s.result().begin(),
      end = s.result().end();
    size = s.result().size();
    while(iter!=end){
      int pos = (*iter)*Coor<double>::DIM;
      cog[0] += c2.coor[pos];
      cog[1] += c2.coor[pos+1];
      cog[2] += c2.coor[pos+2];
      cog1[0] += c1.coor[pos];
      cog1[1] += c1.coor[pos+1];
      cog1[2] += c1.coor[pos+2];
      ++iter;
    }

    for(int i=0;i<3;i++) cog[i] /= size;
    for(int i=0;i<3;i++) cog1[i] /= size;

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;
      c2.coor[pos] -= cog[0];
      c2.coor[pos+1] -= cog[1];
      c2.coor[pos+2] -= cog[2];
    }

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;

      c1.coor[pos] = 
	r[0][0]*c2.coor[pos] +
	r[0][1]*c2.coor[pos+1] +
	r[0][2]*c2.coor[pos+2];

      c1.coor[pos+1] = 
	r[1][0]*c2.coor[pos] +
	r[1][1]*c2.coor[pos+1] +
	r[1][2]*c2.coor[pos+2];

      c1.coor[pos+2] = 
	r[2][0]*c2.coor[pos] +
	r[2][1]*c2.coor[pos+1] +
	r[2][2]*c2.coor[pos+2];
    }

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;
      c1.coor[pos] += cog1[0];
      c1.coor[pos+1] += cog1[1];
      c1.coor[pos+2] += cog1[2];
    }

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;
      //assign
      m2[i][0] = c1.coor[pos];
      m2[i][1] = c1.coor[pos+1];
      m2[i][2] = c1.coor[pos+2];
    }

    return rmsd;
  }


  double fit2(const Molecules & m1, Molecules & m2, 
	      const NameSelection &s1,
	      const NameSelection &s2){

    int size;
    double r[3][3],cog[3],cog1[3];
    Coor<double> c1(m1);
    Coor<double> c2(m2);
    if(s1.result().size()!=s2.result().size()){
      cout<<"\t>> Can't fit molecules with selections of different size "
	  <<s1.result().size()<<" "<<s2.result().size()<<"\n";
      return -1;
    }
    double rmsd = rmsd2_mat(c1,c2,s1.result(),s2.result(),r);
  
    //check selection !!!!
    //now rotate
    for(int i=0;i<3;i++) cog[i] = cog1[i] = 0; 
    set<int>::iterator iter = s2.result().begin(),
      end = s2.result().end();
    size = s1.result().size();
    while(iter!=end){
      int pos = (*iter)*Coor<double>::DIM;
      cog[0] += c2.coor[pos];
      cog[1] += c2.coor[pos+1];
      cog[2] += c2.coor[pos+2];
      ++iter;
    }

    iter = s1.result().begin();
    end = s1.result().end();
    while(iter!=end){
      int pos = (*iter)*Coor<double>::DIM;
      cog1[0] += c1.coor[pos];
      cog1[1] += c1.coor[pos+1];
      cog1[2] += c1.coor[pos+2];
      ++iter;
    }

    for(int i=0;i<3;i++) cog[i] /= size;
    for(int i=0;i<3;i++) cog1[i] /= size;

    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;
      c2.coor[pos] -= cog[0];
      c2.coor[pos+1] -= cog[1];
      c2.coor[pos+2] -= cog[2];
    }

    c1.resize(c2.size());
    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;

      c1.coor[pos] = 
	r[0][0]*c2.coor[pos] +
	r[0][1]*c2.coor[pos+1] +
	r[0][2]*c2.coor[pos+2];

      c1.coor[pos+1] = 
	r[1][0]*c2.coor[pos] +
	r[1][1]*c2.coor[pos+1] +
	r[1][2]*c2.coor[pos+2];

      c1.coor[pos+2] = 
	r[2][0]*c2.coor[pos] +
	r[2][1]*c2.coor[pos+1] +
	r[2][2]*c2.coor[pos+2];
    }

    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;
      c1.coor[pos] += cog1[0];
      c1.coor[pos+1] += cog1[1];
      c1.coor[pos+2] += cog1[2];
    }

    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;
      //assign
      m2[i][0] = c1.coor[pos];
      m2[i][1] = c1.coor[pos+1];
      m2[i][2] = c1.coor[pos+2];
    }

    return rmsd;
  }


  double fit2_coor(const Coor<double> & c1, Coor<double> & c2, 
		   const NameSelection &s1,
		   const NameSelection &s2){

    int size;
    double r[3][3],cog[3],cog1[3];
    //    Coor<double> c1(m1);
    //    Coor<double> c2(m2);
    if(s1.result().size()!=s2.result().size()){
      cout<<"\t>> Can't fit molecules with selections of different size "
	  <<s1.result().size()<<" "<<s2.result().size()<<"\n";
      return -1;
    }
    double rmsd = rmsd2_mat(c1,c2,s1.result(),s2.result(),r);
  
    //check selection !!!!
    //now rotate
    for(int i=0;i<3;i++) cog[i] = cog1[i] = 0; 
    set<int>::iterator iter = s2.result().begin(),
      end = s2.result().end();
    size = s1.result().size();
    while(iter!=end){
      int pos = (*iter)*Coor<double>::DIM;
      cog[0] += c2.coor[pos];
      cog[1] += c2.coor[pos+1];
      cog[2] += c2.coor[pos+2];
      ++iter;
    }

    iter = s1.result().begin();
    end = s1.result().end();
    while(iter!=end){
      int pos = (*iter)*Coor<double>::DIM;
      cog1[0] += c1.coor[pos];
      cog1[1] += c1.coor[pos+1];
      cog1[2] += c1.coor[pos+2];
      ++iter;
    }

    for(int i=0;i<3;i++) cog[i] /= size;
    for(int i=0;i<3;i++) cog1[i] /= size;

    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;
      c2.coor[pos] -= cog[0];
      c2.coor[pos+1] -= cog[1];
      c2.coor[pos+2] -= cog[2];
    }

    Coor<double> c_tmp;
    c_tmp.resize(c2.size());
    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;

      c_tmp.coor[pos] = 
	r[0][0]*c2.coor[pos] +
	r[0][1]*c2.coor[pos+1] +
	r[0][2]*c2.coor[pos+2];

      c_tmp.coor[pos+1] = 
	r[1][0]*c2.coor[pos] +
	r[1][1]*c2.coor[pos+1] +
	r[1][2]*c2.coor[pos+2];

      c_tmp.coor[pos+2] = 
	r[2][0]*c2.coor[pos] +
	r[2][1]*c2.coor[pos+1] +
	r[2][2]*c2.coor[pos+2];
    }

    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;
      c_tmp.coor[pos] += cog1[0];
      c_tmp.coor[pos+1] += cog1[1];
      c_tmp.coor[pos+2] += cog1[2];
    }

    for(int i=0;i<c2.size();i++){
      int pos = i*Coor<double>::DIM;
      //assign
      c2.coor[pos] = c_tmp.coor[pos];
      c2.coor[pos+1] = c_tmp.coor[pos+1];
      c2.coor[pos+2] = c_tmp.coor[pos+2];
    }
    return rmsd;
  }


  double mfit(const Molecules & m1, 
	      Molecules & m2, 
	      const NameSelection &s,
	      const MDB & mdb){

    double r[3][3],com[3],com1[3];
    int size = s.result().size();
    double * mass = new double[size];
    double *mass_sqrt = new double[size];
    double M;
    {
      set<int>::iterator iter = s.result().begin(),
	end = s.result().end();
      int count = 0;
      M = 0;
      while(iter!=end){
	mass[count] = mdb.atomT(m1[*iter].type_name()).mass(); 
	M += mass[count];
	++count;
	++iter;
      }
      iter = s.result().begin();
      count = 0;
      while(iter!=end){
	mass_sqrt[count] = sqrt(mass[count]);
	++count;
	++iter;
      }
    }
    Coor<double> c1(m1);
    Coor<double> c2(m2);

    //check selection !!!!
    //now rotate
    for(int i=0;i<3;i++) com[i] = com1[i] = 0;
    set<int>::iterator iter = s.result().begin(),
      end = s.result().end();
    double * mp = mass;
    while(iter!=end){
      int pos = (*iter)*Coor<double>::DIM;
      com[0] += (*mp)*c2.coor[pos];
      com[1] += (*mp)*c2.coor[pos+1];
      com[2] += (*mp)*c2.coor[pos+2];
      com1[0] += (*mp)*c1.coor[pos];
      com1[1] += (*mp)*c1.coor[pos+1];
      com1[2] += (*mp)*c1.coor[pos+2];
      ++mp;
      ++iter;
    }

    for(int i=0;i<3;i++) com[i] /= M;
    for(int i=0;i<3;i++) com1[i] /= M;

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;
      c2.coor[pos] -= com[0];
      c2.coor[pos+1] -= com[1];
      c2.coor[pos+2] -= com[2];
      c1.coor[pos] -= com1[0];
      c1.coor[pos+1] -= com1[1];
      c1.coor[pos+2] -= com1[2];    
    }

    double rmsd = mrmsd_mat(c1,c2,s.result(),mass_sqrt,r);

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;

      c1.coor[pos] = 
	r[0][0]*c2.coor[pos] +
	r[0][1]*c2.coor[pos+1] +
	r[0][2]*c2.coor[pos+2];

      c1.coor[pos+1] = 
	r[1][0]*c2.coor[pos] +
	r[1][1]*c2.coor[pos+1] +
	r[1][2]*c2.coor[pos+2];

      c1.coor[pos+2] = 
	r[2][0]*c2.coor[pos] +
	r[2][1]*c2.coor[pos+1] +
	r[2][2]*c2.coor[pos+2];
    }

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;
      c1.coor[pos] += com1[0];
      c1.coor[pos+1] += com1[1];
      c1.coor[pos+2] += com1[2];
    }

    for(int i=0;i<c1.size();i++){
      int pos = i*Coor<double>::DIM;
      //assign
      m2[i][0] = c1.coor[pos];
      m2[i][1] = c1.coor[pos+1];
      m2[i][2] = c1.coor[pos+2];
    }

    delete [] mass;
    delete [] mass_sqrt;
    return sqrt(size/M)*rmsd;
  }

  void rmsd_prof(const Molecules & m1, Molecules & m2,string file){
    ofstream out;
    out.open(file.c_str());
    for(int i=0;i<m1.fragment_size();i++){
      vector<int> a = m1.fragment_atoms(i);
      double d2 = 0;
      for(int j=0;j<a.size();j++){
	double d = m1[a[j]][0] - m2[a[j]][0];
	d2 += d*d;
	d = m1[a[j]][1] - m2[a[j]][1];
	d2 += d*d;
	d = m1[a[j]][2] - m2[a[j]][2];
	d2 += d*d;
      }
      out<<i+1<<" "<<sqrt(d2/a.size())<<" "<<a.size()<<" "<<m1.fragment_name(i)<<"\n";
    }
  }

  double subfit(const Molecules & m1, Molecules & m2, const NameSelection &s,
		const NameSelection &s2){
    if(fit(m1,m2,s)<0) return -1;
    Coor<double> c1(m1);
    Coor<double> c2(m2);
    set<int>::iterator iter = s2.result().begin(),
      end = s2.result().end();
    double rmsd = 0;
    while(iter!=end){
      int pos = *iter;
      pos = Coor<double>::DIM*pos;
      double dx = c1.coor[pos] - c2.coor[pos];
      double dy = c1.coor[pos+1] - c2.coor[pos+1];
      double dz = c1.coor[pos+2] - c2.coor[pos+2];
      rmsd += dx*dx+dy*dy+dz*dz;
      ++iter;
    }
    return sqrt(rmsd/s2.result().size());
  }


  void set_torsion(Protein & protein, int s1, int s2,
		   int s3, int s4, double tor){


    BondRotation<double> bond_rotation(protein);
    vector<double*> coor;
    for(int i= protein.atom_offset();
	i<protein.atom_offset()+ protein.atom_size();i++)
      coor.push_back(const_cast<double*>(protein[i].coor()));

    bond_rotation.set_dihedral(coor,s1,s2,s3,s4,tor);
  }

  void set_torsion2(Protein & protein, int s1, int s2,
		    int s3, int s4, double tor,
		    BondRotation<double> & bond_rotation ){
    vector<double*> coor;
    for(int i= protein.atom_offset();
	i<protein.atom_offset()+ protein.atom_size();i++)
      coor.push_back(const_cast<double*>(protein[i].coor()));

    bond_rotation.set_dihedral(coor,s1,s2,s3,s4,tor);
  }

  double ramadist(double phi1, double psi1, 
		  double phi2, double psi2){
    double d1 = abs(phi1-phi2);
    double t;
    t = abs(phi1-phi2+360);     
    if(t<d1) d1 = t;
    t = abs(phi1-phi2-360);
    if(t<d1) d1 = t;


    double d2 = abs(psi1-psi2);
    t = abs(psi1-psi2+360);     
    if(t<d2) d2 = t;
    t = abs(psi1-psi2-360);
    if(t<d2) d2 = t;
    return sqrt(d1*d1+d2*d2);
  }


  double ramarmsd(const Protein & p1,
		  const Protein & p2,
		  int from,
		  int to,
		  double cut){
    vector<double> phi1,psi1,phi2,psi2;
      phi1 = p1.phi();
      phi2 = p2.phi();

      psi1 = p1.psi();
      psi2 = p2.psi();

      if(phi1.size()!=phi2.size()) return -1;

      int size = phi1.size();
      double todeg = (180.0/M_PI);
      for(int i=0;i<size;i++){
	if(phi1[i]!=999) phi1[i] *= todeg;
	if(phi2[i]!=999) phi2[i] *= todeg;
	if(psi1[i]!=999) psi1[i] *= todeg;
	if(psi2[i]!=999) psi2[i] *= todeg;
      }

      from = from -1;

      from = max(from,0);
      to = min(to,(int)phi1.size());
      
      int c = 0;
      int g = 0;
      double drama = 0;
      for(int i=from;i<to;i++){
	if(phi1[i]!=999&&phi2[i]!=999&&psi1[i]!=999&&psi2[i]!=999){
	  double d = ramadist(phi1[i],psi1[i],
			      phi2[i],psi2[i]);
	  drama += 0.5*d*d;
	  aout<<">> "
	      <<i+1<<" "<<phi1[i]<<" "<<psi1[i]<<" "
	      <<phi2[i]<<" "<<psi2[i]<<" "<< d<<"\n";
	  if(d<cut) g++;
	  ++c;
	}
      }
      aout<<">> Good ramachandran amino-acids: "<<(100.0*g)/c<<endl;
      return sqrt(drama/c);
  }

  void set_beta_backbone(Protein & p){
    vector<double> phi = p.phi();
    vector<double> psi = p.psi();
    for(int i=0;i<phi.size();i++){
      cout<<"Rama "<<i+1<<" "<<phi[i]<<" "<<psi[i]<<" "<<"\n";
    }
    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);
    }
    
    double Phi = (-135.0)*(M_PI/180.0);
    double Psi = (135.0)*(M_PI/180.0);
    for(int i=0;i<p.fragment_size();i++){
      if(phi[i]!=999){
	set_torsion2(p,C[i-1],N[i],CA[i],C[i],Phi,br);
      }
      if(psi[i]!=999){
	set_torsion2(p,N[i],CA[i],C[i],N[i+1],Psi,br);
      }
    }    
  }

  void randomize_backbone(Protein & p, int seed){
    srand48(seed);
    vector<double> phi = p.phi();
    vector<double> psi = p.psi();
    for(int i=0;i<phi.size();i++){
      cout<<"Rama "<<i+1<<" "<<phi[i]<<" "<<psi[i]<<" "<<"\n";
    }
    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);
    }
    
    //    double Phi = (-135.0)*(M_PI/180.0);
    //    double Psi = (135.0)*(M_PI/180.0);
    for(int i=0;i<p.fragment_size();i++){
      if(phi[i]!=999){
	double Phi = 180*(2*drand48()-1)*(M_PI/180.0);
	set_torsion2(p,C[i-1],N[i],CA[i],C[i],Phi,br);
      }
      if(psi[i]!=999){
	double Psi = 180*(2*drand48()-1)*(M_PI/180.0);
	set_torsion2(p,N[i],CA[i],C[i],N[i+1],Psi,br);
      }
    }    
  }    

  double rgyr(const Molecules & mols, const NameSelection & s){
    Coor<double> c1(mols);
    return calc_rgyr(c1,s.result());
  }

  void render_tachyon(const Molecules & m, string file){
    ofstream dat;
    dat.open(file.c_str());
    vector<double> vdw;
    for(int i=0;i<m.atom_size();i++){
      string at = m[i].type_name();
      double d = m.protein(0).nbdata().sigma(at);
      vdw.push_back(d);

    }
    dat<<"Begin_Scene\n";
    dat<<"Resolution 600 600\n";
    dat<<"Camera\n";
    dat<<"  Projection Orthographic\n";
    dat<<"  Zoom 0.025\n";
    dat<<"  Aspectratio 1.00000\n";
    dat<<"  Antialiasing 12\n";
    dat<<"  Raydepth 30\n";
    dat<<"  Center  0.000000 0.000000 0.000000\n";
    dat<<"  Viewdir -0.000000 -0.000000 -2.000000\n";
    dat<<"  Updir   0.000000 1.000000 -0.000000\n";
    dat<<"End_Camera\n";
    dat<<"Directional_Light Direction 0.100000 -0.100000 1.000000 Color 1.000 1.000 1.000\n";
    dat<<"Directional_Light Direction -1.000000 -2.000000 0.500000 Color 1.000 1.000 1.000\n";
    dat<<"\n";
    dat<<"Background 1.000 1.000 1.000\n";

    for(int i=0;i<m.atom_size();i++){
      if(m[i].name()[0]=='H') continue;
      dat<<"Sphere\n";
      dat<<"  Center "<<m[i].x()<<" "<<m[i].y()<<" "<<m[i].z()-30<<"\n";
      dat<<"  Rad    "<<vdw[i]<<"\n";
      dat<<"      Texture\n";
      if(m[i].name()[0]=='C'){
	dat<<"  Ambient 0.000 Diffuse 0.700 Specular 0.000 Opacity 1.000\n";
	dat<<"  Phong Plastic 0.000 Phong_size 38.905 Color 0.500 0.490 1.000 TexFunc 0\n";
      }
      else if(m[i].name()[0]=='N'){
	dat<<"  Ambient 0.000 Diffuse 0.700 Specular 0.000 Opacity 1.000\n";
	dat<<"  Phong Plastic 0.000 Phong_size 38.905 Color 0.500  1.000 0.490 TexFunc 0\n";
      }
      else {
	dat<<"  Ambient 0.000 Diffuse 0.700 Specular 0.000 Opacity 1.000\n";
	dat<<"  Phong Plastic 0.000 Phong_size 38.905 Color 1.000 0.500  0.490 TexFunc 0\n";
      }
    }
    dat<<"End_Scene\n";
  }


  void find_clashes(const Molecules & m, MDB & mdb){
    C19ExclusionList excl(m,mdb);
    double q = pow(2.0,1.0/6.0);
    q = 1/q;
    for(int i=0;i<m.atom_size();i++){
      double s1 = mdb.nbondT(m[i].type_name()).sigma();
      for(int j=i+1;j<m.atom_size();j++){
	if(excl.excluded(i,j)) continue;
	double s2 = mdb.nbondT(m[j].type_name()).sigma();
	double d = m.distance(i,j);
	double sigma = s1+s2;
	double sigma2 = sigma*sigma;
	double sigma6 = sigma2*sigma2*sigma2;
	double sigma12 = sigma6*sigma6;
	double r6 = d*d*d*d*d*d;
	double r12 = r6*r6;
	double e = (sigma12/r12-2*sigma6/r6);
	if(d<q*(s1+s2)){
	  aout<<">> Steric clash "<<m.atom_name(i,Molecules::BASE)<<" "<<m.atom_name(j,Molecules::BASE)<<" "
	      <<s1<<" "<<s2<<" "<<d<<" "<<q*(s1+s2)<<" "<<e<<endl;
	}
      }
    }

  }

  double * nsc_sourface(const Molecules & mol, MDB & mdb, int &ndots){
    ndots = 0;
    double * coor;
    coor = new double[mol.atom_size()*3];
    for(int i=0;i<mol.atom_size();i++){
      coor[3*i] =  mol[i].x();
      coor[3*i+1] =  mol[i].y();
      coor[3*i+2] = mol[i].z();
    }

    double * rad = new double[mol.atom_size()];
    for(int i=0;i<mol.atom_size();i++){
      NBondT nbt = mdb.nbondT(mol[i].type_name());
      rad[i] = nbt.sigma()+1.4; //??
    }

    int natoms = mol.atom_size();
    int density = 300;
    int flag = FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS;
    double area = 0;
    double vol = 0;
    double *atomwise_area, *surface_dots; 
    nsc_dclm(coor,rad,natoms,density,flag,&area,&atomwise_area,&vol,
	     &surface_dots, &ndots);

//     for(int i=0;i<ndots;i++){
//       cout<<surface_dots[3*i]<<" "<<surface_dots[3*i+1]<<" "<<surface_dots[3*i+2]<<endl;
//     }

    //free(surface_dots);
    delete [] coor;
    delete [] rad;
    free(atomwise_area);
    return surface_dots;
  }

  int find_disu(const Protein & p, double max_dist){
    vector<int> cysd;
    for(int i=0;i<p.fragment_size();i++){
      int pos = p.fragment_name(i).find("_");
      string base;
      string frg = p.fragment_name(i);
      if(pos!=string::npos){
	base = string(frg.begin(),
		      frg.begin()+pos);
      }
      if(base=="CYS"){
	cysd.push_back(i);
      }
    }
    int count =0;
    for(int i=0;i<cysd.size();i++){
      for(int j=i+1;j<cysd.size();j++){
	int s1 = p.find_atom(p.fragment_name(cysd[i])+"/"+"SG");
	int s2 = p.find_atom(p.fragment_name(cysd[j])+"/"+"SG");
	double d = p.distance(s1,s2);
	if(d<=max_dist) ++count;
      }
    }
    return count;
  }
}

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