#include "Vectors.h" // gives access to vectors and basic text utils

#ifndef __PDBOBJECT_H_INCLUDED__
#define __PDBOBJECT_H_INCLUDED__


class Hbond{
      public:
             
      int donor;
      int h;
      int acceptor;
      int base;
      double E;
      double Ed; //distance-only energy
      double angular; //angular factor
      
      Hbond(){
              //on creation all values are NOT possible atom IDs
              donor = -1;
              h = -1;
              acceptor = -1;
              base = -1;
              E=0.;
              Ed=0.;
              angular=1.;
      }
      string hbond2line();
      
};

class PhobicPair{
	public:
	
	int a1;
	int a2;
	PhobicPair(){
		a1 = -1;
		a2 = -1;
	}
	
};

class PDBConect{
      public:
             
      string recname; //1-6 = CONECT
      int atom1;
      vector< int > bond;
      PDBConect(){
                recname="      ";
                atom1 = 0;
      }
      
      string conect2line();
};


class PDBAtom{
  public:

  string recname; //1-6
  int serial; //7-11
  string name; //13-16
  string altloc; //17
  string resname; //18-20
  //21 blank
  string chainid; //22
  int seqnum; //23-26 #is the resi
  string insert; //27 #insertioncode
  //28-30 blank
  Vector pos; //31-38, 39-46, 47-54 #as r8.3
  double occ; //55-60 #occupancy as r6.2
  double b; //61-66 #bfactor as r6.2
  //67-72 blank
  string segment; //73-76 #segment id, left just
  string element; //77-78 #tail element symbol
  string charge; //79-80 #charge; not for us?

  int thing, file, model, chain, res, myID; //indices in master PDBThing

  bool isAtom, isHet, hasTer, isPolar, isMain, isAlpha;
  bool isStandard; //am I a standard residue?
  
  bool isDonor, isAcceptor, isChargedDonor, isChargedAcceptor;
  bool isPhobic; //suitable for phobic tether formation?
  //these bools are for potential hbond finding routines.
  //isChargedDonor will label positive NH+ nitrogen
  //isChargedAcceptor will label COO- oxygen
  
  int nH; //number of hydrogen covalent neighbours
  bool isImiN; //imadazolate nitrogens
  int inImiWith; //other imi nitrogen in ring
  bool isGuaniC; //guanidinium group C
  bool isGuaniN; //guanidinium group N
  int inGuaniWith; //id of guanidinium C
  bool isCofCOO; //tags the C of COO- groups
  
  int remote; //remoteness indicator for sidechains, not used yet
  int el; // let's translate string element to periodic table entry!
  int q; //translate string charge to numeric q, not used yet


  PDBAtom(){
    recname="      ";
    serial = 0;
    name = "    ";
    altloc = " ";
    resname = "   ";
    chainid = " ";
    seqnum = 0;
    insert = " ";
    pos = Vector(0.,0.,0.);
    occ = 1.;
    b = 0.;
    segment = "    ";
    element = "  ";
    charge = "  ";
    //indices
    thing = 0; file = 0; model = 0; chain = 0; res = 0; myID = 0;
    isAtom = false, isHet = false, hasTer = false;
    isPolar = false, isMain = false, isAlpha = false;
    isStandard = false; isDonor = false; isAcceptor = false;
    isChargedDonor = false; isChargedAcceptor = false;
    isPhobic = false;
    
    nH = 0;
    isImiN = false;
    inImiWith = -1;
    isGuaniC = false;
    isGuaniN = false;
    inGuaniWith = -1;
    isCofCOO = false;    
    
    remote = -2;
    el = -1;
    q=0;
  }

  string atom2line(); // returns string with PDB format of data
  string atom2id(); //returns string with PDB info only up to coords
  string atom2ter(); //returns string with TER card to follow atom
  string atom2verbose(); //chatty output for debug and reporting
  
  bool elH(); // returns true if hydrogen
  bool elO(); //returns true if oxygen
  bool elN(); //returns true if nitrogen
  bool elC(); //returns true if carbon
  bool elCHON(); //returns true if CHON
  bool elCONP(); //returns true if CONP (i.e. standard but not H or S
  bool elS(); // true if sulphur
  bool elSe(); //true if selenium
  bool elP(); //true if phosphorus
  
  //a bank of "special metal" cases plus chlorine ion
  bool elNa();
  bool elCl();
  bool elMg();
  bool elK();
  bool elCa();
  bool elMn();
  bool elFe();
  bool elCu();
  bool elZn();
  bool elHg();
  
  //a bank of categories for completeness and generic cases
  bool elAlkali();
  bool elAlkaliEarth();
  bool elNoble();
  bool elHalogen();
  bool elMetalloid();
  bool elTrans();
  bool elPostTrans();
  bool elMetal(); //catch the common metals
  

}; 

Vector relvec( PDBAtom &a, PDBAtom &b ); //relative vector between two atoms given
double drel( PDBAtom &a, PDBAtom &b ); //distance between two atoms given
bool isSameResidue( PDBAtom &a, PDBAtom &b ); // do they have the same .res index?
bool isResidueWithin( PDBAtom &a, PDBAtom &b, int c ); //same chain, residue difference up to c
bool isResidueAt( PDBAtom &a, PDBAtom &b, int c ); //same chain, residue difference exactly c
bool isAdjacentResidue( PDBAtom &a, PDBAtom &b ); //calls isResidueAt( a,b,1 );

class PDBResidue{
  public:
  
  vector< int > atom ; // serial numbers in master register of atoms
  //items stolen from atom entries, in case needed:
  string resname; //18-20
  //21 blank
  string chainid; //22
  int seqnum; //23-26 #is the resi
  string insert; //27 #insertioncode
  //indices:
  int thing, file, model, chain, myID; //indices in master PDBThing
  bool isStandard; // am I a standard residue?
  
  PDBResidue(){
    resname = "   ";
    chainid = " ";
    seqnum = 0;
    insert = " ";
    //indices
    thing = 0; file = 0; model = 0; chain = 0; myID = 0;
    isStandard = false;              
  }
  
  bool isHistidine(); //needed to trap a special case sometimes
  string res2id(); //to output some identification data in PDB format
                       
};

class PDBChain{
  public:
  
  vector< int > atom ; // serial numbers in master register of atoms
  vector< int > res; // serial in master residue array
  //stolen from atom entries, in case needed:
  string chainid; //22
  //indices:
  int thing, file, model, myID; //indices in master PDBThing
  
  PDBChain(){
    chainid = " ";
    //indices
    thing = 0; file = 0; model = 0; myID = 0;                 
  }
                       
};

class PDBModel{
  public:
  
  vector< int > atom ; // serial numbers in master register of atoms
  vector< int > res; // serial in master residue array
  vector< int > chain; // serial in master chain array
  int modelnum; // model number from PDB input
  bool explicitNumber; // is this an explicitly numbered model or not?
  //single-model files (no MODEL line) has explicit=false
  //indices:
  int thing, file, myID; //indices in master PDBThing
  vector< PDBConect > conect; //array of CONECT record data
  
  PDBModel(){
    modelnum = 0;
    explicitNumber = false;
    //indices
    thing = 0; file = 0; myID = 0;                 
  }
                       
};

class PDBFile{
  public:
  
  vector< int > atom ; // serial numbers in master register of atoms
  vector< int > res; // serial in master residue array
  vector< int > chain; // serial in master chain array
  vector< int > model; //serial in master model array
  string name; // actual filename!
  bool hasModels; // has multiple models?
  //indices:
  int thing, myID; //indices in master PDBThing
  
  PDBFile(){
    name="";
    hasModels = false;
    //indices
    thing = 0; myID = 0;                 
  }
                       
};

class PDBThing{
  public:
  
  vector< PDBAtom > atom ; // master register of atoms
  vector< PDBResidue > res; // master residue array
  vector< PDBChain > chain; // master chain array
  vector< PDBModel > model; // master model array
  vector< PDBFile > file; // master file array
  int myID; // my ID if I'm in an array of PDBThings!
  int myLabel; // numeric label container for program use
  //turtles all the way down
  
  //GRID things need to go here
  //remember it's not periodic!
  vector< vector< int > > grid; // coarse grid atom list
  double gridmin; // minimum spacing in coarse grid 
  int Nx, Ny, Nz; // number of cells in coarse grid
  double gridpad; //pad the edges
  Vector lowcorner; 
  Vector highcorner;

  vector< double > radius; // store some radii here - may depend on function
  vector< vector< int > > closeby; //list of atoms close to each atom per "near"
  vector< vector< int > > cov; //covalent bond list will go here
  vector< vector< int > > bars; // "bar" constraint count for cov constraints!
  vector< bool > sp2; // use sp2 tags to mark atoms for rigidity-unification
  vector< double > q; //structure q, do not confuse with int/string q in the Atom object
  
  vector< int > polarH; // polar hydrogens suitable for hbond information
  vector< Hbond > hb; //my hbond list
  vector< PhobicPair > phobe; //my hydrophobic tether list
  
  PDBThing(){
    myID = 0;
    myLabel = 0;
    gridmin = 5.0;
    Nx=1;
    Ny=1;
    Nz=1;
    gridpad = 6.0; //pad the edges
    lowcorner = Vector(0.,0.,0.);
    highcorner = Vector(0.,0.,0.);
  }
  
  void verboseAtoms(); // chatty output to screen                      
  void setData(); //runs over atoms and residues setting values
  void numberModels(); //sets Models to explicit and numbered
  bool writePDBFile( int whichfile, string filename ); //writes out a file as a file
  void initialiseGrid(); //size the grid, extract Nq values
  void fillGrid(); //place atoms in grid array
  int gridX( Vector& pos );
  int gridY( Vector& pos );
  int gridZ( Vector& pos );
  int gridBox( Vector &pos );
  int gridBox ( PDBAtom &atom ); //overloading
  vector< int > near ( PDBAtom &a, double within ); // get neighbour list
  vector< int > near ( Vector& pos, double within ); //get atoms near a point, overloading
  Vector rel( int i, int j ); // relative vector specified by indices in atom array
  double d( int i, int j ); // distance between atoms by indices
  double d2( int i, int j ); //squared distance between indices ( faster than d because no square root) 
  
  void assignBurialRadii(); //radii for "sexton" burial distance calculation
  void sexton(double prober); //carries out "burial" calculation on heavy atoms
  
  bool writeMap(string filename); //writes out the map between internal index and PDB data
  void getClosebyForCov(); // run covalent-bond identifications
  void clearCov(); //pre-blank the cov array
  void conectToCov(); //work over cov entries as needed
  int IDfromSerial( int start, int N, int serial ); //find the atom ID for a given serial
  void completeCov(); //invoke after getClosebyForCov and conectToCov to finish covalent bonding
  bool covBonded( int i, int j ); // are atoms i,j covalently bonded in this Thing?
  bool writeCov( string filename ); //write out cov file in format matching "map"
  void sp2Check(); //relate bonding to status
  void barCheck(); //assign bars between adjacent sp2s
  void fillPolarH(); //feed the polar H list

  void detectGuani();//scans structure labelling guanidinium groups (important for arginine salt bridges).
  void detectImi(); //scans structure for imidazolate groups (important for histidines).
  void detectCOO(); //scans structure for COO groups (not triggered by COOH!) for salt bridge acceptors.
  void assignPolarityLabels(); //boolean labels for PDBAtom objects
  void seekHbondsByPolarityLabels( int energyFunction = 0, double ECutOff = 0.0 ); //use for hbond assignment, choose energy rating method
  void assignPhobicRadii(); //assign FIRST-like radii for hydrophobic tether finding: C1.7, S1.8
  void assignContactRadii(); //assign "vdW" radii to pass into mobility routine
  void seekPhobicPairs( double margin ); //find matching pairs of phobic atoms
  void filterPhobicPairs(); //remove redundancies from the list of phobic pairs
  void filterHbondsByDistance(); // implement energy version later, then use both
  //void getHbondEnergiesBySimpleLinear(); //use the simple function to assess energy based on distances
  //void energyBySimpleFunctionA( Hbond &hb ); //simple distance function for hbond based on types
  //void energyByFunctionB( Hbond &hb ); //developing function using safe LJ and corrected angular terms
  
  string verboseHbondLine( Hbond &hb); //report atom info in an Hbond object
  string terseHbondLine( Hbond &hb); //report minimal info: hydrogen, acceptor, energy
  string fullHbondLine( Hbond &hb ); //report hbond info with distance and angle terms
  string fullPhobicPairLine( PhobicPair &pp ); //report phobic info with distance
  string fullCovalentLine( int i, int j ); //reports the covalent bond in cov(i)(j) 
  
  bool writeHydrogenBonds( string filename );
  bool writePhobicTethers( string filename , int Cryo );
  
  Vector normalDH( Hbond &hb ); // return the normal to an sp2 donor group
  Vector normalAB( Hbond &hb ); //return the normal to an sp2 acceptor group
  double planeNormalCosine( Hbond &hb ); //return a nonplanarity cosine for sp2 donor and acceptor groups
  double bondPlaneCosine( Hbond &hb ); //return outofplane measure for DH bond versus sp2 acceptor group plane
  
  double LJenergy( Hbond &hb ); //returns the FIRST-style LJ energy based on bond D-A distance
  double safeLJenergy( Hbond &hb ); //returns LJ energy without repulsize penalty at short distance
  double cosSqTheta( Hbond &hb ); // returns cos-squared of the DHA angle
  double FIRSTprefactor( Hbond &hb ); //returns cos squared DHA angle times exp ( - ( pi-theta)^6 ) as in FIRST, bafflingly
  double cosSqPhiPhi0( Hbond &hb, double phi0 ); // returns penalty on HAB angle compared to phi0 ideal in degrees
  double cosSqPhiSp2( Hbond &hb ); //returns appropriate phi penalty for sp2 acceptor
  
  double cosSqPsi( Hbond &hb ); // returns Drieding-style planarity penalty as in FIRST
  double cosSqOmega( Hbond &hb ); //returns my bond/planarity penalty on DH bond vector compared to base plane
  
  void energyFIRST( Hbond &hb ); // gives hb its FIRST energy
  void safetyFIRST( Hbond &hb ); // FIRST with safeLJ
  void energySECOND( Hbond &hb ); //my modified energy function using safeLJ and angle handling
  void energyDREIDING( Hbond &hb ); //using only cos^4 DHA
  
  double FIRSTsp3sp3( Hbond &hb ); //FIRST case 3-3 angular factor 
  double FIRSTsp2sp3( Hbond &hb ); //FIRST case 2-3 angular factor
  double FIRSTsp3sp2( Hbond &hb ); //FIRST case 3-2 angular factor  
  double FIRSTsp2sp2( Hbond &hb ); //FIRST case 2-2 angular factor
  double safetyFIRSTsp3sp3( Hbond &hb ); //FIRST case 3-3 angular factor 
  double safetyFIRSTsp2sp3( Hbond &hb ); //FIRST case 2-3 angular factor
  double safetyFIRSTsp3sp2( Hbond &hb ); //FIRST case 3-2 angular factor  
  double safetyFIRSTsp2sp2( Hbond &hb ); //FIRST case 2-2 angular factor
  double my_sp3sp3( Hbond &hb ); //my case 3-3 angular factor 
  double my_sp2sp3( Hbond &hb ); //my case 2-3 angular factor
  double my_sp3sp2( Hbond &hb ); //my case 3-2 angular factor  
  double my_sp2sp2( Hbond &hb ); //my case 2-2 angular factor  

  //now deprecating all these! moved to Dead folder as DSSProutinesDidntWork.txt :
  //void initialiseCharge(); // size the q array
  //void assignBackboneDSSP(); // DSSP CONH charges
  //void assignChargeLikeDSSP(); //pseudo-DSSP values for sidechains
  //void assignNominalCharges(); //just assign some +/- small values to make hbond search possible
  //void seekHbondsByChargeGeometry(); //hydrogen bond search using charge info and distances
  //void getHbondEnergiesByDSSP(); //use DSSP-like calculation using atomic charges, optional verbose reporting
  //void energyByDSSP( Hbond &hb ); //evaluate an hbond energy, with optional verbose reporting
  
  //void compareDistanceWells(); //just a debugging function to compare outputs
  
  Vector residueCentralPoint( int i ); //returns the "central" point of residue i
  Vector residueAlphaPoint( int i ); //returns the CA position of a standard residue
  vector< Vector > nodesForENM(); // return one node per residue
  bool writePseudoPDBforENM( string filename ); //make a nodesForENM object and write it as a pseudo PDB
  bool writeXYZforENM( string filename ); //make a nodesForENM object and write it as an XYZ list
  bool writePseudoPDBforENM( string filename, vector< Vector> &nodes ); //overload - nodes provided
  //make a nodesForENM object and write it as a pseudo PDB
  bool writeXYZforENM( string filename, vector< Vector> &nodes ); //overload - nodes provided
  //make a nodesForENM object and write it as an XYZ list

  
  vector< Vector > allXYZs(); //all the positions extracted from atom pos
  vector< Vector > fileXYZs( int i ); //all the positions from file i
  vector< Vector > modelXYZs( int i ); //all the positions from model i
  vector< Vector > chainXYZs( int i ); //all the positions from chain i
  vector< Vector > residueXYZs( int i ); //all the positions from residue i
  
  vector< vector< int > > smallRC(); //atoms grouped into small rigid bodies using Garibaldi logic
  bool writeRC( string filename, vector< vector< int > > &RC );

  bool writeFrame( string rider ); //write out my files as file.rider.pdb
 
};

const PDBAtom nullPDBAtom;
const PDBResidue nullPDBResidue;
const PDBChain nullPDBChain;
const PDBModel nullPDBModel;
const PDBFile nullPDBFile;
const PDBThing nullPDBThing;
const PDBConect nullPDBConect;

bool pullData(PDBAtom& atom, string& line);
bool pullConect( PDBConect& con, string& line );
bool readPDBFile(PDBThing& pdbt, string filename); //reads file into a fresh clean Thing
void insertAintoB( PDBThing& A, PDBThing& B ); //sticks contents of A on top of contents of B

int species2number( string species ); // return the atomic number for species
bool isStandardResidue( string resname); //return true if in standard list

double LJ( double r, double r0, double a, double depth ); //Lennard Jones potential for general use
double cubicWell( double r, double r1, double r2, double depth ); //flat-bottom, cubic-side well
double safeLJ( double r, double r0, double a, double depth ); //Lennard Jones potential with short-safe handling

double cosThetaTriangle( double a, double b, double c ); //cosine rule on triangle
double cosSqExp6( double theta ); // FIRST prefactor with exponential dependency


#endif