// Dokholyan's Component Library (DCL).
// Created by Nikolay Dokholyan.
// Copyright 1996-2000
 
//            Part xx:
// PDB manipulaion class


#ifndef _DCLPDB0_H_
#define _DCLPDB0_H_
 
/*
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
*/

#include "header_cpp.h"


#ifndef ABS
#define ABS(x)   ( ((x)<0) ? (-(x)) : (x))
#endif

class vec{
 public:
  vec(){};
  vec(double a, double b, double c){x=a;y=b;z=c;};
  ~vec(){};
  double x,y,z;
  double norm(void){return sqrt(x*x + y*y +z*z);};
  double norm2(void){return x*x + y*y +z*z;};
};



class DCLPDB{ 

  char pfile[100];   // PDB file
  ifstream in;       // stream
  
  
  typedef struct{
    vec r;           // position
//    char *name;
    char name[6];    // actual name
    int CA;          // CA flag, default = 0 (not CA)
    int CB;          // CB flag, default = 0 (not CB)
    int heavy;       // heavy atom (0 if H, 1 otherwise)
    double T;        // temperature factor
  } atom;
  
  class residue{
  public:
    residue(){};
    ~residue(){delete [] a;};
    char name[100];  // name of an amino acid
    int na;          // number of atoms
    atom a[40];         // pointer to an array of atoms belonging to aa
    int pCA;         // position of CA within the aa
    int pCB;         // position of CB within the aa
  };

  
  
 public:
  int N;             // Number of residues
  int Na;            // Number of atoms
  int ipos;          // initial number of the first residue
  residue res[1000];      // residues
  int CM[1000][1000];          // Contact matrix
  vec COM;           // Center of Mass vector
  char **Seq3;       // Sequence in a 3-letter representation
  char *Seq;         // Sequence in a 1-letter representation



  DCLPDB();
  DCLPDB(char *);
  ~DCLPDB();


  void Read(void);                  // reads the PDB file 
  void PrintCheckPdb();


  void GetCAContMap(double);        // returns contact map based on CAs
  void GetCBContMap(double);        // returns contact map based on CBs
  void GetCenterContMap(double);    // returns contact map based on sidechain centers

  double D(int,int,int,int);  // returns the distance between two atoms
  double D2(int,int,int,int); // returns the distance square between two atoms
  double DCA(int,int);    // returns the distance between two CAs
  double DCA2(int,int);   // returns the distance square between two CAs
  double DCB(int,int);    // returns the distance between two CBs
  double DCB2(int,int);   // returns the distance square between two CBs
  double Dist2(vec *);        // returns the square of vector size
  double Dist2(vec *, vec *); // returns the square of diff between 2 vectors

  void GetSeq(void);          // initializes a seq string
  vec *GetSidechainCenter(residue &r);
  double GetDistance(vec *, vec*);

 private:
  void squeeze(char [],int);
  void clean(char s[], int n){for(int i=n;i<strlen(s);i++)s[i]='\0';};
  void clear(char s[]){for(int i=0;i<strlen(s);i++)s[i]='\0';};
  void clear(char s[], int n){for(int i=0;i<n;i++)s[i]='\0';};
  void Error(char *msg){cerr<<msg<<endl<<flush; exit(1);};
  char tmp[100];
  char aa1[20], aa3[20][3];
};


DCLPDB::DCLPDB(){
  N=Na=0;
  aa1[0]='A';  aa1[1]='C';  aa1[2]='D';  aa1[3]='E';
  aa1[4]='F';  aa1[5]='G';  aa1[6]='H';  aa1[7]='I';
  aa1[8]='K';  aa1[9]='L';  aa1[10]='M'; aa1[11]='N';
  aa1[12]='P'; aa1[13]='Q'; aa1[14]='R'; aa1[15]='S';
  aa1[16]='T'; aa1[17]='V'; aa1[18]='W'; aa1[19]='Y';

  strcpy(aa3[0],"ALA");  strcpy(aa3[1],"CYS");
  strcpy(aa3[2],"ASP");  strcpy(aa3[3],"GLU");
  strcpy(aa3[4],"PHE");  strcpy(aa3[5],"GLY");
  strcpy(aa3[6],"HIS");  strcpy(aa3[7],"ILE");
  strcpy(aa3[8],"LYS");  strcpy(aa3[9],"LEU");
  strcpy(aa3[10],"MET"); strcpy(aa3[11],"ASN");
  strcpy(aa3[12],"PRO"); strcpy(aa3[13],"GLN");
  strcpy(aa3[14],"ARG"); strcpy(aa3[15],"SER");
  strcpy(aa3[16],"THR"); strcpy(aa3[17],"VAL");
  strcpy(aa3[18],"TRP"); strcpy(aa3[19],"TYR");

};


DCLPDB::DCLPDB(char *afname){
  strcpy(pfile,afname);
  N=Na=0;

  aa1[0]='A';  aa1[1]='C';  aa1[2]='D';  aa1[3]='E';
  aa1[4]='F';  aa1[5]='G';  aa1[6]='H';  aa1[7]='I';
  aa1[8]='K';  aa1[9]='L';  aa1[10]='M'; aa1[11]='N';
  aa1[12]='P'; aa1[13]='Q'; aa1[14]='R'; aa1[15]='S';
  aa1[16]='T'; aa1[17]='V'; aa1[18]='W'; aa1[19]='Y';

  strcpy(aa3[0],"ALA");  strcpy(aa3[1],"CYS");
  strcpy(aa3[2],"ASP");  strcpy(aa3[3],"GLU");
  strcpy(aa3[4],"PHE");  strcpy(aa3[5],"GLY");
  strcpy(aa3[6],"HIS");  strcpy(aa3[7],"ILE");
  strcpy(aa3[8],"LYS");  strcpy(aa3[9],"LEU");
  strcpy(aa3[10],"MET"); strcpy(aa3[11],"ASN");
  strcpy(aa3[12],"PRO"); strcpy(aa3[13],"GLN");
  strcpy(aa3[14],"ARG"); strcpy(aa3[15],"SER");
  strcpy(aa3[16],"THR"); strcpy(aa3[17],"VAL");
  strcpy(aa3[18],"TRP"); strcpy(aa3[19],"TYR");

};





DCLPDB::~DCLPDB(){
};




//Reads the pdb file 
void DCLPDB::Read(void){
  char line[100];
  N=Na=0;

  // determining number of residues
  in.open(pfile,ios::in);
  if(!in) {
     cout << "pdb file cannot be read" << endl;
     exit(0);
  }
  int atomcount = 0;
  int residuecount = -1;
  char residuenum[10], oldresiduenum[10];
  strcpy(residuenum, "1"); strcpy(oldresiduenum, "0");
  while(!in.eof()){
    clear(line);
    in.getline(line,100,'\n');
    if(!strncmp(line,"ENDMDL",6))break; // select only one model
    if(!strncmp(line,"TER",3))break;    // select only one chain
    if(!strncmp(line,"END",3))break;    // select only one chain
    if((!strncmp(line,"ATOM",4))){
	strncpy(residuenum, line+22, 6);
	residuenum[6] = '\0';
	if(strcmp(residuenum, oldresiduenum)) {
	   res[residuecount].na = atomcount;
	   atomcount = 0;
	   residuecount++;
	   res[residuecount].pCB = -1;
	   //cout << residuecount<< endl;
	   strncpy(res[residuecount].name,line+17, 3);
           res[residuecount].name[3] = '\0';
	   strcpy(oldresiduenum, residuenum);
	}
	strncpy(res[residuecount].a[atomcount].name,line+13, 3);
	res[residuecount].a[atomcount].name[3] = '\0';
	if(!strncmp(res[residuecount].a[atomcount].name,"CA",2) ) {
	    res[residuecount].a[atomcount].CA = 1; 
	    res[residuecount].pCA = atomcount;
	}
	else res[residuecount].a[atomcount].CA = 0;
	if(!strncmp(res[residuecount].a[atomcount].name,"CB",2) ) {
	    res[residuecount].a[atomcount].CB = 1; 
	    res[residuecount].pCB = atomcount;
	}
	else res[residuecount].a[atomcount].CB  = 0;

	// reading coordinates
	sscanf(line+30, "%lf", &res[residuecount].a[atomcount].r.x);
	sscanf(line+38, "%lf", &res[residuecount].a[atomcount].r.y);
	sscanf(line+46, "%lf", &res[residuecount].a[atomcount].r.z);

	atomcount ++ ;
    }
  }
  res[residuecount].na = atomcount;
  N = residuecount+1;
 
  in.close();
};


void DCLPDB::PrintCheckPdb() {
  int i,j;
  
  cout << "There are " << N << " residues"<< endl;
  for(i=0;i<N;i++) {
    for(j=0;j<res[i].na;j++) {
	cout << i << "\t" << res[i].a[j].name <<" " << res[i].name << " " << res[i].a[j].r.x << "\t" << res[i].a[j].r.y << "\t" << res[i].a[j].r.z << endl;
    }
  }

} 




// GetCAContMap(double) return the contact matrix 
// at threshold distance thr with all CAs of the protein
void DCLPDB::GetCAContMap(double thr){
  int i,j;
  //CM = new int*[N];
  //for(i=0;i<N;i++) CM[i] = new int[N];
  if(!CM)Error("Cannot allocate memory for CM");
  for(i=0;i<N;i++)
    for(j=0;j<N;j++)CM[i][j]=0;

  thr*=thr;

  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      if(i==j)j++;
      if(j==N)break;
      if(Dist2(&res[j].a[res[j].pCA].r,&res[i].a[res[i].pCA].r)<=thr){
		CM[i][j]=CM[j][i]=1;
      }
    }
  }
};



// GetCBContMap(double) return the contact matrix 
// at threshold distance thr with all CBs of the protein
// if CB does not exist (GLY), use CA 
void DCLPDB::GetCBContMap(double thr){
  int i,j;
  int itmp=0,jtmp=0;
  //CM = new int*[N];
  //for(i=0;i<N;i++) CM[i] = new int[N];
  if(!CM)Error("Cannot allocate memory for CM");
  for(i=0;i<N;i++)
    for(j=0;j<N;j++)CM[i][j]=0;

  thr*=thr;

  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      if(i==j)j++;
      if(j==N)break;

      if(res[i].pCB==-1){res[i].pCB=res[i].pCA;itmp=1;}
      if(res[j].pCB==-1){res[j].pCB=res[j].pCA;jtmp=1;}
	
      //cout << "i" << i << "j" << j<< endl;
      //cout << res[i].pCB << "\t" << res[j].pCB << endl;
      //cout << res[j].a[res[j].pCB].r.x << "\t" <<  res[j].a[res[j].pCB].r.y<< "\t"<<  res[j].a[res[j].pCB].r.z << endl;
       //cout << res[i].a[res[i].pCB].r.x << "\t" <<  res[i].a[res[i].pCB].r.y<< "\t"<<  res[i].a[res[i].pCB].r.z << endl;
      if(Dist2(&res[j].a[res[j].pCB].r,&res[i].a[res[i].pCB].r)<=thr)
		CM[i][j]=CM[j][i]=1;
      if(itmp){res[i].pCB=-1;itmp=0;}
      if(jtmp){res[j].pCB=-1;jtmp=0;}
    }
  }
};

void DCLPDB::GetCenterContMap(double thr){
  int i,j;
  int itmp=0,jtmp=0;
  //CM = new int * [N];
  //for(i=0;i<N;i++) {CM[i] = new int [N]; for(j=0;j<N;j++) CM[i][j] = 0;cout << CM[i][75]<< endl;;if(!CM[i]) cout << i << endl;}
  if(!CM)Error("Cannot allocate memory for CM");
  for(i=0;i<N;i++)
    for(j=0;j<N;j++) {CM[i][j]=0;}

  double dist;
  
  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      if(i==j) {CM[i][j] = 0; continue;}
   
      vec *v1, *v2;
      v1 = GetSidechainCenter(res[i]);
      v2 = GetSidechainCenter(res[j]);
      dist = GetDistance(v1, v2);
      if(dist>thr) CM[i][j] = 0;
      else CM[i][j] = 1;
      delete v1; delete v2;
    }
  }  
};


// initializes a seq string
void DCLPDB::GetSeq(void){
  int i,j;
  Seq3 =  new char*[N];
  for(i=0;i<N;i++)Seq3[i] = new char[3];

  Seq = new char [N];

  //for(i=0;i<20;i++) cout << aa1[i] << endl;
  for(i=0;i<N;i++){
    strcpy(Seq3[i],res[i].name);
    //cout << Seq3[i] << " " << res[i].name << endl;
    for(j=0;j<20;j++){
      if(!strncmp(res[i].name,aa3[j],3)){
	Seq[i]=aa1[j];
	//cout << j<< " " << aa1[j] <<" " <<  Seq[i] << endl;
	break;
      }
    }
  }
};


// returns the distance between atoms j and m of amino acids i and k 
double DCLPDB::D(int i, int j, int k, int m){
  if(j<res[i].na && m<res[k].na){
    return sqrt((res[i].a[j].r.x-res[k].a[m].r.x)*\
		(res[i].a[j].r.x-res[k].a[m].r.x)+\
		(res[i].a[j].r.y-res[k].a[m].r.y)*\
		(res[i].a[j].r.y-res[k].a[m].r.y)+\
		(res[i].a[j].r.z-res[k].a[m].r.z)*\
		(res[i].a[j].r.z-res[k].a[m].r.z));
  }
  else{cerr<<"Atoms are out of bounds\n"<<flush;return -1;}
};


// returns the square of distance between atoms j and m of amino acids i and k 
double DCLPDB::D2(int i, int j, int k, int m){
  if(j<res[i].na && m<res[k].na){
    return ((res[i].a[j].r.x-res[k].a[m].r.x)*\
	    (res[i].a[j].r.x-res[k].a[m].r.x)+\
	    (res[i].a[j].r.y-res[k].a[m].r.y)*\
	    (res[i].a[j].r.y-res[k].a[m].r.y)+\
	    (res[i].a[j].r.z-res[k].a[m].r.z)*\
	    (res[i].a[j].r.z-res[k].a[m].r.z));
  }
  else{cerr<<"Atoms are out of bounds\n"<<flush;return -1;}
};


// returns the distance between two CAs
double DCLPDB::DCA(int i, int k){
  int j=res[i].pCA;
  int m=res[k].pCA;
  return sqrt((res[i].a[j].r.x-res[k].a[m].r.x)*\
	      (res[i].a[j].r.x-res[k].a[m].r.x)+\
	      (res[i].a[j].r.y-res[k].a[m].r.y)*\
	      (res[i].a[j].r.y-res[k].a[m].r.y)+\
	      (res[i].a[j].r.z-res[k].a[m].r.z)*\
	      (res[i].a[j].r.z-res[k].a[m].r.z));
};


// returns the square distance between two CAs
double DCLPDB::DCA2(int i, int k){
  int j=res[i].pCA;
  int m=res[k].pCA;
  return ((res[i].a[j].r.x-res[k].a[m].r.x)*\
	  (res[i].a[j].r.x-res[k].a[m].r.x)+\
	  (res[i].a[j].r.y-res[k].a[m].r.y)*\
	  (res[i].a[j].r.y-res[k].a[m].r.y)+\
	  (res[i].a[j].r.z-res[k].a[m].r.z)*\
	  (res[i].a[j].r.z-res[k].a[m].r.z));
};


// returns the distance between two CBs
double DCLPDB::DCB(int i, int k){
  int j=res[i].pCB;
  if(j==-1)j=res[i].pCA;
  int m=res[k].pCB;
  if(m==-1)m=res[k].pCA;
  return sqrt((res[i].a[j].r.x-res[k].a[m].r.x)*\
	      (res[i].a[j].r.x-res[k].a[m].r.x)+\
	      (res[i].a[j].r.y-res[k].a[m].r.y)*\
	      (res[i].a[j].r.y-res[k].a[m].r.y)+\
	      (res[i].a[j].r.z-res[k].a[m].r.z)*\
	      (res[i].a[j].r.z-res[k].a[m].r.z));
};


// returns the square distance between two CBs
double DCLPDB::DCB2(int i, int k){
  int j=res[i].pCB;
  if(j==-1)j=res[i].pCA;
  int m=res[k].pCB;
  if(m==-1)m=res[k].pCA;
  return ((res[i].a[j].r.x-res[k].a[m].r.x)*\
	  (res[i].a[j].r.x-res[k].a[m].r.x)+\
	  (res[i].a[j].r.y-res[k].a[m].r.y)*\
	  (res[i].a[j].r.y-res[k].a[m].r.y)+\
	  (res[i].a[j].r.z-res[k].a[m].r.z)*\
	  (res[i].a[j].r.z-res[k].a[m].r.z));
};


vec * DCLPDB::GetSidechainCenter(residue &res) {
     int Na = res.na;
     int i;
     double posx=0, posy = 0, posz =0;
     int count=0;
     for(i=0;i<Na;i++ ) {
          if( (strcmp(res.a[i].name, "N")==0) || (strcmp(res.a[i].name, "C")==0) || (strcmp(res.a[i].name, "O") ==0 ) || (strcmp(res.a[i].name, "CA")==0 ) ) continue;
          if(res.a[i].heavy) {
              count ++ ;
              posx += res.a[i].r.x;
              posy += res.a[i].r.y;
              posz += res.a[i].r.z;
          }
     }
     if(count>0) {
      posx = posx/count;
      posy = posy/count;
      posz = posz/count;
     }
     else {
      posx = res.a[res.pCA].r.x;
      posy = res.a[res.pCA].r.y;
      posz = res.a[res.pCA].r.z;     
     }

     vec *v = new vec(posx, posy, posz);
     return v;
};

double DCLPDB::GetDistance( vec *v1, vec *v2) {

    double value;

    value = (v1->x-v2->x)* (v1->x-v2->x);
    value +=  (v1->y - v2->y)*(v1->y-v2->y);
    value += (v1->z-v2->z)*(v1->z-v2->z);
    return (sqrt(value));
};

void DCLPDB::squeeze(char s[], int c){
  int i, j;
  for(i = j= 0; s[i] != '\0'; i++)
    if (s[i] != c)
      s[j++] = s[i];
  s[j] = '\0';
};



// returns the square of vector size
double DCLPDB::Dist2(vec *a){
  return (a->x*a->x+a->y*a->y+a->z*a->z);
};

// returns the square of diff between 2 vectors
double DCLPDB::Dist2(vec *a, vec *b){
  return ((a->x - b->x)*(a->x - b->x) + \
	  (a->y - b->y)*(a->y - b->y) + \
	  (a->z - b->z)*(a->z - b->z));
};

#endif /* _DCLPDB0_H_ */










