/* EGAD: dihedral_cartesian.cpp
   
Navin Pokala and Tracy Handel
Dept. of Molecular and Cell Biology
University of California, Berkeley
Copyright (C) 2003 Regents of the University of California
GNU Public License
Aug 12 2003

Absolutely no warranties are made or are implied with the use of this program or its parts.

This file contains functions for converting dihedral representations of molecules into coordinates, and vice versa.
   
*/

#include "dihedral_cartesian.h"


/* based on the Makeatom function in ROC (Fortran 77) by John R. Desjarlais */

/* returns the coordinates of atom P  
   P---contactAtom             
             \             
           twobondAtom------dihedAtom   
  bondlength = P-b bondlength  
  bondangle = orbital hybridization angle of P-b-c  
  dihedral = P-d dihedral angle, w/ respect to d, clockwise  
*/

CARTESIAN chi2xyz(CARTESIAN *contactAtom, CARTESIAN *twobondAtom, 
		  CARTESIAN *dihedAtom, double *bondlength, 
		  double *bondangle, double *dihedral)
{
  CARTESIAN r1,r2,r3,r4;
  CARTESIAN P;
  double one80_minus_bondangle;
  double cos_180_minus_bondangle;
  double sin_180_minus_bondangle;
  double cos_dihedral;
  double sin_dihedral;
  
  r1 = unitvector((*contactAtom),(*twobondAtom));
  r2 = unitvector((*dihedAtom),(*twobondAtom));
  r3 = unitcrossprod(r1,r2);
  r4 = unitcrossprod(r3,r1);
  
  one80_minus_bondangle = 180.0 - (*bondangle);
  cos_180_minus_bondangle = cosd(one80_minus_bondangle);
  sin_180_minus_bondangle = sind(one80_minus_bondangle);
  cos_dihedral = cosd(*dihedral);
  sin_dihedral = sind(*dihedral);

  P.x = contactAtom->x + 
    (*bondlength)*((r1.x)*cos_180_minus_bondangle + ((r3.x)*sin_dihedral+ (r4.x)*cos_dihedral)*sin_180_minus_bondangle);
  P.y = contactAtom->y + 
    (*bondlength)*((r1.y)*cos_180_minus_bondangle + ((r3.y)*sin_dihedral+ (r4.y)*cos_dihedral)*sin_180_minus_bondangle);
  P.z = contactAtom->z + 
    (*bondlength)*((r1.z)*cos_180_minus_bondangle + ((r3.z)*sin_dihedral+ (r4.z)*cos_dihedral)*sin_180_minus_bondangle);

  return P;

}

/* This function returns the dihedral angle between atoms A and D */
/* This is based on the "dihedral" subroutine from "phipsi" by John R. Desjarlais */

double dihedral(CARTESIAN *A, CARTESIAN *B, CARTESIAN *C, CARTESIAN *D)
{
  double dihed;
  double q1, q2, q3;
  double sign;
  CARTESIAN r1, r2, r3, r4, W;

  r1 = unitvector((*C),(*B));
  r2 = unitvector((*A),(*B));
  r3 = unitcrossprod(r2,r1);
  r4 = unitcrossprod(r1,r3);

  W.x = D->x - C->x;
  W.y = D->y - C->y;
  W.z = D->z - C->z;

  q1 = dotproduct(W,r1);
  
  W.x = W.x - q1*r1.x;
  W.y = W.y - q1*r1.y;
  W.z = W.z - q1*r1.z;

  q2 = dotproduct(W,r3);
  q3 = dotproduct(W,r4);

  if(q2>0)
    sign = -1;
  else
    sign = 1;

  dihed = sign*arccos(q3/(sqrt(q2*q2 + q3*q3)));

  return dihed;
}

/* 
   This function generates coordinates for a BACKBONE array that has been loaded w/ dihedrals
   bkbn[1] must contain backbone coords for the first residue 
*/

void build_BACKBONE(BACKBONE *bkbn)
{
  int i, chain_number;
  double len_CB_CA, angle_CB_CA_N, len_HA_CA, angle_HA_CA_N, len_C_CA, angle_C_CA_N, len_O_C, angle_O_C_CA,
    len_N_C, angle_N_C_CA, len_H_N, angle_H_N_C, len_CA_N, angle_CA_N_C;
  double dihed_HA, dihed_C, dihed_N, dihed_CA, dihed_CB, dihed_O, dihed_H;
  div_t chain_id;

  len_CB_CA=1.525; angle_CB_CA_N=111.10; 
  len_HA_CA=1.090; angle_HA_CA_N=109.50; 
  len_C_CA=1.522;  angle_C_CA_N=111.10;
  len_O_C=1.229;   angle_O_C_CA=120.50; 
  len_N_C=1.335;   angle_N_C_CA=116.60; 
  len_H_N=1.010;   angle_H_N_C=119.80;
  len_CA_N=1.449;  angle_CA_N_C=121.90;
    
  i=1;
  while(bkbn[i].seq_position!=ENDFLAG)
  {
	chain_id = div(bkbn[i].seq_position, 1000);
     chain_number = chain_id.quot;
  i=i+1;
  while(chain_id.quot == chain_number && bkbn[i].seq_position!=ENDFLAG)
    {
      dihed_N = bkbn[i-1].psi;
      bkbn[i].N = chi2xyz(&bkbn[i-1].C, &bkbn[i-1].CA, &bkbn[i-1].N, &len_N_C, &angle_N_C_CA, &dihed_N);
      dihed_CA = bkbn[i-1].omega;
      bkbn[i].CA = chi2xyz(&bkbn[i].N, &bkbn[i-1].C, &bkbn[i-1].CA, &len_CA_N, &angle_CA_N_C, &dihed_CA);
      dihed_C =  bkbn[i].phi;
      bkbn[i].C = chi2xyz(&bkbn[i].CA, &bkbn[i].N, &bkbn[i-1].C, &len_C_CA, &angle_C_CA_N, &dihed_C);
      dihed_O = bkbn[i].psi + 180;
      bkbn[i].O = chi2xyz(&bkbn[i].C, &bkbn[i].CA, &bkbn[i].N, &len_O_C, &angle_O_C_CA, &dihed_O);
      
      if(bkbn[i].O_CTE != NULL)
      {
	  dihed_O = bkbn[i].psi + 180;
	  bkbn[i].O = chi2xyz(&bkbn[i].C, &bkbn[i].CA, &bkbn[i].N, &len_O_C, &angle_O_C_CA, &dihed_O);
	  
	  dihed_O = bkbn[i].psi;
	  *bkbn[i].O_CTE = chi2xyz(&bkbn[i].C, &bkbn[i].CA, &bkbn[i].N, &len_O_C, &angle_O_C_CA, &dihed_O);
      }
      
      dihed_CB = bkbn[i].phi - 120;
      bkbn[i].CB = chi2xyz(&bkbn[i].CA, &bkbn[i].N, &bkbn[i-1].C, &len_CB_CA, &angle_CB_CA_N, &dihed_CB);
      dihed_HA = bkbn[i].phi + 120;
      bkbn[i].HA = chi2xyz(&bkbn[i].CA, &bkbn[i].N, &bkbn[i-1].C, &len_HA_CA, &angle_HA_CA_N, &dihed_HA);
      dihed_H  = bkbn[i-1].omega + 180;
      bkbn[i].H = chi2xyz(&bkbn[i].N, &bkbn[i-1].C, &bkbn[i-1].CA, &len_H_N, &angle_H_N_C, &dihed_H);

      if(bkbn[i].H_NTE != NULL)
      {
	dihed_H  = bkbn[i-1].omega + 60;
        bkbn[i].H = chi2xyz(&bkbn[i].N, &bkbn[i-1].C, &bkbn[i-1].CA, &len_H_N, &angle_H_N_C, &dihed_H);
	dihed_H  = bkbn[i-1].omega + 180;
	bkbn[i].H_NTE[1] = chi2xyz(&bkbn[i].N, &bkbn[i-1].C, &bkbn[i-1].CA, &len_H_N, &angle_H_N_C, &dihed_H);
	dihed_H  = bkbn[i-1].omega + 300;
	bkbn[i].H_NTE[2] = chi2xyz(&bkbn[i].N, &bkbn[i-1].C, &bkbn[i-1].CA, &len_H_N, &angle_H_N_C, &dihed_H);

      }

      ++i;
	chain_id = div(bkbn[i].seq_position, 1000);
    }
    }
}
    
/* given the position, the resparam, and the rotamers, this  
  function will generate the cartesian coords of the sidechain atoms into side */

void build_a_SideChain(SIDECHAIN *side, BACKBONE *bkbn, RESPARAM *resparam, int *res_position, double chi[])
{
  int sideAtomNumber;
  int endSideAtomNumber;
  double dihedral, dihed2;
  int i,n;
  static int *flag = NULL;
  static CARTESIAN *dum;

  if(flag==NULL)
  {
  flag = (int *)calloc(MAX_RES_SIZE,sizeof(int));
  dum = (CARTESIAN *)calloc(MAX_RES_SIZE,sizeof(CARTESIAN));
  }
   
  for(i=0;i<MAX_RES_SIZE;++i)
	flag[i]=0;

  side->seq_position = *res_position;
  side->resparam_ptr = resparam;

  if(resparam->ligand_flag==0)
  {
  /* the AMBER file contains info on N,H,CA,HA,CB,C,O which we consider as backbone */
 /* numSideAtoms = resparam->numberofAtoms-7; */
 
  /* AMBER file has atomnumber 4->8 as N->CB and the last two C and O */
  /* the first atom we consider sidechain starts at atomnumber 9 */
  /* endSideAtomNumber = numSideAtoms+8; */
	endSideAtomNumber = resparam->numberofAtoms+1;

  /* assign coords from BACKBONE to the foundation atoms */
  dum[8] = bkbn->CB; dum[6] = bkbn->CA; dum[4] = bkbn->N;
 
  n=0;
  for(sideAtomNumber=9;sideAtomNumber<=endSideAtomNumber;++sideAtomNumber)
    {
      if(resparam->atom[sideAtomNumber].dihedral<0)   /* variable chi angle */
	{
	  ++n; flag[sideAtomNumber]=0;
	  dihedral=chi[n]; 
	  /* change dihedrals of other atoms that will change with this one */
	  for(i=9;i<=endSideAtomNumber;++i) 
	    if(resparam->atom[i].contactAtomNumber==resparam->atom[sideAtomNumber].contactAtomNumber && 
	       i!=sideAtomNumber)
	      {
		flag[i]=1; /* mark those atoms w/ changed dihedrals so they won't be changed when sideAtomNumber=i */
 		dihed2=dihedral+(resparam->atom[i].dihedral -(-1)*(resparam->atom[sideAtomNumber].dihedral));
		dum[i] = chi2xyz(&dum[resparam->atom[i].contactAtomNumber],
				&dum[resparam->atom[i].twobondAtomNumber],
				&dum[resparam->atom[i].dihedAtomNumber],
				&(resparam->atom[i].bondlength),
				&(resparam->atom[i].bondangle),
				&dihed2);
	      }
	}
      else
	  dihedral = resparam->atom[sideAtomNumber].dihedral;
      if(flag[sideAtomNumber]!=1) /* won't calc. coords for those atoms we already took care of above */ 
	dum[sideAtomNumber] = chi2xyz(&dum[resparam->atom[sideAtomNumber].contactAtomNumber],
				      &dum[resparam->atom[sideAtomNumber].twobondAtomNumber],
				      &dum[resparam->atom[sideAtomNumber].dihedAtomNumber],
				      &(resparam->atom[sideAtomNumber].bondlength),
				      &(resparam->atom[sideAtomNumber].bondangle),
				      &dihedral);
     
      side->atom[sideAtomNumber].atom_ptr =  &resparam->atom[sideAtomNumber];
    }
  side->atom[endSideAtomNumber+1].atom_ptr = NULL;
  for(sideAtomNumber=9;sideAtomNumber<=endSideAtomNumber;++sideAtomNumber)  /* copy sidechain atom coords from dum */
    side->atom[sideAtomNumber].coord = dum[sideAtomNumber];
  }
  else
	{
		build_a_Ligand(side, resparam, chi);
	}

  
}

/* extract BACKBONE from atom (array of mini_pdbATOM)
	bkbn must be allocated by the calling function  */

void extract_mini_to_bkbn(mini_pdbATOM *atom, BACKBONE bkbn[])
{  
  int i, j,chain_number;
  mini_pdbATOM *atom0;
  double len_C_CA, angle_C_CA_N, len_N_C, angle_N_C_CA, len_CA_N, angle_CA_N_C, dihed;
  div_t chain_id;

  len_C_CA=1.522;  angle_C_CA_N=111.10;
  len_N_C=1.335;   angle_N_C_CA=116.60; 
  len_CA_N=1.449;  angle_CA_N_C=121.90;
  
  atom0 = atom;
  
  /* get coordinates */
  ++atom; i=0; j=0;
  while(atom->seq_position!=ENDFLAG)
    {	
      if(atom->seq_position>i)
      {
	i = atom->seq_position;
	++j;
	bkbn[j].H_NTE = NULL;
	bkbn[j].O_CTE = NULL;
      }	
      if(strcmp(atom->atom_ptr->atomname,"N")==0)
	{ 
	  bkbn[j].N = atom->coord;
	  bkbn[j].seq_position=atom->seq_position;
	}
      if(strcmp(atom->atom_ptr->atomname,"CA")==0)
	bkbn[j].CA = atom->coord;
      if(strcmp(atom->atom_ptr->atomname,"C")==0)
	bkbn[j].C = atom->coord;
      if(strcmp(atom->atom_ptr->atomname,"O")==0)
	bkbn[j].O = atom->coord;
      if(strcmp(atom->atom_ptr->atomname,"CB")==0 || strcmp(atom->atom_ptr->atomname,"2HA")==0)
	bkbn[j].CB = atom->coord;
      if(strcmp(atom->atom_ptr->atomname,"H")==0 || strcmp(atom->atom_ptr->atomname,"1H")==0) 
	bkbn[j].H = atom->coord;
      if(strcmp(atom->atom_ptr->atomname,"HA")==0 || strcmp(atom->atom_ptr->atomname,"1HA")==0)
	bkbn[j].HA = atom->coord;
      
      if(strcmp(atom->atom_ptr->atomname,"1H")==0 || strcmp(atom->atom_ptr->atomname,"2H")==0  || 
		strcmp(atom->atom_ptr->atomname,"3H")==0)
      {   
	if(bkbn[j].H_NTE == NULL)
		bkbn[j].H_NTE = (CARTESIAN *)calloc(3, sizeof(CARTESIAN));
	if(strcmp(atom->atom_ptr->atomname,"2H")==0)   
	    bkbn[j].H_NTE[1] = atom->coord;
	else if(strcmp(atom->atom_ptr->atomname,"3H")==0)
	    bkbn[j].H_NTE[2] = atom->coord;     
      }
      if(strcmp(atom->atom_ptr->atomname,"OT")==0 || strcmp(atom->atom_ptr->atomname,"OXT")==0)
	{
	    bkbn[j].O_CTE = (CARTESIAN *)malloc(sizeof(CARTESIAN));
	    *bkbn[j].O_CTE = atom->coord;
	}
	
      ++atom; 
    }
    
  bkbn[j+1].seq_position=ENDFLAG;    
  
 	/* get dihedrals */
  i=1; chain_number = -10;
  while(bkbn[i].seq_position!=ENDFLAG)
    {
        chain_id = div(bkbn[i].seq_position, 1000);
	if(chain_number == chain_id.quot)
	    bkbn[i].phi = dihedral(&bkbn[i-1].C,&bkbn[i].N,&bkbn[i].CA,&bkbn[i].C);
        else
	{
	    bkbn[i].phi = 0;  /* marker for non-existant phi value for the first residue of a chain */       
	    chain_number = 	chain_id.quot;
	}

	if(bkbn[i+1].seq_position!=ENDFLAG)
	{    
	    bkbn[i].psi = dihedral(&bkbn[i].N,&bkbn[i].CA,&bkbn[i].C,&bkbn[i+1].N);
	    bkbn[i].omega = dihedral(&bkbn[i].CA,&bkbn[i].C,&bkbn[i+1].N,&bkbn[i+1].CA);
	}
	++i;
    }  
    
    /* generate bkbn[0]; coords for C and CA for a residue zero, if it existed */
    bkbn[0].omega = 180; bkbn[0].phi = 180; bkbn[0].psi = 180; 
    bkbn[0].O_CTE = NULL;  bkbn[0].H_NTE = NULL;
    dihed = 180.0;
    angle_CA_N_C = 180 + angle_CA_N_C;
    bkbn[0].C = chi2xyz(&bkbn[1].N, &bkbn[1].CA, &bkbn[1].C, &len_N_C, &angle_CA_N_C, &dihed);
    dihed = 180.0;
    angle_N_C_CA = 180 - angle_N_C_CA;
    bkbn[i].CA = chi2xyz(&bkbn[0].C, &bkbn[1].N, &bkbn[1].CA, &len_C_CA, &angle_N_C_CA, &dihed);

    atom = atom0;
}

/* extract BACKBONE from atom (array of pdbATOM)
	bkbn must be allocated by the calling function  */

void extract_bkbn(pdbATOM *atom, BACKBONE bkbn[])
{ 
  mini_pdbATOM *mini_pdb;
  
  mini_pdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM));
  
  pdbATOM_to_mini_pdbATOM(atom, mini_pdb);
  
  extract_mini_to_bkbn(mini_pdb, bkbn);
  
  free_memory(mini_pdb);
}

/* This function extracts the dihedrals from a pdbATOM array containing the atoms for a single residue
   The dihedrals are returned in chi array of double. resparam must be for the residuetype of interest.
   This is not very fast or efficient, but there is no need to be as it is done infrequently */

void extract_dihedrals(pdbATOM *pdb, double chi[], RESPARAM *resparam)
{
  int atomnumber, seq_position, ok_flag;
  int pdb_atom, pdb_contactatom, pdb_twobondatom, pdb_dihedatom;  /* indexes for pdbATOM elements */ 
  int chi_index;
  double *dummychi;
  extern int *MISSING_SIDECHAIN_LIST;


  for(chi_index=0;chi_index<=resparam->rotamerlib_ptr->numChi+1;++chi_index)
    chi[chi_index]=0; 

  if(resparam->ligand_flag==0)
  {
	if((*resparam).num_rotatable_bonds==0)
      		return;
  }
  else
  {
	if(resparam->rotamerlib_ptr->numChi==1)
	{
		extract_transform_matrix_for_ligand(pdb,chi,resparam);
		return;
	}
  }
	
  seq_position = pdb[1].seq_position;
  
  atomnumber=8;  /* atomnumbering starts w/ 4 for N; CB is 8 */
  for(chi_index=1;chi_index<=(*resparam).num_rotatable_bonds;++chi_index)
    {
      ok_flag = 1;
      
      while((*resparam).atom[atomnumber].dihedral>=0)
	++atomnumber;
      
      /* find atom in the pdbATOM array */
      pdb_atom = 1;
      while(strcmp(pdb[pdb_atom].atomname,(*resparam).atom[atomnumber].atomname)!=0 && pdb[pdb_atom].seq_position!=ENDFLAG)
	++pdb_atom;
    
      if(pdb[pdb_atom].seq_position == ENDFLAG)
      {
	if(QUIET_FLAG==0)
	{
	fprintf(stderr, "WARNING missing atom: %d %s %s\n", seq_position,  (*resparam).residuetype, (*resparam).atom[atomnumber].atomname);
	fprintf(stderr, "will rebuild with chi%d = 180.0\n", chi_index);
	}

	if(MISSING_SIDECHAIN_LIST == NULL)
	{
		MISSING_SIDECHAIN_LIST = (int *)calloc(MAX_RESIDUES,sizeof(int));
		MISSING_SIDECHAIN_LIST[1]=ENDFLAG; MISSING_SIDECHAIN_LIST[2]=ENDFLAG; 
	}
	ok_flag=1;
	while(MISSING_SIDECHAIN_LIST[ok_flag] != ENDFLAG)
		++ok_flag;
	MISSING_SIDECHAIN_LIST[ok_flag] = seq_position;
	++ok_flag;
	MISSING_SIDECHAIN_LIST[ok_flag] = ENDFLAG;

	ok_flag = 0;
      }	
      
      if(ok_flag == 1)
      {
	/* find contactatom in the pdbATOM array */
	pdb_contactatom = 1;
	while(strcmp(pdb[pdb_contactatom].atomname,
		   (*resparam).atom[(*resparam).atom[atomnumber].contactAtomNumber].atomname)!=0)
	    ++pdb_contactatom;

	/* find twobondatom in the pdbATOM array */
	pdb_twobondatom = 1;
	while(strcmp(pdb[pdb_twobondatom].atomname,
		   (*resparam).atom[(*resparam).atom[atomnumber].twobondAtomNumber].atomname)!=0)
	    ++pdb_twobondatom;

	/* dihedatom in the pdbATOM array */
	pdb_dihedatom = 1;
	while(strcmp(pdb[pdb_dihedatom].atomname,
		   (*resparam).atom[(*resparam).atom[atomnumber].dihedAtomNumber].atomname)!=0)
	    ++pdb_dihedatom;

	chi[chi_index] = dihedral(&pdb[pdb_atom].coord, &pdb[pdb_contactatom].coord, &pdb[pdb_twobondatom].coord, &pdb[pdb_dihedatom].coord);
      }
      else
	chi[chi_index] = 180.0;
	
      ++atomnumber; 
    }

    if(resparam->ligand_flag==1) /* ligand w/ tranform matrix */
    {
	dummychi = (double *)calloc(resparam->num_rotatable_bonds+2,sizeof(double));
	
	chi_index=1;
	while(chi_index <= resparam->num_rotatable_bonds)	/* copy bond diheds */
	{
		dummychi[chi_index] = chi[chi_index];
		++chi_index;
	}

	extract_transform_matrix_for_ligand(pdb,chi,resparam);
	
	/* put the bond dihedrals at the end of the chi array */
	chi_index=13;
	while(chi_index<=resparam->rotamerlib_ptr->numChi)
	{
		chi[chi_index] = dummychi[chi_index-12];
		++chi_index;
	}
		

	free_memory(dummychi);
    }	


}


/* returns the coordinates of the centroid for an array of mini_pdbATOM */
CARTESIAN get_centroid(mini_pdbATOM *pdb)
{
	int i,j;
	double one_over_i;
	CARTESIAN centroid;	

	CARTESIAN_zero(centroid);
	i=1; j=1;
	while(pdb[j].seq_position!=ENDFLAG)
	{
		if(pdb[j].atom_ptr->atom_ptr->atom_class!=ENDFLAG)
		{
			centroid = add_CARTESIAN(centroid, pdb[j].coord);
			++i;
		}
		++j;
	}
	i = i-1;

	one_over_i = 1.0/((double)i);
	centroid = scalar_multiply_CARTESIAN(one_over_i, centroid);

	return(centroid);
}


/* reduces an array pdb of pdbATOM to an array minipdb of mini_pdbATOM */

void pdbATOM_to_mini_pdbATOM(pdbATOM *pdb, mini_pdbATOM *minipdb)
{
  pdbATOM *dummypdb;
  mini_pdbATOM *dummymini;

  dummypdb = pdb;
  dummymini=minipdb;
  minipdb->seq_position=0;
  ++pdb; ++minipdb;

  while(pdb->seq_position!=ENDFLAG)
    {
      minipdb->coord = pdb->coord;
      minipdb->seq_position = pdb->seq_position;
      minipdb->atom_ptr = pdb->atom_ptr;
      minipdb->sasa = pdb->sasa; 
      minipdb->hbond_satisfied = pdb->hbond_satisfied;
      minipdb->born_radius = pdb->born_radius;
      minipdb->box=0;
      ++pdb; ++minipdb;
    }
  minipdb->sasa = ENDFLAG;
  minipdb->seq_position=ENDFLAG;
  minipdb->atom_ptr = NULL;
  pdb->sasa = ENDFLAG;
  pdb->seq_position=ENDFLAG;
  pdb->atom_ptr = NULL;
  strcpy(pdb->residuetype, "END");

  minipdb=dummymini;
  pdb=dummypdb;

}