/***************************************************************************
                          mkbetapair.cpp  -  description
                             -------------------
    begin                : Tue Nov 18 01:45:35 2003
    copyright            : (C) 1999-2006 by Cavalli Andrea
    author               : : cavalli $
    date                 : : 2003/12/09 11:10:33 $
    id                   : : coor.h,v 1.1.2.2 2003/12/09 11:10:33 cavalli Exp $
    email                : cavalli@bioc.unizh.ch amc82@cam.ac.uk
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <almost.h>
#include <algorithm/molalgo.h>
#include <algorithm/mkbetapair.h>
namespace Almost {
  bool MakeBetaPair::apply(Protein & p1,
			   Protein & p2){
    //Assume build from scratch else
    // fit ...
    if(parallel_)
      makeparallel(p1,p2);
    else makeantiparallel(p1,p2);

    return true;
  }

  bool MakeBetaPair::build(const vector<string> & s1,
			   const vector<string> & s2,
			   Molecules & mols, MDB & mdb){

    //build prot 1
    Protein p1("TRG");
    Protein p2("IHB");
    p1.build_sequ(s1,mdb,"DEFAULT","DEFAULT");
    p2.build_sequ(s2,mdb,"DEFAULT","DEFAULT");

    apply(p1,
	  p2);
    mols.add_protein(p1);
    mols.add_protein(p2);

    return true;
  }


  void MakeBetaPair::makeparallel(Protein & p1,
				  Protein & p2){
    vector<string> sequ1 = p1.sequence(Protein::LONG);
    vector<string> sequ2 = p2.sequence(Protein::LONG);
    if(up_){
      swap(p1,sequ1);
      swap(p2,sequ2);
    }
    align_y(p1,sequ1);

    if(betabb_){
      set_beta_backbone(p1);
    }

    //Applay offset
    if(offset_>0){
      for(int q=0;q<offset_;q++){
	Coor3D c,c1;
	int a1 = p2.find_atom(sequ1[0]+"/CA");
	int a2 = p2.find_atom(sequ1[1]+"/CA");
	int a3 = p2.find_atom(sequ1[2]+"/CA");
	c[0]=p2[a2][0] - p2[a1][0];
	c[1]=p2[a2][1] - p2[a1][1];
	c[2]=p2[a2][2] - p2[a1][2]; 
	p2.translate(c);
	c[0]=p2[a3][0] - p2[a1][0];
	c[1]=p2[a3][1] - p2[a1][1];
	c[2]=p2[a3][2] - p2[a1][2]; 
	double n= 0;
	for(int k=0;k<3;k++) n += c[k]*c[k];
	n =sqrt(n);
	for(int k=0;k<3;k++) c[k] = c[k]/n;
	c1[0] = p2[a1][0];
	c1[1] = p2[a1][1];
	c1[2] = p2[a1][2];
	p2.rotate(M_PI,c,c1);
      }
    }

    Coor3D c;
    int a1 = p2.find_atom(sequ2[0]+"/CA");
    int a2 = p2.find_atom(sequ2[1]+"/CA");
    c[0]=p2[a2][1] - p2[a1][1];
    c[1]=-(p2[a2][0] - p2[a1][0]);
    c[2]=p2[a2][2] - p2[a1][2]; 
	
    //scale c
    double n= 0;
    for(int k=0;k<3;k++) n += c[k]*c[k];
    n =sqrt(n);
    for(int k=0;k<3;k++) c[k] = 5.0*c[k]/n;
    p2.translate(c);
    align_y(p2,sequ2);
    if(betabb_)
      set_beta_backbone(p2);
  }

  void MakeBetaPair::makeantiparallel(Protein & p1,
				      Protein & p2){
    vector<string> sequ1 = p1.sequence(Protein::LONG);
    vector<string> sequ2 = p2.sequence(Protein::LONG);
    if(up_){
      swap(p1,sequ1);
      //    swap(p2,sequ2);
    }
    align_y(p1,sequ1);
  
    if(betabb_){
      set_beta_backbone(p1);
    }
    //Do anti ...

    Coor3D c,c1,ct;
    int a1 = p2.find_atom(sequ2[0]+"/CA");
    int a2 = p2.find_atom(sequ2[sequ2.size()-1]+"/CA");

    c[0]=p2[a2][1] - p2[a1][1];
    c[1]=-(p2[a2][0] - p2[a1][0]);
    c[2]=p2[a2][2] - p2[a1][2]; 

    //scale c
    double n= 0;
    for(int k=0;k<3;k++) n += c[k]*c[k];
    n =sqrt(n);
    for(int k=0;k<3;k++) c[k] = c[k]/n;


    c1[0] = p2[a1][0];
    c1[1] = p2[a1][1];
    c1[2] = p2[a1][2];

    ct[0]=p2[a2][0]-p2[a1][0]+5*c[0];
    ct[1]=p2[a2][1]-p2[a1][1]+5*c[1];
    ct[2]=p2[a2][2]-p2[a1][2]+5*c[2];


    //rotate around c
    p2.rotate(M_PI,c,c1);


    if(offset_>0){
      for(int q=0;q<offset_;q++){
	Coor3D c,c1;
	int a1 = p2.find_atom(sequ1[0]+"/CA");
	//	int a2 = p2.find_atom(sequ1[1]+"/CA");
	int a3 = p2.find_atom(sequ1[2]+"/CA");
	c[0]=0.5*(p2[a1][0] - p2[a3][0]);
	c[1]=0.5*(p2[a1][1] - p2[a3][1]);
	c[2]=0.5*(p2[a1][2] - p2[a3][2]); 
	p2.translate(c);
	c[0]=p2[a3][0] - p2[a1][0];
	c[1]=p2[a3][1] - p2[a1][1];
	c[2]=p2[a3][2] - p2[a1][2]; 
	double n= 0;
	for(int k=0;k<3;k++) n += c[k]*c[k];
	n =sqrt(n);
	for(int k=0;k<3;k++) c[k] = c[k]/n;
	c1[0] = p2[a1][0];
	c1[1] = p2[a1][1];
	c1[2] = p2[a1][2];
	// p2.rotate(M_PI,c,c1);
      }
    }

    int diff = sequ1.size()-sequ2.size()+offset_;
    diff = abs(diff);
    if(up_&&((offset_%2)==0)){
      swap(p2,sequ2);
    }
    if(sequ2.size()%2) swap(p2,sequ2);
    p2.translate(ct);
    align_y(p2,sequ2);
  }


  void MakeBetaPair::swap(Protein & p, const vector<string> sequ){
    Coor3D c,c1,ct;
    int a1 = p.find_atom(sequ[0]+"/CA");
    int a2 = p.find_atom(sequ[sequ.size()-1]+"/CA");
    c[0]=p[a2][0]-p[a1][0];
    c[1]=p[a2][1]-p[a1][1];
    c[2]=p[a2][2]-p[a1][2];
    c1[0] = p[a1][0];
    c1[1] = p[a1][1];
    c1[2] = p[a1][2];
    double n= 0;
    for(int k=0;k<3;k++) n += c[k]*c[k];
    n =sqrt(n);
    for(int k=0;k<3;k++) c[k] = c[k]/n;
    p.rotate(M_PI,c,c1);
  }

  void MakeBetaPair::align_y(Protein & p, const vector<string> sequ){
    Coor3D c,c1,ct;
    int a1 = p.find_atom(sequ[0]+"/CA");
    int a2 = p.find_atom(sequ[sequ.size()-1]+"/CA");
    c[0]=p[a2][0]-p[a1][0];
    c[1]=p[a2][1]-p[a1][1];
    c[2]=p[a2][2]-p[a1][2];
    c1[0] = p[a1][0];
    c1[1] = p[a1][1];
    c1[2] = p[a1][2];
    double n= 0;
    for(int k=0;k<3;k++) n += c[k]*c[k];
    n =sqrt(n);
    for(int k=0;k<3;k++) c[k] = c[k]/n;
    p.rotate(30.0*(M_PI/180.0),c,c1);
  }
}
  // Local Variables:
  // mode:C++
  // End: