#include "PDBobject.h"


string Hbond::hbond2line(){
       stringstream record;
       record << "HBOND ids: " << donor+1 << " " << h+1 << " " << acceptor+1 << " " << base+1;
       record << " Ed: " << Ed << " angular: " << angular << " E: " << E;
       return record.str();
}

string PDBConect::conect2line(){
  stringstream record;
  
  record << setw(6) << recname;
  record << setw(5) << atom1;
  for ( int i = 0 ; i < bond.size(); i++ ){
      int j = bond.at(i);
      if ( j < 1 ){
         record << setw(5) << "     ";     
      }
      else{
           record << setw(5) << j;
      } 
  }
  //allows for blanks if CONECT is used properly, which it probably isn't
  return record.str();       
}

string PDBAtom::atom2line(){
  stringstream record;

  record << setw(6) << recname;
  record << setw(5) << serial << " ";
  record << setw(4) << name;
  record << altloc;
  record << setw(3) << resname << " ";
  record << chainid;
  record << setw(4) << seqnum;
  record << insert << "   ";
  record << fixed << showpoint << setw(8) << setprecision(3) << pos.x;
  record << fixed << showpoint << setw(8) << setprecision(3) << pos.y;
  record << fixed << showpoint << setw(8) << setprecision(3) << pos.z;
  record << fixed << showpoint << setw(6) << setprecision(2) << occ;
  record << fixed << showpoint << setw(6) << setprecision(2) << b;
  record << "      ";
  record << setw(4) << segment;
  record << setw(2) << element;
  record << setw(2) << charge;
  
  return record.str();      
         
}

string PDBAtom::atom2id(){
  stringstream record;

  record << setw(6) << recname;
  record << setw(5) << serial << " ";
  record << setw(4) << name;
  record << altloc;
  record << setw(3) << resname << " ";
  record << chainid;
  record << setw(4) << seqnum;
  record << insert;
  
  return record.str();      
         
}

string PDBResidue::res2id(){
  stringstream record;

  record << setw(3) << resname << " ";
  record << chainid;
  record << setw(4) << seqnum;
  record << insert;  
  return record.str();      
         
}

string PDBAtom::atom2ter(){
  stringstream record;

  record << setw(6) << "TER   ";
  record << setw(5) << serial+1 << " ";
  record << setw(4) << name;
  record << altloc;
  record << setw(3) << resname << " ";
  record << chainid;
  record << setw(4) << seqnum;
  record << insert;
  
  return record.str();       
}



string PDBAtom::atom2verbose(){
  stringstream verb;
  
  verb << atom2line() << endl; //outputting the PDB line itself...
  verb << "Serial " << serial << " name " << name << " alt " << altloc;
  verb << " resname " << resname << " chainid " << chainid;
  verb << " seqnum " << seqnum;
  if ( insert != " " ) verb << " insert " << insert << endl;
  verb << endl;
  verb << " X " << pos.x << " Y " << pos.y << " Z " << pos.z;
  verb << " occ " << occ << " b " << b << " segment " << segment;
  verb << " element " << element;
  if ( charge != "  " ) verb << " charge " << charge <<endl;
  verb << endl;

  verb << "Additional properties:" << endl;
  verb << "Thing " << thing << " File " << file << " Model " << model;
  verb << " Chain " << chain << " Res " << res << " MyID " << myID << endl;
  if ( isAtom ) verb << "isAtom ";
  if ( isHet ) verb << "isHet ";
  if ( hasTer ) verb << "hasTer ";
  if ( el == 1 ) verb << "isHydrogen ";
  if ( isPolar ) verb << "isPolar ";
  if ( el == 16 ) verb << "isSulfur ";
  if ( isMain ) verb << "isMain ";
  if ( !isMain ) verb << "notMain ";
  if ( isAlpha ) verb << "isAlpha ";
  if ( isStandard ) verb << "isStandard ";
  if ( !isStandard ) verb << "NonStandard ";
  verb << endl;
  
  verb << "Remote " << remote << " El " << el << " q " << q << endl << endl;

  return verb.str();
}

bool PDBAtom::elH(){
     if ( el == 1 ) return true;
     return false;
}

bool PDBAtom::elC(){
     if ( el == 6 ) return true;
     return false;
}

bool PDBAtom::elO(){
     if ( el == 8 ) return true;
     return false;
}

bool PDBAtom::elN(){
     if ( el == 7 ) return true;
     return false;
}

bool PDBAtom::elCHON(){
     if ( el == 1 || el == 6 || el == 7 || el == 8 ) return true;
     return false;
}

bool PDBAtom::elCONP(){
     if ( el == 6 || el == 7 || el == 8 || el == 15 ) return true;
     return false;
}

bool PDBAtom::elS(){
     if ( el == 16 ) return true;
     return false;
}

bool PDBAtom::elSe(){
     if ( el == 34 ) return true;
     return false;
}

bool PDBAtom::elP(){
     if ( el == 15 ) return true;
     return false;
}

bool PDBAtom::elNa(){
     if ( el == 11 ) return true;
     return false;
}

bool PDBAtom::elCl(){
     if ( el == 17 ) return true;
     return false;     
}

bool PDBAtom::elMg(){
     if ( el == 12 ) return true;
     return false;     
}

bool PDBAtom::elK(){
     if ( el == 19 ) return true;
     return false;     
}

bool PDBAtom::elCa(){
     if ( el == 20 ) return true;
     return false;     
}

bool PDBAtom::elMn(){
     if ( el == 25 ) return true;
     return false;     
}

bool PDBAtom::elFe(){
     if ( el == 26 ) return true;
     return false;     
}

bool PDBAtom::elCu(){
     if ( el == 29 ) return true;
     return false;     
}

bool PDBAtom::elZn(){
     if ( el == 30 ) return true;
     return false;     
}

bool PDBAtom::elHg(){
     if ( el == 80 ) return true;
     return false;     
}

bool PDBAtom::elAlkali(){
     if ( el == 3 || el == 11 || el == 19 || el == 37 || el == 55 || el == 87 ) return true;
     return false;     
}

bool PDBAtom::elAlkaliEarth(){
     if ( el == 4 || el == 12 || el == 20 || el == 38 || el == 56 || el == 88 ) return true;
     return false;     
}

bool PDBAtom::elNoble(){
     if ( el == 2 || el == 10 || el == 18 || el == 36 || el == 54 || el == 86 ) return true;
     return false;     
}

bool PDBAtom::elHalogen(){
     if ( el == 9 || el == 17 || el == 35 || el == 53 || el == 85 ) return true;
     return false;     
}

bool PDBAtom::elMetalloid(){
     if ( el == 5 || el == 14 || el == 32 || el == 33 || el == 51 || el == 52 || el == 84 ) return true;
     return false;     
}

bool PDBAtom::elPostTrans(){
     if ( el == 13 || el == 31 || el == 49 || el == 50 || el == 81 || el == 82 || el == 83 ) return true;
     return false;     
}

bool PDBAtom::elTrans(){
     if ( el > 20 && el< 31 ) return true; //first row trans
     if ( el > 38 && el < 49 ) return true; // second row trans
     if ( el > 56 && el < 81 ) return true; // third row trans
     if ( el > 88 ) return true; // fourth row trans
     return false;     
}


bool PDBAtom::elMetal(){
     //return true on general metals
     if ( elAlkali() ) return true;
     if ( elTrans() ) return true;
     if ( elAlkaliEarth() ) return true;
}


//pulldata gets PDB information out of a line and returns true if successful.
//returns false if too short or broken
bool pullData (PDBAtom& atom, string& line){
     int linelen = line.length();
     if ( linelen < 78 ){
        cerr << "Too short PDB line: " << line << endl;
        //elements must be missing at least
        return false;     
     }
     
     string temp;
     temp.assign( line, 0,6);
     atom.recname = temp;
     temp.assign( line, 6, 5);
     atom.serial = atoi( temp.c_str() );
     temp.assign( line, 12, 4 );
     atom.name = temp;
     temp.assign( line, 16, 1 );
     atom.altloc = temp; //one char
     temp.assign( line, 17,3);
     atom.resname = temp;
     temp.assign( line, 21, 1 );
     atom.chainid = temp;
     temp.assign( line, 22, 4);
     atom.seqnum = atoi( temp.c_str() );
     temp.assign( line, 26, 1 );
     atom.insert = temp;
     temp.assign( line, 30, 8 );
     atom.pos.x = atof( temp.c_str() );
     temp.assign( line, 38, 8 );
     atom.pos.y = atof( temp.c_str() );
     temp.assign( line, 46, 8 );
     atom.pos.z = atof( temp.c_str() );
     temp.assign( line, 54, 6);
     atom.occ = atof( temp.c_str() );
     temp.assign( line, 60, 6 );
     atom.b = atof( temp.c_str() );
     temp.assign( line, 72, 4 );
     atom.segment = temp;
     temp.assign( line, 76, 2 );
     atom.element = temp;
     if ( linelen < 79 ) return true; // might have clipped off nonexistent charge column
     temp.assign( line, 78, 2 );
     atom.charge = temp;
     return true;
}

bool pullConect( PDBConect& con, string& line ){
     int linelen = line.length();
     if ( linelen < 16 ) return false; //too short to contain anything
     string temp;
     temp.assign(line, 0, 6);
     if ( temp != "CONECT" ){
          cerr << temp << " isn't CONECT in pullConect?" << endl;
          return false;
     }
     con.recname = temp;
     temp.assign(line, 6, 5);
     con.atom1 = atoi( temp.c_str() ); // first in line;
     //now a loop to get other atoms...
     con.bond.clear();
     for ( int i = 0; i < 10 ; i++ ){
         //set up start and finish field positions
         int s = 11 + 5*i; //12, 17, 22... minus 1 each time
         int f = 15 + 5*i; //16, 21, 26... minus 1 each time
         //know when to stop...
         if ( linelen > f ){
              temp.assign( line, s, 5); // a field
              int a = atoi( temp.c_str() );
              if ( a == 0 ) continue; //skip blanks which become zero!
              con.bond.push_back( a ); // value goes in
         }
         else{
              //no more length
              break; //leave loop     
         }
     }
     return true; //all done
}


bool readPDBFile(PDBThing& pdbt, string filename){
     bool status = false; // return status success/fail of read

     bool explicitModels = false;
     bool alreadyReading = false;// use to trap MODEL after ATOMs have already started!
     bool inExplicit = false; //use to check for matching MODEL//ENDMDL
     
     pdbt = nullPDBThing; // clear everything at the start
     //some null objects for reading into
     PDBFile infile = nullPDBFile;
     PDBModel inmodel = nullPDBModel;
     PDBChain inchain = nullPDBChain;
     PDBResidue inres = nullPDBResidue;
     PDBAtom inatom = nullPDBAtom;
     PDBAtom lastatom = inatom; //we will maintain lastatom to detect new (res|chain)
     PDBConect incon = nullPDBConect;
     
     int runmodel = 0; //running indices
     int runchain = 0;
     int runres = 0;
     int runatom = 0;
     //remember, file is file 0 in pdbt by construction
     pdbt.file.push_back( infile );
     pdbt.file.at(0).myID = 0;
     pdbt.file.at(0).name = filename;

     infile.name = filename;    
     
     ifstream inf;
     inf.open( filename.c_str(), ios::in );
     if ( inf.fail() ){
          cerr << "ERROR: cannot open file " << filename << endl;
          inf.close();
          return false;
     }
     
     string foo; // line for reading
     while( getline( inf, foo ) ){
            //line parsing time!
            //cerr << foo << endl; //just echo line first
            int foolen = foo.length();
            if ( foolen < 6 ) continue; //not even a record here!
            
            string temp;
            temp.assign( foo, 0, 6 ); // the record name...
            
            //catch error cases first
            if ( temp == "ENDMDL" && !inExplicit ){
                 cerr << "ERROR: read ENDMDL when not in a MODEL." << endl;
                 inf.close();
                 return false; // bad!
            }
            if ( temp == "CONECT" && !alreadyReading ){
                 cerr << "ERROR: read CONECT before meeting atoms." << endl;
                 inf.close();
                 return false;
            }
            if ( temp == "TER   " && !alreadyReading ){
                 cerr << "ERROR: read TER before meeting atoms." << endl;
                 inf.close();
                 return false;
            }
            if ( temp == "MODEL " && inExplicit ){
                 cerr << "ERROR: read MODEL when already in MODEL." << endl;
                 inf.close();
                 return false;
            }
            if ( temp == "MODEL " && alreadyReading && !explicitModels ){
                 cerr << "ERROR: read MODEL when already started reading atoms." << endl;
                 inf.close();
                 return false;
            }
            
            if ( temp == "ENDMDL" && inExplicit ){
                 inExplicit = false; // just left a MODEL
            }
            if ( temp == "MODEL " ){
                 explicitModels = true;
                 inExplicit = true;
                 if ( foolen < 11 ){
                      //bad MODEL record
                      cerr << "ERROR: MODEL record line too short." << endl;
                      cerr << foo << endl;
                      cerr << "MODEL should have serial in characters 11-14." << endl;
                      inf.close();
                      return false;
                 }
                 string temp2;
                 temp2.assign( foo, 10, 4);
                 int ser = atoi( temp2.c_str() );
                 inmodel = nullPDBModel; // new blank Model
                 inmodel.modelnum = ser; //model number in record
                 inmodel.explicitNumber = true; //explicitly numbered
                 cerr << "Start of MODEL     " << setw(4) << ser << endl;
                 pdbt.model.push_back( inmodel ); // on master array
                 runmodel = pdbt.model.size() -1; //its serial
                 cerr << "This MODEL is entry " << runmodel << " in master array." << endl;
                 pdbt.file.at(0).model.push_back( runmodel );
                 pdbt.file.at(0).hasModels = true;
                 pdbt.model.at(runmodel).file = 0;
                 pdbt.model.at(runmodel).myID = runmodel;// I know who I am
            }
            if ( temp == "ATOM  " || temp == "HETATM" ){
                 if ( !alreadyReading ){
                      //cerr << "Met first atom: " << foo << endl;
                      cerr << "Starting to read atoms." << endl;
                      alreadyReading = true;
                      
                      if ( !explicitModels ){
                         //set up a blank MODEL with EXPLICIT = false;
                         inmodel = nullPDBModel; // new blank Model
                         inmodel.modelnum = -1; //no model number in record
                         inmodel.explicitNumber = false; //not explicitly numbered
                         cerr << "Single-Model file." << endl;
                         pdbt.model.push_back( inmodel ); // on master array
                         runmodel = 0; //its serial
                         pdbt.file.at(0).model.push_back( runmodel );
                         pdbt.file.at(0).hasModels = false;
                         pdbt.model.at(runmodel).file = 0;
                         pdbt.model.at(runmodel).myID = 0;
                      }                      
                      
                 }
                 //set up lastatom for previous case...
                 lastatom = inatom;
                 inatom = nullPDBAtom; //new blank
                 //try to read atom data
                 bool goodatom = false;
                 goodatom = pullData( inatom, foo );
                 if (!goodatom ){
                               cerr << "ERROR: bad atom read: " << foo << endl;
                               inf.close();
                               return false;
                 }
                 //okay atom data is in there now...
                 //test for (i) new chain (ii) new residue
                 bool changeChain = false;
                 bool changeRes = false;
                 
                 if ( inatom.chainid != lastatom.chainid ){
                      changeChain = true;
                      changeRes = true;
                      //new chain, hence also new res
                 }
                 if ( inatom.seqnum != lastatom.seqnum ){
                      changeRes = true;
                      //new sequence number
                 }
                 if ( inatom.insert != lastatom.insert ){
                      changeRes = true;
                      //new insertion code
                 }                 
                 
                 if ( changeChain ){
                      inchain = nullPDBChain;
                      inchain.chainid = inatom.chainid;
                      inchain.file = 0;
                      //current model:
                      runmodel = pdbt.model.size() -1;
                      runchain = pdbt.chain.size(); //not -1, we're adding here
                      inchain.model = runmodel;
                      inchain.myID = runchain;
                      
                      pdbt.chain.push_back( inchain ); //new chain entry
                      pdbt.file.at(0).chain.push_back( runchain);
                      pdbt.model.at(runmodel).chain.push_back( runchain );
                      
                 }
                 if ( changeRes ){
                      inres = nullPDBResidue;
                      inres.chainid = inatom.chainid;
                      inres.resname = inatom.resname;
                      inres.seqnum = inatom.seqnum;
                      inres.insert = inatom.insert;
                      inres.file = 0;
                      
                      runmodel = pdbt.model.size() -1;
                      runchain = pdbt.chain.size() -1;
                      runres = pdbt.res.size();
                      
                      inres.myID = runres;
                      inres.model = runmodel;
                      inres.chain = runchain;
                      
                      pdbt.res.push_back( inres );
                      pdbt.file.at(0).res.push_back( runres );
                      pdbt.model.at( runmodel ).res.push_back( runres );
                      pdbt.chain.at( runchain ).res.push_back( runres );
                 }
                 
                 //finally actually put the atom into place!
                 runmodel = pdbt.model.size() - 1;
                 runchain = pdbt.chain.size() -1;
                 runres = pdbt.res.size() -1;
                 runatom = pdbt.atom.size(); //inserting here
                 
                 inatom.file = 0;
                 inatom.model = runmodel;
                 inatom.chain = runchain;
                 inatom.res = runres;
                 inatom.myID = runatom;
                 
                 pdbt.atom.push_back( inatom );
                 pdbt.file.at(0).atom.push_back( runatom );
                 pdbt.model.at( runmodel ).atom.push_back( runatom );
                 pdbt.chain.at( runchain ).atom.push_back( runatom );
                 pdbt.res.at( runres ).atom.push_back( runatom );
                 
            }
            //TER case
            if ( temp == "TER   " ){
                 runatom = pdbt.atom.size() -1;
                 pdbt.atom.at( runatom).hasTer = true; //label previous atom with ter
            }
            //CONECT case and we're done!
            if ( temp == "CONECT" ){
                 //extract connect information
                 incon = nullPDBConect;
                 bool goodcon = false;
                 goodcon = pullConect( incon, foo );
                 if (!goodcon ){
                              cerr << "Bad CONECT line: " << foo << endl;
                              inf.close();
                              return false;
                 }
                 runmodel = pdbt.model.size() -1;
                 pdbt.model.at( runmodel ).conect.push_back( incon );
                 //cerr << "Read CONECT information: " << incon.conect2line() << endl;
            }
            
            
     }
     
     inf.close(); //file closed

     return true; //return good at the end           
}

//sticks contents of A on top of contents of B
void insertAintoB( PDBThing& A, PDBThing& B ){
     //first collect the relevant values for Thing B
     int nbatom, nbres, nbchain, nbmodel, nbfile;
     //insert A's atoms on top of B's;
     //likewise res, chain, model, file;
     nbatom = B.atom.size();
     nbres = B.res.size();
     nbchain = B.chain.size();
     nbmodel = B.chain.size();
     nbfile = B.file.size();

     cerr << "Starting insertion: B values: ";
     cerr << nbatom << " ";
     cerr << nbres << " ";
     cerr << nbchain << " ";
     cerr << nbmodel << " ";
     cerr << nbfile << endl;
    
     //all serial references update by +nb values
     //these references are the file model chain res myID uprefs;
     //and the contents of the model chain res atom array downrefs;
     
     //meanwhile thing must be set to B's thing value
     int bthing = B.myID;
     PDBAtom aatom; //use to transfer without damaging A
     PDBResidue ares;
     PDBChain achain;
     PDBModel amodel;
     PDBFile afile;
     
     for ( int i = 0; i < A.atom.size(); i++ ){
         aatom = A.atom.at(i);
         aatom.thing = bthing;
         aatom.file += nbfile;
         aatom.model += nbmodel;
         aatom.chain += nbchain;
         aatom.res += nbres;
         aatom.myID += nbatom;
         B.atom.push_back( aatom );
     }
     for ( int i = 0; i < A.res.size(); i++ ){
         ares = A.res.at(i);
         ares.thing = bthing;
         ares.file += nbfile;
         ares.model += nbmodel;
         ares.chain += nbchain;
         ares.myID += nbres;
         for ( int j = 0; j < ares.atom.size(); j++ ){
             ares.atom.at(j) += nbatom; // update atom numbers
         }
         B.res.push_back( ares );
     }     
     for ( int i = 0; i < A.chain.size(); i++ ){
         achain = A.chain.at(i);
         achain.thing = bthing;
         achain.file += nbfile;
         achain.model += nbmodel;
         achain.myID += nbchain;
         for ( int j = 0; j < achain.res.size(); j++ ){
             achain.res.at(j) += nbres; // update res numbers
         }         
         for ( int j = 0; j < achain.atom.size(); j++ ){
             achain.atom.at(j) += nbatom; // update atom numbers
         }
         B.chain.push_back( achain );
     }     
     for ( int i = 0; i < A.model.size(); i++ ){
         amodel = A.model.at(i);
         amodel.thing = bthing;
         amodel.file += nbfile;
         amodel.myID += nbmodel;
         
         for ( int j = 0; j < amodel.chain.size(); j++ ){
             amodel.chain.at(j) += nbchain; // update res numbers
         }
         for ( int j = 0; j < amodel.res.size(); j++ ){
             amodel.res.at(j) += nbres; // update res numbers
         }         
         for ( int j = 0; j < amodel.atom.size(); j++ ){
             amodel.atom.at(j) += nbatom; // update atom numbers
         }
         B.model.push_back( amodel );
     } 
     for ( int i = 0; i < A.file.size(); i++ ){
         afile = A.file.at(i);
         afile.thing = bthing;
         afile.myID += nbfile;
         
         for ( int j = 0; j < afile.model.size(); j++ ){
             afile.model.at(j) += nbmodel; // update res numbers
         }         
         for ( int j = 0; j < afile.chain.size(); j++ ){
             afile.chain.at(j) += nbchain; // update res numbers
         }
         for ( int j = 0; j < afile.res.size(); j++ ){
             afile.res.at(j) += nbres; // update res numbers
         }         
         for ( int j = 0; j < afile.atom.size(); j++ ){
             afile.atom.at(j) += nbatom; // update atom numbers
         }
         B.file.push_back( afile );
     }

     nbatom = B.atom.size();
     nbres = B.res.size();
     nbchain = B.chain.size();
     nbmodel = B.chain.size();
     nbfile = B.file.size();

     cerr << "Ending insertion: B values: ";
     cerr << nbatom << " ";
     cerr << nbres << " ";
     cerr << nbchain << " ";
     cerr << nbmodel << " ";
     cerr << nbfile << endl;     
     
}

void PDBThing::numberModels(){
     cerr << "Running numberModels." << endl;
     for ( int i = 0 ; i < file.size() ; i++ ){
         cerr << "File " << i+1 << "( " << file.at(i).name << ")" << endl;
         if ( file.at(i).model.size() == 1 ){              
               int k = file.at(i).model.at(0);
               model.at(k).explicitNumber = false;   
               cerr << "Single model file." << endl;
         }
         else{
              for ( int j = 0 ; j < file.at(i).model.size() ; j++ ){
                  int k = file.at(i).model.at(j);
                  model.at(k).explicitNumber = true;
                  model.at(k).modelnum = j+1;
                  cerr << "Model " << j+1 << endl;
              }
         }
     }
     cerr << "Finished numberModels." << endl;
}


bool PDBThing::writePDBFile( int whichfile, string filename ){
     //open the file for writing;
     ofstream outf;
     outf.open( filename.c_str(), ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file " << filename << endl;
          return false;
     }
     //loop over the entities in the given File
     //write out, taking MODEL information into account
     //outer loop is over Models, just trip on Explicit
     for ( int i = 0 ; i < file.at(whichfile).model.size() ; i++ ){
         int k = file.at(whichfile).model.at(i);
         if ( model.at(k).explicitNumber ){
              outf << "MODEL     " << setw(4) << model.at(k).modelnum << endl; //MODEL start
         }
         for ( int j = 0 ; j < model.at(k).atom.size(); j++ ){
             int p = model.at(k).atom.at(j);
             string aline = atom.at(p).atom2line();
             outf << aline << endl;
             if ( atom.at(p).hasTer ){
                  aline = atom.at(p).atom2ter();
                  outf << aline << endl;
             }
         }
         for ( int j = 0 ; j < model.at(k).conect.size(); j++ ){
             string aline = model.at(k).conect.at(j).conect2line();
             outf << aline << endl;
         }
         if ( model.at(k).explicitNumber ){
              outf << "ENDMDL" << endl; //MODEL end
         }         
         
     }
     outf << "END   " << endl;
     outf.close(); //closed stream at the end
     return true;
}


void PDBThing::verboseAtoms(){
     for ( int i = 0 ; i < atom.size(); i++ ){
         cerr << atom.at(i).atom2verbose();
     }
}

int species2number( string species ){
    chomp( species ); //whitespace handling
    to_upper( species ); //uppering
    if ( species == "H" || species == "D" ){
         return 1;
    }
    else if ( species == "C" ){
         return 6;
    }
    else if ( species == "O" ){
         return 8;
    }
    else if ( species == "N" ){
         return 7;
    }
    else if ( species == "S" ){
         return 16;
    }
    else if ( species == "NA" ){
         return 11;
    }
    else if ( species == "CL" ){
         return 17;
    }
    else if ( species == "P" ){
         return 15;
    }
    else if ( species == "HE" ){
         return 2;
    }
    else if ( species == "LI" ){
         return 3;
    }
    else if ( species == "BE" ){
         return 4;
    }
    else if ( species == "B" ){
         return 5;
    }
    else if ( species == "F" ){
         return 9;
    }
    else if ( species == "NE" ){
         return 10;
    }
    else if ( species == "MG" ){
         return 12;
    }    
    else if ( species == "AL" ){
         return 13;
    }    
    else if ( species == "SI" ){
         return 14;
    }
    else if ( species == "AR" ){
         return 18;
    }
    else if ( species == "K" ){
         return 19;
    }
    else if ( species == "CA" ){
         return 20;
    }
    else if ( species == "SC" ){
         return 21;
    }
    else if ( species == "TI" ){
         return 22;
    }
    else if ( species == "V" ){
         return 23;
    }
    else if ( species == "CR" ){
         return 24;
    }
    else if ( species == "MN" ){
         return 25;
    }
    else if ( species == "FE" ){
         return 26;
    }
    else if ( species == "CO" ){
         return 27;
    }
    else if ( species == "NI" ){
         return 28;
    }
    else if ( species == "CU" ){
         return 29;
    }
    else if ( species == "ZN" ){
         return 30;
    }
    else if ( species == "GA" ){
         return 31;
    }
    else if ( species == "GE" ){
         return 32;
    }
    else if ( species == "AS" ){
         return 33;
    }
    else if ( species == "SE" ){
         return 34;
    }
    else if ( species == "BR" ){
         return 35;
    }
    else if ( species == "KR" ){
         return 36;
    }
    else if ( species == "RB" ){
         return 37;
    }
    else if ( species == "SR" ){
         return 38;
    }
    else if ( species == "Y" ){
         return 39;
    }
    else if ( species == "ZR" ){
         return 40;
    }
    else if ( species == "NB" ){
         return 41;
    }
    else if ( species == "MO" ){
         return 42;
    }
    else if ( species == "TC" ){
         return 43;
    }
    else if ( species == "RU" ){
         return 44;
    }
    else if ( species == "RH" ){
         return 45;
    }
    else if ( species == "PD" ){
         return 46;
    }
    else if ( species == "AG" ){
         return 47;
    }
    else if ( species == "CD" ){
         return 48;
    }
    else if ( species == "IN" ){
         return 49;
    }
    else if ( species == "SN" ){
         return 50;
    }
    else if ( species == "SB" ){
         return 51;
    }
    else if ( species == "TE" ){
         return 52;
    }
    else if ( species == "I" ){
         return 53;
    }
    else if ( species == "XE" ){
         return 54;
    }
    else if ( species == "CS" ){
         return 55;
    }
    else if ( species == "BA" ){
         return 56;
    }
    else if ( species == "LA" ){
         return 57;
    }
    else if ( species == "CE" ){
         return 58;
    }
    else if ( species == "PR" ){
         return 59;
    }
    else if ( species == "ND" ){
         return 60;
    }
    else if ( species == "PM" ){
         return 61;
    }
    else if ( species == "SM" ){
         return 62;
    }    
    else if ( species == "EU" ){
         return 63;
    }
    else if ( species == "GD" ){
         return 64;
    }
    else if ( species == "TB" ){
         return 65;
    }
    else if ( species == "DY" ){
         return 66;
    }
    else if ( species == "HO" ){
         return 67;
    }
    else if ( species == "ER" ){
         return 68;
    }
    else if ( species == "TM" ){
         return 69;
    }
    else if ( species == "YB" ){
         return 70;
    }
    else if ( species == "LU" ){
         return 71;
    }
    else if ( species == "HF" ){
         return 72;
    }
    else if ( species == "TA" ){
         return 73;
    }
    else if ( species == "W" ){
         return 74;
    }
    else if ( species == "RE" ){
         return 75;
    }
    else if ( species == "OS" ){
         return 76;
    }
    else if ( species == "IR" ){
         return 77;
    }
    else if ( species == "PT" ){
         return 78;
    }
    else if ( species == "AU" ){
         return 79;
    }
    else if ( species == "HG" ){
         return 80;
    }
    else if ( species == "TL" ){
         return 81;
    }
    else if ( species == "PB" ){
         return 82;
    }
    else if ( species == "BI" ){
         return 83;
    }
    else if ( species == "PO" ){
         return 84;
    }
    else if ( species == "AT" ){
         return 85;
    } 
    else if ( species == "RN" ){
         return 86;
    }
    else if ( species == "FR" ){
         return 87;
    }
    else if ( species == "RA" ){
         return 88;
    }           
    else if ( species == "AC" ){
         return 89;
    }
    else if ( species == "TH" ){
         return 90;
    }
    else if ( species == "PA" ){
         return 91;
    }
    else if ( species == "U" ){
         return 92;
    }
    else{
         return 0; //report problem case!
    }
}

bool PDBResidue::isHistidine(){
     //
     if ( resname == "HIS" ) return true;
     if ( resname == "His" ) return true;
     if ( resname == "his" ) return true;
     return false;
     
}

bool isStandardResidue( string resname){
     to_upper(resname); //just to make sure!
     if ( resname == "ALA" || resname == "ARG" || resname == "ASN" ){
          return true;
     }
     else if ( resname == "ASP" || resname == "CYS" || resname == "LEU" ){
          return true;
     }
     else if ( resname == "LYS" || resname == "MET" || resname == "PHE" ){
          return true;
     }
     else if ( resname == "PRO" || resname == "GLN" || resname == "GLU" ){
          return true;
     }
     else if ( resname == "GLY" || resname == "HIS" || resname == "ILE" ){
          return true;
     }
     else if ( resname == "SER" || resname == "THR" || resname == "TRP" ){
          return true;
     }
     else if ( resname == "TYR" || resname == "VAL" || resname == "PYL" ){
          return true;
     }
     else if ( resname == "SEC" || resname == "ASX" || resname == "GLX" ){
          return true;
     }

     return false;
}

void PDBThing:: setData(){
     //turn PDB entry into logic
     //first use element setting...
     for ( int i = 0; i < res.size(); i++ ){
         //working on the ith residue
         res.at(i).isStandard = isStandardResidue( res.at(i).resname );
         
         for ( int j = 0; j < res.at(i).atom.size(); j++ ){
             int k = res.at(i).atom.at(j);
             //working on kth atom of PDBThing
             //jth atom of this residue
             
             atom.at(k).isStandard = res.at(i).isStandard; //inherits
             atom.at(k).el = species2number( atom.at(k).element );
             string mytype = atom.at(k).recname;
             to_upper( mytype );
             if ( mytype == "ATOM  " ){
                  atom.at(k).isAtom = true;
                  atom.at(k).isHet = false;
             }
             else {
                  atom.at(k).isAtom = false;
                  atom.at(k).isHet = true;
             }
             
             bool foundCAYet = false;
             to_upper( atom.at(k).name ); //just in case
             if ( atom.at(k).name == " CA " ){
                  if (foundCAYet ){
                     continue; // this is the pyrrolysine trap for sidechain CA
                  }
                  else{
                       foundCAYet = true;
                       atom.at(k).isMain = true;
                       atom.at(k).isAlpha = true;
                       //cerr << "Setting isAlpha on atom " << k+1 << ": " << atom.at(k).atom2line() << endl;
                  }
             }
             else if ( atom.at(k).name == " N  " ){
                  atom.at(k).isMain = true;
             }
             else if ( atom.at(k).name == " C  " ){
                  atom.at(k).isMain = true;
             }
             else if ( atom.at(k).name == " H  " ){
                  atom.at(k).isMain = true;
             }                          
             else if ( atom.at(k).name == " O  " ){
                  atom.at(k).isMain = true;
             } //and don't forget the weird terminal oxygen cases this time     
             else if ( atom.at(k).name == " OXT" ){
                  atom.at(k).isMain = true;
             }
             else if ( atom.at(k).name == " OT " ){
                  atom.at(k).isMain = true;
             }
             else if ( atom.at(k).name == " OT1" ){
                  atom.at(k).isMain = true;
             }
             else if ( atom.at(k).name == " OT2" ){
                  atom.at(k).isMain = true;
             }
 
             if ( atom.at(k).elO() || atom.at(k).elN() ) atom.at(k).isPolar = true;
             if ( atom.at(k).elS() || atom.at(k).elSe() ) atom.at(k).isPolar = true;
             if (atom.at(k).elP() ) atom.at(k).isPolar = true;
             //later: label polar hydrogens in bonding
             
         }
     }
     
     
}

/*
  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
*/
void PDBThing::initialiseGrid(){
     cerr << "Initialising coarse grid array." << endl;
     //first set corners
     if ( atom.size() < 1 ) return; //no atoms!
     cerr << "Gridding " << atom.size() << " atoms." << endl;
     lowcorner = atom.at(0).pos;
     highcorner = atom.at(0).pos;
     for ( int i = 0 ; i < atom.size() ; i++ ){
         Vector tempos = atom.at(i).pos;
         if ( tempos.x < lowcorner.x ) lowcorner.x = tempos.x;
         if ( tempos.y < lowcorner.y ) lowcorner.y = tempos.y;         
         if ( tempos.z < lowcorner.z ) lowcorner.z = tempos.z;
         if ( tempos.x > highcorner.x ) highcorner.x = tempos.x;
         if ( tempos.y > highcorner.y ) highcorner.y = tempos.y;         
         if ( tempos.z > highcorner.z ) highcorner.z = tempos.z;
     }
     cerr << "Low corner: " << lowcorner << endl;
     cerr << "High corner: " << highcorner << endl;
     //padding
     Vector padder = Vector( gridpad, gridpad, gridpad );
     lowcorner -= padder;
     highcorner += padder;
     //padded!
     cerr << "After padding:" << endl;
     cerr << "Low corner: " << lowcorner << endl;
     cerr << "High corner: " << highcorner << endl;     
     //how many cells:
     Vector extent = highcorner - lowcorner;
     Nx = 1+floor( extent.x / gridmin );
     Ny = 1+floor( extent.y / gridmin );     
     Nz = 1+floor( extent.z / gridmin );
     cerr << "Grid extent: X " << Nx << " " << Ny << " " << Nz << endl;     
     //Now size the grid appropriately
     grid.resize( Nx * Ny * Nz );
     cerr << "Grid size: " << grid.size() << endl;
     //and we're done
     return;
}


void PDBThing::fillGrid(){
     //run on clean grid after initialise
     for ( int i = 0 ; i < atom.size(); i++ ){
         int c = gridBox( atom.at(i) );
         grid.at(c).push_back( i );
     }
     return; // done!
}
int PDBThing::gridX( Vector& pos ){
    int r = floor( (pos.x - lowcorner.x ) / gridmin );
    if ( r < 0 ){
         cerr << "FATAL ERROR: grid cell < 0 (" << r <<")" << endl;
         exit(1);
    }
    if ( r > Nx-1 ){
         cerr << "FATAL ERROR: grid cell > Nx-1 (" << r <<")" << endl;
         exit(1);
    }
    return r;
}
int PDBThing::gridY( Vector& pos ){
    int r = floor( (pos.y - lowcorner.y ) / gridmin );
    if ( r < 0 ){
         cerr << "FATAL ERROR: grid cell < 0 (" << r <<")" << endl;
         exit(1);
    }
    if ( r > Ny-1 ){
         cerr << "FATAL ERROR: grid cell > Ny-1 (" << r <<")" << endl;
         exit(1);
    }
    return r;
}
int PDBThing::gridZ( Vector& pos ){
    int r = floor( (pos.z - lowcorner.z ) / gridmin );
    if ( r < 0 ){
         cerr << "FATAL ERROR: grid cell < 0 (" << r <<")" << endl;
         exit(1);
    }
    if ( r > Nz-1 ){
         cerr << "FATAL ERROR: grid cell > Nz-1 (" << r <<")" << endl;
         exit(1);
    }
    return r;
}
int PDBThing::gridBox( Vector &pos ){
    int x = gridX( pos );
    int y = gridY( pos );
    int z = gridZ( pos );
    int r = z + y*Nz + x*Ny*Nz;
    return r;
}

int PDBThing::gridBox ( PDBAtom &atom ){
    return gridBox( atom.pos );
}

vector< int > PDBThing::near ( PDBAtom &a, double within ){
        Vector apos = a.pos;
        Vector opos;
        int ax, ay, az, ox, oy, oz;
        int acell, ocell;
        int oat;
        
        ax = gridX( apos );
        ay = gridY( apos );
        az = gridZ( apos );
        acell = gridBox( apos );
        vector< int > result;
        double w = within*within;
        //double nwithin = -1.*within;
        
        for ( ox =ax -1 ; ox < ax+2 ; ox++ ){
            if ( ox < 0 ) continue; //don't search out of grid
            if ( ox > Nx -1 ) continue;
            
            for ( oy = ay-1; oy < ay+2; oy++ ){
                if ( oy < 0 ) continue;
                if ( oy > Ny-1 ) continue;
                
                for ( oz = az-1; oz < az+2; oz++ ){
                    if ( oz < 0 ) continue;
                    if ( oz > Nz-1 ) continue;
                    
                    ocell = oz + oy*Nz + ox*Ny*Nz; // other cell!
                    for ( int i = 0; i < grid.at(ocell).size(); i++ ){
                        oat = grid.at(ocell).at(i);
                        if ( oat == a.myID ) continue; //skip self
                        opos = relvec( a, atom.at(oat) );
                        //if ( opos.x > within || opos.y > within || opos.z > within ) continue; //test before squaring for speed?
                        //if ( opos.x < nwithin || opos.y < nwithin || opos.z < nwithin ) continue;
                        //but six comparison are probably more expensive than one mult and a comparison! Don't bother
                        double dist = opos.sq();
                        if ( dist > w ) continue;
                        result.push_back( oat );
                    }
                    
                }
            }
            
        }
        return result; //array of nearbys
}
vector< int > PDBThing::near ( Vector& pos, double within ){
        Vector opos;
        int ax, ay, az, ox, oy, oz;
        int acell, ocell;
        int oat;
        
        ax = gridX( pos );
        ay = gridY( pos );
        az = gridZ( pos );
        acell = gridBox( pos );
        vector< int > result;
        double w = within*within;
        
        for ( ox =ax -1 ; ox < ax+2 ; ox++ ){
            if ( ox < 0 ) continue; //don't search out of grid
            if ( ox > Nx -1 ) continue;
            
            for ( oy = ay-1; oy < ay+2; oy++ ){
                if ( oy < 0 ) continue;
                if ( oy > Ny-1 ) continue;
                
                for ( oz = az-1; oz < az+2; oz++ ){
                    if ( oz < 0 ) continue;
                    if ( oz > Nz-1 ) continue;
                    
                    ocell = oz + oy*Nz + ox*Ny*Nz; // other cell!
                    for ( int i = 0; i < grid.at(ocell).size(); i++ ){
                        oat = grid.at(ocell).at(i);
                        opos = atom.at(oat).pos - pos; //relative to point...
                        double dist = opos.sq();
                        if ( dist > w ) continue;
                        result.push_back( oat );
                    }
                    
                }
            }
            
        }
        return result; //array of nearbys
}

Vector PDBThing::rel( int i, int j ){
       return relvec( atom.at(i), atom.at(j) );
}

double PDBThing::d( int i, int j ){
       return drel( atom.at(i), atom.at(j) );
}

double PDBThing::d2( int i, int j ){
       Vector v;
       v = rel( i, j );
       return v.sq(); //squared value
}

void PDBThing::assignPhobicRadii(){
	radius.resize( atom.size() ); //refresh and size the radius array
	//assign radii to the main phobic atom cases
	for ( int i = 0 ; i < atom.size(); i++ ){
		if ( atom.at(i).elC() ){
			radius.at(i) = 1.7;
		}
		else if ( atom.at(i).elS() ){
			radius.at(i) = 1.8;
		}
	}
	
}

void PDBThing::assignContactRadii(){
	radius.resize( atom.size() ); //refresh and size the radius array
	//assign radii to the recognisable elements
	for ( int i = 0 ; i < atom.size(); i++ ){
		if ( atom.at(i).elC() ){
			radius.at(i) = 1.7;
		}
		else if ( atom.at(i).elS() ){
			radius.at(i) = 1.8;
		}
		else if ( atom.at(i).elN() ){
			radius.at(i) = 1.5; //discuss...
		}
		else if ( atom.at(i).elO() ){
			radius.at(i) = 1.6; //discuss...
		}
		else if ( atom.at(i).elH() ){
			radius.at(i) = 1.0; //discuss...
		}
		else if ( atom.at(i).elP() ){
			radius.at(i) = 1.8; //discuss...
		}
		else {
			radius.at(i) = 1.5; //discuss...
		}
	}
	
}


void PDBThing::assignBurialRadii(){
     //give the radii used in the "sexton" routine
     //only heavy atoms are calculated for
     //radii for elements;
     //generic metal radius?
     
     int nat = atom.size(); //number of atoms in this Thing
     radius.resize( nat );
     for ( int i = 0 ; i < nat ; i++ ){
         if ( atom.at(i).el == 1 ){
              //hydrogen case
              //fictitious radius as not used
              radius.at(i) = 1.0;
         }
         else if ( atom.at(i).el == 6 ){
              //carbon
              //specialcase the mainchain carboxyl
              if ( atom.at(i).isMain ){
                   radius.at(i) = 1.5;
              }
              else{
                   radius.at(i) = 2.0; //includes some hydrogen cover!
              }
         }
         else if ( atom.at(i).el == 7 ){
              radius.at(i) == 1.5; //nitrogen case
         }
         else if ( atom.at(i).el == 8 ){
              radius.at(i) == 1.5; //nitrogen case
         }
         else if ( atom.at(i).el == 16 ){
              radius.at(i) == 1.9; //sulphur case
         }
         else radius.at(i) == 1.5; //generic case, rearrange later?
     }
}

void PDBThing::sexton(double prober){
     //note: usual prober is 1.4 for water!
     cerr << "Carrying out SEXTON routine for burial calculation." << endl;
     cerr << "Calculates burial distance for heavy atoms." << endl;
     cerr << "Burial distance written into occupancy column as dummy." << endl;
     cerr << "Using probe radius: " << prober << endl;
     
     assignBurialRadii(); // give radii to atoms, crude heavy-only version
     
     //loop over atoms;
     
     int nat = atom.size();
     vector< int > surfaceAtom;
     vector< int > buriedAtom;
     
     for ( int i = 0 ; i < nat ; i++ ){
         if ( atom.at(i).el == 1 ){
              atom.at(i).occ = 99.99; //fictitious suspicious value!
              continue; // skipping hydrogens!
         }
         bool exposed = false;
         //candidate collision atoms are those near the current atom under test;
         vector< int > nearby;
         double wit = 3.0 + prober + radius.at(i);
         nearby = near( atom.at(i), wit ); 
         //build solvent points in a spiral around the atom
         int Npoints = 400;
         double rad = prober + radius.at(i);
         double inc = mypi * (3.0 - sqrt(5.0)); //golden section
         double y,r,phi;
         for ( int k = 0; k < Npoints; k++){
             bool blocked = false;// is this point blocked by another atom?
             y = -1.0 + ( 2.0*k/(Npoints-1) ); // -1 to -1
             r = sqrt( 1.0 - y*y ); //r,y make a unit vector...
             phi = k*inc; //angle increment
             Vector dum = Vector( r*cos(phi), y, r*sin(phi) );
             dum *= rad; //point now represents a solvent probe touching this atom
             Vector pos = atom.at(i).pos + dum;
             //test against neighbours of atom i
             int nn = nearby.size();
             for ( int o = 0; o < nn; o++ ){
                 int oat = nearby.at(o); //other atom
                 if ( atom.at(oat).el == 1 ) continue; //skip hydrogens again
                 Vector opos = atom.at(oat).pos;
                 Vector dr = opos - pos; //relative;
                 double crit = prober + radius.at(oat); //criterion
                 double crit2 = crit*crit;
                 double dr2 = dr.sq(); // size of relative
                 if ( dr2 < crit2 ) blocked = true;
                 if ( blocked ) break; //leave the nearby loop; no need to check other atoms
             }
             if (blocked ){
                         continue; //check another point
             } 
             else{
                  exposed = true; // this atom can touch solvent!
             }
             if ( exposed ) break; // no need to check more solvent points
         }       
         if ( exposed ){
              surfaceAtom.push_back( i ) ; //goes into surface collection
              atom.at(i).occ = 0.01; //"exposed" value
         }
         else{
              buriedAtom.push_back( i );
         }    
     }
     
     //now loop over buried atoms checking for the closest surface atom
     cerr << "Found " << surfaceAtom.size() << " surface atoms." << endl;
     cerr << "Found " << buriedAtom.size() << " buried atoms." << endl;
     
     for ( int which = 0 ; which < buriedAtom.size(); which++ ){
         int bat = buriedAtom.at(which);
         double mind2 = 1e10; // ridiculous length
         for ( int other = 0 ; other < surfaceAtom.size(); other++ ){
             int sat = surfaceAtom.at(other);
             Vector v = rel( bat, sat); //relative vector;
             double d2 = v.sq();
             if ( d2 < mind2 ) mind2 = d2; // keep shortest bury distance
         }
         double deep = sqrt( mind2 );
         atom.at(bat).occ = deep; //value goes into occ      
     }
     
}

bool PDBThing::writeMap(string filename){
     cerr << "Writing MAP." << endl;
     //open the file
     ofstream outf;
     outf.open( filename.c_str(), ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file " << filename << endl;
          return false;
     }
     
     outf << "FILES" << endl;    
     //write the list of files in it iwth model numbers
     for ( int f = 0 ; f < file.size(); f++ ){
         outf << f << " " << file.at(f).name;
         if ( file.at(f).hasModels ){
              outf << " " << file.at(f).model.size();
         }
         outf << endl;
     }
     
     outf << "MAP" << endl;
     
     //for each serial
     //write index; file-model-serial; human-readable id
     for ( int i = 0 ; i < atom.size(); i++ ){
         outf << i+1 << " ";
         stringstream foo;
         int f = atom.at(i).file;
         if ( file.at(f).hasModels ){
              int m = atom.at(i).model;
              int p = model.at(m).modelnum;
              foo << f << "-" << p << "-" << atom.at(i).serial;
         }
         else{
              foo << f << "-" << atom.at(i).model << "-" << atom.at(i).serial;
         }
         
         outf << foo.str() << " ";
         outf << atom.at(i).atom2id() << endl;
     }
     outf.close(); 
     return true;
}



Vector relvec( PDBAtom &a, PDBAtom &b ){
       Vector result = b.pos - a.pos; // convention is From -> To = pos(To)-pos(From).
       return result;
}

double drel( PDBAtom &a, PDBAtom &b ){
       Vector r = relvec( a, b );
       return sqrt( r.sq() );
}

bool isSameResidue( PDBAtom &a, PDBAtom &b ){
	//use the internal "res" index for simplicity
     if ( a.res != b.res ) return false;
     return true;
}

bool isResidueWithin( PDBAtom &a, PDBAtom &b, int c ){
	if ( a.chain != b.chain ) return false; // not same chain, drop it
	int d = abs( a.res - b.res ); //absolute difference
	if ( d > c ) return false; //this will also return false if c < 0, handily
	return true;
} //same chain, residue difference up to c

bool isResidueAt( PDBAtom &a, PDBAtom &b, int c ){
	if ( a.chain != b.chain ) return false; // not same chain, drop it
	int d = abs( a.res - b.res ); //absolute difference
	if ( d != c ) return false; //this will also return false if c < 0, handily
	return true;	
}//same chain, residue difference exactly c
//note, isResidueAt( a,b,1) is effectively "isAdjacentResidue( a,b )
//so I'll overload that :)
bool isAdjacentResidue( PDBAtom &a, PDBAtom &b ){
	return isResidueAt( a, b, 1 );
}


void PDBThing::clearCov(){
     cov.clear(); // not difficult :)
     int covsize = atom.size(); //number of atoms in this Thing
     cov.resize( covsize ); //is right size now
     return; // done
}

int PDBThing::IDfromSerial( int start, int N, int serial ){
    //start should be the beginning of a MODEL's atom list
    //N should be its size
    int result;
    result = -1;
    for ( int a = start ; a < start+N ; a++ ){
        int s = atom.at(a).serial;
        if ( s == serial ){
             result = a;
             return result;
        }
    }
    
    return result;
}

void PDBThing::conectToCov(){
     //loop over models;
     cerr << "Bonding phase 0: checking CONECT." << endl;
     int msize = model.size();
     for ( int m = 0 ; m < msize; m++ ){
         int consize = model.at(m).conect.size();
         if ( consize == 0 ) continue; // just for clarity; we'll move on
         
         for ( int c = 0 ; c < consize ; c++ ){
             int atom1, atom2, atomA, atomB;
             atom1 = model.at(m).conect.at(c).atom1; //this is a serial number
             cerr << "CONECT serial " << atom1 << " ";
             int start = model.at(m).atom.at(0); //first atom in model
             int N = model.at(m).atom.size(); //number of atoms in model
             
             atomA = IDfromSerial( start, N, atom1 );
             cerr << " ID " << atomA+1 << "  ";
             
             if (atomA == -1 ){
                       //this is not good
                       cerr << endl;
                       continue; //let's pass on
             }
             //okay now let's do the next atom
             
             int nb = model.at(m).conect.at(c).bond.size();
             for ( int b = 0 ; b < nb ; b++ ){
                 int atom2 = model.at(m).conect.at(c).bond.at(b); //other atom serial
                 atomB = IDfromSerial( start, N, atom2 );
                 cerr << " other serial " << atom2 << " ";
                 cerr << " ID " << atomB+1 << endl;
                 if ( atomB == -1 ){
                      continue; //not good, skip
                 }
                 
                 //okay so now we have atomA and atomB
                 //put them into each other's cov arrays.
                 cerr << "Covalent connection from CONECT: atom IDs ";
                 cerr << atomA+1 << " " << atomB+1 << endl;
                 insertIfNotIn( atomA, cov.at(atomB) );
                 insertIfNotIn( atomB, cov.at(atomA) );
             }
             
                       
         }
         
     }
     
}

void PDBThing::getClosebyForCov(){
     int nat = atom.size();
     //first step: get a nearby list for each atom
     //use different cuts for each type
     cerr << "Getting CLOSEBY list before checking bonding." << endl;
     closeby.resize( nat );
     for ( int i = 0 ; i < nat; i++ ){
         double l = 2.5; // look within a long bond range for starters, filter later
         closeby.at(i) = near( atom.at(i), l );

/*
         for ( int j = 0 ; j < closeby.at(i).size() ; j ++ ){
             int k = closeby.at(i).at(j);
             double r = d( i, k );
             //cerr << j << " " << setprecision(3) << r << " " ;
             //cerr  << k << " " << atom.at(k).atom2line() << endl;
         }
*/
     }
     
}
     //okay now what do I DO with this info?
     //CONECT intervention:
     //if there's CONECT info
     
     //if an atom is in a CONECT and it's a HET,
     // then the CONECT replaces all previous candidates
     //but if an atom is in a CONECT and it's an ATOM,
     // then the CONECT is ADDED to previous candidates
     //note that the CONECT is a compulsory addition not a suggestion.
     
//completeCov is to be invoked after conectToCov
//for HET, if it already has cov, don't do further bonding because CONECT has already provided it
//but if no cov, bond as normal
//for ATOM, do normal bonding as well as CONECT bonding
void PDBThing::completeCov(){
     int nat = atom.size(); //numberof atoms in list; we'll cycle repeatedly
     //set up the sp2 array for later
     sp2.clear();
     for ( int i = 0 ; i < nat ; i++ ){
         sp2.push_back( false );
     }
     
     //first pass: get bonding for hydrogens. Needs only one each
     //must be same residue!
     //consider bonding to CONP; S; Se; "other"
     //bond limits: CONP 0.9-1.2, S 1.1-1.4, Se 1.4-1.7, "other" 0.8-1.5
     //function: bondH
     
     cerr << "Bonding phase 1: checking HYDROGENs." << endl;
     for ( int i = 0 ; i < nat ; i++ ){
         //working on atom index i
         if ( !atom.at(i).elH() ) continue; //skip everything that is not H
         int nc = cov.at(i).size(); //how many cov so far? in case conect
         if ( nc > 0 ) continue; //skip bonded H
         int nn = closeby.at(i).size(); // how many within 2.5A?
         for ( int j = 0 ; j < nn ; j++ ){
             int o = closeby.at(i).at(j); // index of "other" atom
             if ( atom.at(o).elH() ) continue; //don't bond H to H
             if ( !isSameResidue( atom.at(i), atom.at(o) ) ) continue; //skip not same residue
             double d = drel( atom.at(i), atom.at(o) ); //distance...
             if ( atom.at(o).elCONP() ){
                  if ( d > 0.8 && d < 1.2 ){
                       //good case: same residue CONP in right range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       if ( atom.at(o).isPolar ) atom.at(i).isPolar = true; // polar H on polar O or N or S (or P!)
                       //because we are hydrogen, one good case is enough
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;                      
                       break; // leave the loop
                  } 
             } //end isCONP case
             else if ( atom.at(o).elS() ){
                  if ( d > 1.1 && d < 1.4 ){
                       //good case: same residue S in right range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //because we are hydrogen, one good case is enough
                       atom.at(i).isPolar = true;
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       break; // leave the loop
                  } 
             } //end isS case
             else if ( atom.at(o).elSe() ){
                  if ( d > 1.3 && d < 1.7 ){
                       //good case: same residue Se in right range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //because we are hydrogen, one good case is enough
                       atom.at(i).isPolar = true;
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       break; // leave the loop
                  } 
             } //end isSe case
             else{
                  //odd case? check for problems...
                  if ( d > 0.8 && d < 1.5 ){
                       //good case: same residue other atom in rational range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //because we are hydrogen, one good case is enough
                       if ( atom.at(o).isPolar ) atom.at(i).isPolar = true;
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       break; // leave the loop
                  }
             } //done
         }
         
         if ( cov.at(i).size() == 0 ){
              cerr << "WARNING: unbonded hydrogen, index " << i << endl;
              cerr << atom.at(i).atom2line() << endl; //for debugging
         }
             
     }
     
     //second pass: examine conventional atoms, CONP
     //consider bonding to CONP; S; Se; "other" (?metal?)
     //bond limits: CONP 1.2-1.7, S 1.7-2.0, Se 1.7-2.2, "other" variable...
     //watch particularly for polar (O,N) coordination of metals at 2-2.5 Angstroms
     //if isStandard, should bond only to:
     // ATOMs in SAME residue, not other residue;
     //HETATMS;
     // if backbone, to backbone ATOMs in previous or next residue
     //function: bondStandardCONP
     //if not Standard, or HET: bond generally
     //function: bondGeneralCONP
     cerr << "Bonding phase 2: checking CONP." << endl;
     for ( int i = 0 ; i < nat ; i++ ){
         if ( !atom.at(i).elCONP() ) continue; //skip not CONP
         //remember hydrogens already done!
         int maxbond;
         if ( atom.at(i).elO() ){
              maxbond = 2; //oxygen case
         } 
         else if ( atom.at(i).elN() ){
              maxbond = 4; //nitrogen... up to 4? but allow 2 or 3
         }
         else {
              maxbond = 4; //carbon or P
         }
         if ( cov.at(i).size() >= maxbond ) continue; //skip fully bonded
         //this should trap HETs where COV has already done it for us

         bool sameOnly = true;
         if ( atom.at(i).isHet ) sameOnly = false; //generalise bonding for HET if no CONECT
         if ( atom.at(i).isMain && atom.at(i).elC() && !atom.at(i).isAlpha ) sameOnly = false; // mainchain C
         if ( atom.at(i).isMain && atom.at(i).elN() ) sameOnly = false; // mainchain N.
         int nn = closeby.at(i).size(); // how many within 2.5A?
         for ( int j = 0 ; j < nn ; j++ ){
             bool gechecked = false;
             int o = closeby.at(i).at(j); // index of "other" atom
             if ( atom.at(o).elH() ) continue; // done hydrogens already!
             if ( sameOnly && !isSameResidue( atom.at(i), atom.at(o) ) ) continue; //skip not same residue
             if ( cov.at(i).size() >= maxbond ) continue; //already filled?
             double d = drel( atom.at(i), atom.at(o) ); //distance...
             if ( atom.at(o).elCONP() ){
                  gechecked = true;
                  if ( d > 1.0 && d < 1.8 ){
                       //good case: same residue CONP in right range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isCONP case
             else if ( atom.at(o).elS() ){
                  gechecked = true;
                  if ( d > 1.4 && d < 2.0 ){
                       //good case: same residue S in right range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isS case             
             else if ( atom.at(o).elSe() ){
                  gechecked = true;
                  if ( d > 1.5 && d < 2.2 ){
                       //good case: same residue Se in right range
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isSe case             
             //further check for polar groups coordinating metals?
             if ( gechecked ) continue;
             if ( !atom.at(i).isPolar ) continue; //don't further check nonpolar
             //if here we are a polar atom talking to something not HCONPSSe
             if ( d > 1.2 && d < 2.4 ){
                  insertIfNotIn( i, cov.at(o) );
                  insertIfNotIn( o, cov.at(i) );
                  //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                  continue;
             } 
             
         }
         
     }
     
     //third pass: examine S
     //consider bonding to S; Se; "other"
     //disulfide: 2 to 2.4
     //to selenium ? 2 to 2.5?
     //"other": particularly metals around 2 to 2.5
     //function: bondS
     cerr << "Bonding phase 3: checking SULPHUR." << endl;
     for ( int i = 0 ; i < nat ; i++ ){
         if ( !atom.at(i).elS() ) continue; // only S this time
         int maxbond = 4; // for now
         if ( cov.at(i).size() >= maxbond ) continue; //skip fully bonded
         int nn = closeby.at(i).size(); // how many within 2.5A?
         for ( int j = 0 ; j < nn ; j++ ){
             if ( cov.at(i).size() >= maxbond ) continue; //skip fully bonded
             int o = closeby.at(i).at(j); // index of "other" atom
             if ( atom.at(o).elH() ) continue; // done hydrogens already!
             if ( atom.at(o).elCONP() ) continue; // done CONP already!
             double d = drel( atom.at(i), atom.at(o) ); //distance...
             if ( atom.at(o).elS() ){
                  if ( d > 1.8 && d < 2.3 ){
                       //good case: probably a disulphide!
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isS case             
             else if ( atom.at(o).elSe() ){
                  if ( d > 1.9 && d < 2.4 ){
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isSe case
             //What about sulphur coordinating metals?
             //check on elTrans() for transition metals
             else if ( atom.at(o).elTrans() ){
                  if ( d > 1.9 && d < 2.4 ){
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isTransitionMetal case             
             
         }
         
     }
     
     //fourth pass: examine Se (if any)
     //consider bonding to Se; "other"
     //Se-Se is around 2.1-2.5
     //"other": particularly metals around 2 to 2.5
     //function: bondSe
     cerr << "Bonding phase 4: checking SELENIUM." << endl;
     for ( int i = 0 ; i < nat ; i++ ){
         if ( !atom.at(i).elSe() ) continue; // only S this time
         int maxbond = 4; // for now
         if ( cov.at(i).size() >= maxbond ) continue; //skip fully bonded
         int nn = closeby.at(i).size(); // how many within 2.5A?
         for ( int j = 0 ; j < nn ; j++ ){
             if ( cov.at(i).size() >= maxbond ) continue; //skip fully bonded
             int o = closeby.at(i).at(j); // index of "other" atom
             if ( atom.at(o).elH() ) continue; // done hydrogens already!
             if ( atom.at(o).elCONP() ) continue; // done CONP already!
             if ( atom.at(o).elS() ) continue; // done S already!             
             double d = drel( atom.at(i), atom.at(o) ); //distance...
             if ( atom.at(o).elSe() ){
                  if ( d > 1.9 && d < 2.4 ){
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isSe case
             //check on elTrans() for transition metals
             else if ( atom.at(o).elTrans() ){
                  if ( d > 1.9 && d < 2.4 ){
                       insertIfNotIn( i, cov.at(o) );
                       insertIfNotIn( o, cov.at(i) );
                       //cerr << "Bonding index " << i << " to " << o << " at " << d << endl;
                       continue;
                  } 
             } //end isTransitionMetal case             
             
         }
         
     }
     
     //populate the nH property by counting hydrogen neighbours of heavy atoms
     
     for ( int i = 0 ; i < atom.size() ; i++ ){
		 int myNH = 0;
		 for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			 int k = cov.at(i).at(j); //bonded atom;
			 if ( atom.at(k).elH() ) myNH++;
		 }
		 atom.at(i).nH = myNH; //inform the troops
	 }
     
     sp2Check();
     barCheck();
     
     //Problem reporting: report uncoordinated atoms and under-coordinated carbons
     for ( int i = 0; i < nat ; i++ ){
         bool unco = false;
         bool underco = false;
         if ( cov.at(i).size() == 0 ) unco = true;
         if ( atom.at(i).elC() && ( cov.at(i).size() < 3 ) ) underco = true;
         
         if ( unco ){
              cerr << "UNCOORDINATED atom, index " << i+1 << " ";
              cerr << atom.at(i).atom2id() << endl;
         }
         else if ( underco ){
              cerr << "UNDERCOORDINATED atom, index " << i+1 << " ";
              cerr << atom.at(i).atom2id() << endl;
         }
         
         //for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			// cerr << fullCovalentLine( i, j ) << endl;
		 //}
     }

	
    
     
}

string PDBThing::fullCovalentLine( int i, int j ){
	stringstream record;
	if ( i > cov.size() - 1 || i < 0 ) return record.str(); // safety!
	if ( j < 0 || j > cov.at(i).size() - 1 ) return record.str();
	int o = cov.at(i).at(j);
	record << i+1 << " " << o+1;
	record << " " << bars.at( i ).at( j ) << " :COVALENT: ";
	record << atom.at( i ).atom2id();
	record << " :: " << atom.at( o ).atom2id();
	record << " :d: " << d( i, o );

	return record.str();
}

bool PDBThing::writeCov( string filename ){
     cerr << "Writing COV to file " << filename << endl;
     //open the file
     ofstream outf;
     
     outf.open( filename.c_str(), ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file " << filename << endl;
          return false;
     }
     
     for ( int i = 0 ; i < cov.size() ; i++ ){
		 for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			 int o = cov.at(i).at(j);
			 if ( o < i ) continue; //skipping dupes
			 outf << fullCovalentLine( i, j ) << endl;
		 }
     }
     
     outf.close();
     
     outf.open("cov.in", ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file cov.in" << endl;
          return false;
     }
     
     for ( int i = 0 ; i < cov.size() ; i++ ){
		 for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			 int o = cov.at(i).at(j);
			 if ( o < i ) continue; //skipping dupes
			 outf << atom.at( i ).serial << " " << atom.at( o ).serial << " " << bars.at( i ).at( j ) << endl;
		 }
	 }
     
     outf.close();
     
     return true;
}

bool PDBThing::writeHydrogenBonds( string filename ){
     cerr << "Writing HB to file " << filename << endl;
     //open the file
     ofstream outf;
     
     outf.open( filename.c_str(), ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file " << filename << endl;
          return false;
     }
     
     for ( int i = 0 ; i < hb.size() ; i++ ){
		 outf << fullHbondLine( hb.at(i) ) << endl;
     }
     
     outf.close();
     
     outf.open( "hbonds.in", ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file hbonds.in" << endl;
          return false;
     }
     
     for ( int i = 0 ; i < hb.size() ; i++ ){
		 outf << atom.at( hb.at( i ).h ).serial << " " << atom.at( hb.at( i ).acceptor ).serial << " " << hb.at( i ).E << " 5" << endl;
     }
     
     outf.close();
     
     return true;	
}

bool PDBThing::writePhobicTethers( string filename , int Cryo ){
     cerr << "Writing PH to file " << filename << endl;
     //open the file
     ofstream outf;
     
     outf.open( filename.c_str(), ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file " << filename << endl;
          return false;
     }
     
     for ( int i = 0 ; i < phobe.size() ; i++ ){
		 outf << fullPhobicPairLine( phobe.at(i) ) << endl;
     }
     
     outf.close();
     
     outf.open( "hphobes.in", ios::out );
     if ( outf.fail() ){
          cerr << "ERROR: cannot open file hphobes.in" << endl;
          return false;
     }
     
     for ( int i = 0 ; i < phobe.size() ; i++ ){
		 outf << atom.at( phobe.at( i ).a1 ).serial << " " << atom.at( phobe.at( i ).a2 ).serial << " ";
		 if(Cryo) { outf << "1" << endl; }
		 else { outf << "2" << endl; }
     }
     
     outf.close();
     
     return true;	
}


void PDBThing::sp2Check(){
     //scan over the cov array by atoms
     //skip all atoms except C,O,N
     //assign sp2 status
     //C,N labelled sp2 if less than 4 neighbours
     //O labelled sp2 if one neighbour
     cerr << "Checking sp2 status." << endl;
     sp2.clear();
     for ( int i = 0 ; i < cov.size() ; i++ ){
         if ( atom.at(i).isAlpha ){
              sp2.push_back( false );
         }
         else if ( atom.at(i).elO() && cov.at(i).size() == 1 ){
              sp2.push_back( true );
         }
         else if ( atom.at(i).elC() && cov.at(i).size() < 4 ){
              sp2.push_back( true );
         }
         else if ( atom.at(i).elN() && cov.at(i).size() < 4 ){
              sp2.push_back( true );
         }
         else{
              sp2.push_back ( false );
         }
     }
     
     for ( int i = 0 ; i < sp2.size() ; i++ ){
         if ( sp2.at(i) ){
              //cerr << "Identified sp2 atom, index " << i << " " << atom.at(i).atom2id() << endl;
              //also: remember to label polar C
              //if sp2 and has a neighbour that isPolar
              if ( !atom.at(i).elC() ) continue; //skip not carbon
              for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
                  int k = cov.at(i).at(j); //a bonded atom
                  if ( atom.at(k).elO() && sp2.at(k) ){
                       atom.at(i).isPolar = true; //polarise C bonded to carboxyl O
                       //cerr << "Index " << i << " is polar carbon, " << atom.at(i).atom2id() << endl ;
                       break; //one is enough
                  }
              }
              
         }
     }
     
}

void PDBThing::barCheck(){
     //roll over atoms
     //roll over cov
     //if both atoms are sp2
     //then bars = 6
     //or if one atom is a dangler
     //then bars = 6
     //this is to guarantee good cluster building
     //otherwise
     //bars = 5
     bars.clear();
     bars.resize( atom.size() );
     
     for ( int i = 0 ; i < cov.size() ; i++ ){
         for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
             int o = cov.at(i).at(j); //id of other atom
             if (  sp2.at(i) && sp2.at(o) ){
                   //is locked
                   bars.at(i).push_back( 6 );
                   //if ( i < o ){
                        //cerr << "Found LOCKED bond between " << i << " and " << o << endl;
                        //cerr << atom.at(i).atom2id() << "-------" << atom.at(o).atom2id() << endl;
                   //}
             }
             else if ( cov.at(i).size() == 1 || cov.at(o).size() == 1 ){
				 //dangling atom
				 bars.at(i).push_back( 6 );
			 }
             else{
                  bars.at(i).push_back( 5 );
             }
             
         }     
     }
}



void PDBThing::assignPolarityLabels(){
	//bool isDonor, isAcceptor, isChargedDonor, isChargedAcceptor;
	//recall isPolar was set on H atoms during cov bond testing phase...
	//and was set on N.O.S.Se.P during setData phase
	//set isDonor if I'm a polar atom with hydrogens to donate
	//set isAcceptor if I'm a polar atom that could accept hbonds
	//set Charged labels on nitrogen of +charged groups and oxygen of carboxylate-
	//I have special cases for guani (arginine) and imi (histidine) group nitrogens
	//this means a "true base" will be a non-hydrogen polar that is not an acceptor OR donor;
	//only carbonyl C should receive this label
	
	//hey, also assign isPhobic label to appropriate cases:
	//non-Main non-Polar carbon;
	//sulphur (note sulphur both isPhobic and isPolar!)
	//carbon and sulphur should not be Phobic if they have a strongly polar neighbour
	//i.e. if their neighbour is not H and not C and not S and isPolar, they are not Phobic.

	cerr << "Assigning polarity labels." << endl;
	
	for ( int i = 0 ; i < atom.size(); i++){
		if ( atom.at(i).elC() ){
			//carbon case
			if ( atom.at(i).isAlpha ) continue;
			if ( atom.at(i).isMain ) atom.at(i).isPolar = true; //mainchain carbonyl labelled
			if ( atom.at(i).isCofCOO ) atom.at(i).isPolar = true; //base for COO- groups
			if ( !atom.at(i).isPolar ) atom.at(i).isPhobic = true; //become hydrophobic!
			for ( int j = 0 ; j < cov.at(i).size(); j++ ){
				int k = cov.at(i).at(j); //other atom
				if ( atom.at(k).elH() || atom.at(k).elC() || atom.at(k).elS() ) continue; //not relevant
				if ( atom.at(k).isPolar ){
					atom.at(i).isPhobic = false; //has strong polar neighbour
					break; //done checking for strongly polar neighbour
				}
			}
			
		}
		if ( atom.at(i).elS() ){
			//sulphur case
			atom.at(i).isPhobic = true;
			if ( atom.at(i).nH > 0 ) atom.at(i).isDonor = true;
			//is at an acceptor? Could be
			if ( cov.at(i).size() < 4 ) atom.at(i).isAcceptor = true;
			for ( int j = 0 ; j < cov.at(i).size(); j++ ){
				int k = cov.at(i).at(j); //other atom
				if ( atom.at(k).elH() || atom.at(k).elC() || atom.at(k).elS() ) continue; //not relevant
				if ( atom.at(k).isPolar ){
					atom.at(i).isPhobic = false; //has strong polar neighbour
					break; //done checking for strongly polar neighbour
				}
			} //generally I wouldn't get this in standard residues but het groups might have something odd.
		}
		if  ( atom.at(i).elO() ){
			//oxygen case
			atom.at(i).isAcceptor = true;
			if ( atom.at(i).nH > 0 ){
				atom.at(i).isDonor = true;
			} 
			else if ( cov.at(i).size() == 1 ){
				int k = cov.at(i).at(0);//my only friend
				if ( atom.at(k).elP() || atom.at(k).isCofCOO ) atom.at(i).isChargedAcceptor = true; // I'm an O-
			}
		}
		if ( atom.at(i).elN() ){
			if ( atom.at(i).nH > 0 ){
				atom.at(i).isDonor = true; //NHsomething
				if ( cov.at(i).size() == 4 ) atom.at(i).isChargedDonor = true; //Something like NHx+
			} 
			if ( cov.at(i).size() == 3 && atom.at(i).nH == 3 ) atom.at(i).isAcceptor = true; //ammonia molecule :)
			if ( atom.at(i).isGuaniN && atom.at(i).nH > 0 ) atom.at(i).isChargedDonor = true;//arginine probably
			if ( atom.at(i).isImiN ){
				if ( atom.at(i).nH == 0 ){
					atom.at(i).isAcceptor = true; //neutral histidine probably
				}
				else{
					int on = atom.at(i).inImiWith; //my imi friend
					if ( atom.at(on).nH > 0 ){
						//I and my friend both have hydrogens, we are a charged histidine
						atom.at(i).isChargedDonor = true;
					}
				}
			}
			if ( cov.at(i).size() < 3 && !atom.at(i).isGuaniN && !atom.at(i).isImiN ) atom.at(i).isAcceptor = true; //catchall
		}
		
	}
	
     //interim report!
     cerr << "Reporting full polarity label assignments." << endl;
     for ( int i = 0 ; i < atom.size() ; i++ ){
		 if ( atom.at(i).isPhobic ){
			 if ( !atom.at(i).elS() ){
				 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as HYDROPHOBIC." << endl;
				 continue; //next!
			 }
			 if ( atom.at(i).isDonor && atom.at(i).isAcceptor ){
				 //complex sulphur case
				 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as DONOR AND ACCEPTOR AND HYDROPHOBIC." << endl;
			 }
			 else if ( atom.at(i).isDonor ){
				 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as DONOR AND HYDROPHOBIC." << endl;
			 }
			 else if ( atom.at(i).isAcceptor ){
				 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as ACCEPTOR AND HYDROPHOBIC." << endl;
			 }
			 continue; //avoid duplication of S entries
		 }
         if ( !atom.at(i).isPolar ) continue; //skip the nonpolar cases
         if ( atom.at(i).isChargedAcceptor ){
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as CHARGED ACCEPTOR." << endl;
		 } 
         else if ( atom.at(i).isChargedDonor ){
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as CHARGED DONOR." << endl;
		 }
         else if ( atom.at(i).isDonor && atom.at(i).isAcceptor ){
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as DONOR AND ACCEPTOR." << endl;
		 }		 
         else if ( atom.at(i).isDonor ){
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as DONOR only." << endl;
		 }
         else if ( atom.at(i).isAcceptor ){
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as ACCEPTOR only." << endl;
		 }
         else if ( atom.at(i).elH() ){
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as POLAR HYDROGEN." << endl;
			 //continue; //boring
		 }		 
         else {
			 cerr << "Labelled atom " << i+1 << " " << atom.at(i).atom2id() << " as POLAR NOT DONOR OR ACCEPTOR." << endl;
		 }		 

     }  
	
}


void PDBThing::seekPhobicPairs( double margin ){
	//find phobic atom 1 that is near phobic atom 2 and within cutoff and they're not bonded or same residue
	phobe.clear();
	cerr << "Identifying potential hydrophobic tethers." << endl;
	cerr << "Seeking among phobic atoms within sum of radii plus " << margin << endl;


	for ( int i = 0 ; i < atom.size() ; i++ ){
		if (!atom.at(i).isPhobic) continue; 
		vector< int > n = near ( atom.at(i), 4.5 ); // generous
		for ( int j = 0 ; j < n.size(); j++ ){
			int k = n.at(j); //other atom
			if ( !atom.at(k).isPhobic ) continue; //unsuitable
			if ( !( k > i ) ) continue; //avoid double counting
			if ( isSameResidue( atom.at(i), atom.at(k) ) ) continue; //no self dealing
			double cutoff = radius.at(i) + radius.at(k) + margin;
			double dis = d( i, k );
			if ( dis > cutoff ) continue; //too far
			PhobicPair pp;
			pp.a1 = i;
			pp.a2 = k;
			phobe.push_back( pp ); //on the list
			//cerr << "Identified possible pair: " << i+1 << " " << atom.at(i).atom2id();
			//cerr << " :# " << dis << " #: " << k+1 << " " << atom.at(k).atom2id() << endl;
		}
		
	}
	cerr << "Found a total of " << phobe.size() << " possible phobic pairs." << endl;
	
	filterPhobicPairs(); //remove redundant cases
	
	//for ( int i = 0 ; i < phobe.size() ; i++ ){
	//	cerr << fullPhobicPairLine( phobe.at( i ) ) << endl;		
	//} //this reports in form comparable to Hbond for easier cutting into FIRST-like format
		
	
}

string PDBThing::fullPhobicPairLine( PhobicPair &pp ){
	stringstream record;
	double dis = d( pp.a1, pp.a2 );
	record << pp.a1+1 << " " << pp.a2+1 << " ";
	record << " :d: " << dis;
	record << " :PHOBIC: " << atom.at( pp.a1).atom2id();
	record << " :: " << atom.at( pp.a2 ).atom2id();
	return record.str();
}

void PDBThing::filterHbondsByDistance(){
	//loop over the hb array
	//from end to beginning
	//for each case i
	//loop over the array from end to i for otherbond j
	//if j and i have the same donor and acceptor
	//assess the h-a distance for both bonds
	//and retain the case with the shorter distance
	//remembering to loop properly depending on whether i or j died.
	
	cerr << "Filtering Hbonds by distance." << endl;
	int startingSize = hb.size();
	cerr << "Starting with " << startingSize << " h-bonds." << endl;
	
	for ( int i = startingSize - 1 ; i > -1 ; i-- ){
		for ( int j = hb.size() -1 ; j > i ; j-- ){
			//check for identical donor and acceptor
			if ( hb.at(i).donor != hb.at(j).donor ) continue;
			if ( hb.at(i).acceptor != hb.at(j).acceptor ) continue;
			//cerr << "DETECTED REDUNDANT HBONDS: ";
			//cerr << "Donor " << hb.at(i).donor+1 << " " << atom.at( hb.at(i).donor ).atom2id() << " ";
			//cerr << "Acceptor " << hb.at(i).acceptor+1 << " " << atom.at( hb.at(i).acceptor ).atom2id() << endl;			
			double di = d( hb.at(i).h, hb.at(i).acceptor );
			double dj = d( hb.at(j).h, hb.at(j).acceptor );
			bool killi = false;
			if ( dj < di ) killi = true;
			
			if ( killi ){
				/*
				cerr << "REMOVING  bond through H: ";
				cerr << hb.at(i).h+1 << " " << atom.at( hb.at(i).h ).atom2id();
				cerr << " at dHA= " << di << endl;
				cerr << "RETAINING bond through H: ";
				cerr << hb.at(j).h+1 << " " << atom.at( hb.at(j).h ).atom2id();
				cerr << " at dHA= " << dj << endl;				
				*/
				hb.erase( hb.begin() + i );
				break; //skip to next i
			}
			else{
				/*
				cerr << "RETAINING bond through H: ";
				cerr << hb.at(i).h+1 << " " << atom.at( hb.at(i).h ).atom2id();
				cerr << " at dHA= " << di << endl;
				cerr << "REMOVING  bond through H: ";
				cerr << hb.at(j).h+1 << " " << atom.at( hb.at(j).h ).atom2id();
				cerr << " at dHA= " << dj << endl;				
				*/
				hb.erase( hb.begin() + j );				
			}
			
		}
		
	}
	
	cerr << "Finished filtering hydrogen bonds by distance." << endl;
	cerr << "Retained a total of " << hb.size() << " h-bonds." << endl;
	
	
	 
}

void PDBThing::filterPhobicPairs(){
	//loop over the phobe array
	//from end to beginning
	//for each case i
	//loop over the remaining array from end to i
	//examine the residue identity of the two atoms
	//if the other pair involves the same two residues as this pair
	//and if one of the atoms is the same
	//then delete the longer of the two pairs
	//because it is a redundant case
	
	//remember to set up an iterator using vector.begin() for the erase command
	
	cerr << "Filtering for redundant phobic pairs." << endl;
	int startingSize = phobe.size();
	int a11, a12, a21, a22;
	
	//test version: do it backwards!
	reverse( phobe.begin() , phobe.end() );
	
	for ( int i = startingSize - 1 ; i > -1 ; i-- ){
		//looping BACKWARDS over the array
		a11 = phobe.at(i).a1;
		a12 = phobe.at(i).a2;
		int r11 = atom.at( a11 ).res; 
		int r12 = atom.at( a12 ).res;
		int sizeNow = phobe.size();
		//bool iAmDead = false;; // activate if "killi" happens in the next loop
		for ( int j = sizeNow - 1; j > i ; j-- ){
			//looping BACKWARDS over the tail of the array
			a21 = phobe.at(j).a1;
			a22 = phobe.at(j).a2;
			bool match = false;
			//test for one atom in common
			if (  a11 == a21 || a11 == a22 ){
				match = true;
			}
			else if ( a12 == a21 || a12 == a22 ){
				match = true;
			}
			if ( !match ) continue; //no common atom
			
			//now test for residue commonality
			match = false;
			int r21 = atom.at( a21 ).res;
			int r22 = atom.at( a22 ).res;
			if ( r11 == r21 && r12 == r22 ){
				match = true;
			}
			else if ( r11 == r22 && r12 == r21 ){
				match = true;
			}
			if ( !match ) continue; //non-matching residues
			
			//okay we have a match
			//test distances
			double di, dj;
			di = d( a11, a12 );
			dj = d( a21, a22 );
			bool killi = false;
			if ( di > dj ) killi = true; //i is the one to remove
			/*
			cerr << "DETECTED REDUNDANCY: ";
			cerr << a11+1 << " " << atom.at(a11).atom2id() << " and ";
			cerr << a12+1 << " " << atom.at(a12).atom2id() << " matches ";
			cerr << a21+1 << " " << atom.at(a21).atom2id() << " and ";
			cerr << a22+1 << " " << atom.at(a22).atom2id() << endl;
			*/
			
			//now fix it;
			if ( killi ){
				/*
				cerr << "REMOVING  " << a11+1 << " " << atom.at(a11).atom2id() << " and ";
				cerr << a12+1 << " " << atom.at(a12).atom2id() << " at ";
				cerr << di << endl;
				cerr << "RETAINING " << a21+1 << " " << atom.at(a21).atom2id() << " and ";
				cerr << a22+1 << " " << atom.at(a22).atom2id() << " at ";
				cerr << dj << endl;
				*/ 
				phobe.erase( phobe.begin() + i );
				//iAmDead = true;
				break; //do not cycle to next j, rather to next i
				
			}
			else {
				/*
				cerr << "REMOVING  " << a21+1 << " " << atom.at(a21).atom2id() << " and ";
				cerr << a12+1 << " " << atom.at(a22).atom2id() << " at ";
				cerr << dj << endl;
				cerr << "RETAINING " << a11+1 << " " << atom.at(a11).atom2id() << " and ";
				cerr << a12+1 << " " << atom.at(a12).atom2id() << " at ";
				cerr << di << endl;
				*/
				phobe.erase( phobe.begin() + j );
			}
			
		}
		
	}
	
	//test: reverse the reversal!
	reverse( phobe.begin() , phobe.end() );
	cerr << "Finished filtering for redundant phobic pairs." << endl;
	cerr << "Retaining a total of " << phobe.size() << " hydrophobic tethers." << endl;
	
}

void PDBThing::seekHbondsByPolarityLabels( int energyFunction /* = 0 */, double ECutOff ){
	//finds hbonds using labels on atoms
	//optional choice of energy function
	//0 = native, 1 = "FIRST", 2 = "safety FIRST"
	//default is  0 through optional argument
	if ( energyFunction == 0 ){
		cerr << "Using NATIVE energy function( option 0 )." << endl;
	}
	else if ( energyFunction == 1 ){
		cerr << "Using FIRST's energy function( option 1 )." << endl;
	}
	else if ( energyFunction == 2 ){
		cerr << "Using safety-FIRST energy function( option 2 )." << endl;
	}
	else{
		cerr << "ERROR:  unrecognized energy function selection " << energyFunction << endl;
		return; //no hbonds found!	
	}

	
	
    hb.clear(); //null the hbonds record
    cerr << "Identifying potential hbonds by polar geometry rules." << endl;
    //establish some binned arrays to avoid massive memory inefficiency    
    vector< vector< Hbond > > hbin;
    hbin.resize( 20 ); //bins by half step from 0 to 10 
     
    //first we loop over the polarH array getting their near neighbours.
	
	for ( int p = 0 ; p < polarH.size() ; p++ ){
		int hyd = polarH.at(p); //a polar hydrogen
		int don = cov.at( hyd ).at(0); //my bonded parent
		vector< int > n = near( atom.at(hyd), 4.0 ); //potential partners
		for ( int pp = 0 ; pp < n.size() ; pp++ ){
			int acc = n.at( pp ); //possible acceptor
			if ( !atom.at(acc).isAcceptor ) continue; //ignore everything that can't accept hbonds
			if ( acc == don ) continue; //ignore my parent
			//cerr << "Polar H " << hyd+1 << " " << atom.at(hyd).atom2id() << " near acceptor ";
			//cerr << acc+1 << " " << atom.at(acc).atom2id() << " with donor ";
			//cerr << don+1 << " " << atom.at(don).atom2id() << endl;
			if ( isSameResidue( atom.at(don), atom.at(acc) ) ){
				if ( atom.at(don).isMain && atom.at(acc).isMain ){
					//cerr << "Rejecting same-residue mainchain-mainchain interaction." << endl;
					continue; //rejecting this case
				}
				//cerr << "WARNING: same-residue interaction." << endl;
			}
			//set a "salt bridge" flag for future ease of use
			bool isSB = false;
			if ( atom.at(don).isChargedDonor && atom.at(acc).isChargedAcceptor ) isSB = true;
			//if ( isSB ) cerr << "Potential SALT BRIDGE; relaxing geometry rules." << endl;
			
            //sanity check on donor-acceptor
            //now implementing three "blocking" rules:
            //(i) the D->H vector must "face" A; we can test this before looking for a base atom.
            //Failure of rule i means the donor is "blocking" H--A.
            //(ii) the B->A vector must "face" D; otherwise base is "blocking" D--A.
            //(iii) the B->A vector must "face" H; otherwise base is "blocking" H--A.
            //further note: rule (i) can be relaxed for a salt bridge?
            double lha = d( hyd, acc ) ; //H to A distance
            double lda = d( don, acc ); //D to A distance
            double lhb = 0.; //will be H to B distance later;
            double ldb = 0.; //will be D to B distance later
            double slack = 0.;
            if (isSB ) slack = 0.1;
            //cerr << "Geometry: lHA = " << lha << " , lDA = " << lda << endl;
            if ( lda + slack < lha ){
				//cerr << "Rejecting on rule (i), donor blocking hydrogen" << endl;
				continue; //dead
			}            
            //now we identify the "base" atom as the heavy covalent partner of acc which is closest to H
            //special cases:
            //if acc has only hydrogen covs (i.e. HOH) then hydrogen is an acceptable base
            //if acc has no covs at all then base is -1 //this should probably never happen
            //but of course we have those HOH residues with only O... so it could.
            int base = -1;
            if ( cov.at( acc ).size() > 0 ){
				bool  ignoreH = true; //generally
				if ( atom.at(acc).nH == cov.at(acc).size() ) ignoreH = false; //HOH case
				int mayB;
				double mayL = 0.;
				double bestL = 1000.; //ridiculous
				for ( int ppp = 0 ; ppp < cov.at(acc).size() ; ppp++){
					mayB = cov.at(acc).at(ppp);
					if ( ignoreH && atom.at(mayB).elH() ) continue;
					mayL = d( don, hyd );
					if ( mayL < bestL ){
						bestL = mayL;
						base = mayB;
					}
				}	
			}
			
			//if we're here we either have a base or we don't
			if ( base < 0 ){
				//cerr << "WARNING: no base atom; cannot check base geometry." << endl;
			}
			else{
				//we have a base
				lhb = d( hyd, base );
				ldb = d( don, base );
				//cerr << "Base " << base+1 << " " << atom.at(base).atom2id() << " ";
				//cerr << "lHB = " << lhb << " , lDB = " << ldb << endl;
				//rule ii on b,a,d
				if ( ( ldb + slack ) < lda ){
					//cerr << "Rejecting on rule (ii), base blocking donor." << endl;
					continue;
				}
				//rule iii on b,a,h
				if ( ( lhb + slack ) < lha ){
					//cerr << "Rejecting on rule (iii), base blocking hydrogen." << endl;
					continue;
				}
			}
			
			
			Hbond thisHB;
			thisHB.donor = don;
			thisHB.h = hyd;
			thisHB.acceptor = acc;
			thisHB.base = base; //could be -1 remember!
			//energyBySimpleFunctionA( thisHB ); //sets energy
			//select energy function to use...
			if ( energyFunction == 0 ){
				energySECOND( thisHB ); //using native (post-FIRST) function
				//this is probably not going to beat "safetyFIRST"
				//sois no longer the default :)				
			}
			else if ( energyFunction == 1 ){
				energyFIRST( thisHB ); //use FIRST's energy functional as best I can reconstruct it
			}
			else if ( energyFunction == 2 ){
				safetyFIRST( thisHB ); //use FIRST's energy functional as best I can reconstruct it
				//use safe LJ version in repulsive slope
			}


			cerr << "FOUND potential hydrogen bond by polar geometry with energy " << thisHB.E << endl;
		
			if (thisHB.E > ECutOff)
			{
				cerr << "Energy rejected by user cut off " << ECutOff << endl;
			}
			else
			{
			//which hbin entry should this live in?
			int mybin = - thisHB.E / 0.5; // bins are half steps in value
			//should range from 0 to 19
			if ( mybin < 0 ) mybin = 0; //catch slop around zero
			if ( mybin > 19 ) mybin = 19; //top of range slop
			//insert into appropriate bin
			int binsize = hbin.at(mybin).size();
			if ( binsize == 0 ){
				hbin.at(mybin).push_back( thisHB );
			}
			else{
				vector< Hbond >:: iterator it;
				for ( int j = 0 ; j < binsize ; j++ ){
					it = hbin.at(mybin).begin();
					if ( thisHB.E > hbin.at(mybin).at(j).E ){
						hbin.at(mybin).insert( it+j, thisHB );
						break;
					}
					//catch if not inserted already
					if ( j == binsize - 1 ) hbin.at(mybin).push_back( thisHB );
				}
			}
			}	
			//hb.push_back( thisHB );
			//cerr << terseHbondLine( thisHB ) << endl;
		}
		
	}
	
	//place bins into main array
	for ( int j = 0 ; j < 20 ; j++ ){
		double bin1 = j*0.5;
		double bin2 = bin1+0.5;
		//cerr << "Hbond strength " << -bin1 << " to " << -bin2 << " kcal/mol: " << hbin.at(j).size() << " hbonds." << endl;
		for ( int k = 0; k < hbin.at(j).size() ; k++ ){
			hb.push_back( hbin.at(j).at(k) ); //fill the main array
		} 
	}
	
	filterHbondsByDistance(); //carry out redundant distance check on main array
	
	//cerr << "Full HB list: " << hb.size() << " hbonds." << endl;
	//for ( int j = 0 ; j < hb.size() ; j++ ){
	//	cerr << fullHbondLine( hb.at( j ) ) << endl;
	//}
	
} 

/*
void PDBThing::compareDistanceWells(){
	for ( int i = 0 ; i < hb.size(); i++){
		bool isSB = false;
		if ( atom.at( hb.at(i).donor).isChargedDonor && atom.at( hb.at(i).acceptor ).isChargedAcceptor ) isSB = true;
		
		double depth, r0, a, rda;
		rda = d( hb.at(i).donor, hb.at(i).acceptor );
		
		if ( isSB ){
			//r1 = 3.0;
			//r2 = 6.0;
			depth = 10.;
			r0 = 3.2;
			a = 0.375;
		}
		else{
			//r1 = 3.0;
			//r2 = 5.0;
			depth = 8.;
			r0 = 2.8;
			a = 0.;			
		}
		
		cerr << "rDA: " << rda << " safeLJ: " << safeLJ( rda, r0, a, depth ) << " LJ: " << LJ( rda, r0, a, depth );
		cerr << " D: " << atom.at( hb.at(i).donor ).atom2id() << " H: " << atom.at( hb.at(i).h ).atom2id();
		cerr << " A: " << atom.at( hb.at(i).acceptor ).atom2id() << endl;
		
	}
}
*/

/*
void PDBThing::energyBySimpleFunctionA( Hbond &hb ){
	//assign r1, r2, depth based on bond type
	//energy = depth when lDA < r1
	//then quadratic interpolates to zero over r1->r2
	//depth is 8 for HBs in general
	//depth is 10 for salt bridge
	double r1, r2, depth;
	double range, diff; //width of well, position in well, for calculation
	double lda, lha, ldh; //for angles
	double cosT; // cos theta for dha angle
	double lab, lhb; //for base angles
	double cosP; //cos phi for hab angle
	double cos0; //ideal cos value when needed
	bool isSB = false; //
	double E; //energy of this one
	
	bool chatty = true;
	if ( chatty ) cerr << "HB energy fxn A: ";
	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ){
		isSB = true;
		if ( chatty ) cerr << "SB: ";
	}
	
	if ( isSB ){
		depth = 10.;
		r1 = 3.0;
		r2 = 6.0;
	}
	else{
		depth = 8.;
		r1 = 3.0;
		r2 = 5.0;
	}
	
	lda = d( hb.donor, hb.acceptor );
	E = cubicWell( lda, r1, r2, depth );
	
	if ( chatty ) cerr << "E(dist)= " << E << " ;";
	//then angular check if not salt bridge and got base
	if (isSB){
		hb.E = E;
		if ( chatty ) cerr << " SB done." << endl;
		return; //we are done.
	}
	
	//angular check on D-H-A angle
	lha = d( hb.h, hb.acceptor );
	ldh = d( hb.donor, hb.h ); //bond length!
	cosT = ( (ldh*ldh) + (lha*lha) - (lda*lda) ) / ( 2* ldh * lha );
	//by construction this must be negative...
	E *= cosT * cosT; // quadratic penalty reaching *0 at ninety degrees
	if ( chatty ) cerr << " cosT= " << cosT << ", E(T)= " << E << "; ";
	
	//now, if base exists, angular check on H-A-B angle
	//noting that this should ideally be the sp3 or sp2 bond angle!
	//do using difference of cosines.
	if ( hb.base < 0 ){
		hb.E = E;
		if ( chatty ) cerr << "No base." << endl;
		return; //done, no base checking as no base!
	}
	
	if ( sp2.at( hb.acceptor ) ){
		cos0 = -0.5; // -1/2 for sp2 bond angle of 120 degrees
	}
	else{
		cos0 = -0.333333333; // -1/3 for sp3 bond angle 0f 109.5 degrees
	}
	
	//got lha already, need lab, lhb
	lab = d( hb.acceptor, hb.base ); //bond length
	lhb = d( hb.h, hb.base );
	cosP = ( (lab*lab) + (lha*lha) - (lhb*lhb) ) / ( 2*lha*lab );
	if ( cosP > 0.1 ){
		E = 0; //too tight angle at base! around eighty degrees, bad
	}
	else{
		diff = cosP - cos0 ;//difference of cosine from ideal
		E *= ( 1 - ( diff * diff ) ); //quadratic penalty on cosine difference
	}
	hb.E = E;
	if ( chatty ) cerr << "cos0= " << cos0 << ", cosP= " << cosP << ", E(P)= " << E << endl;
	return; //all done!
	
} // evaluate an hbond energy by simple geometry
*/


/*
void PDBThing::energyByFunctionB( Hbond &hb ){
	//assign r0, a, depth based on bond type
	//depth is 8 for HBs in general
	//depth is 10 for salt bridge
	//then call safeLJ function to avoid FIRST short-distance fail
	double r0, a, depth, diff;
	double lda, lha, ldh; //for angles
	double cosT; // cos theta for dha angle
	double lab, lhb; //for base angles
	double cosP; //cos phi for hab angle
	double cos0; //ideal cos value when needed
	bool isSB = false; //
	double E; //energy of this one
	
	bool chatty = true;
	if ( chatty ) cerr << "HB energy fxn A: ";
	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ){
		isSB = true;
		if ( chatty ) cerr << "SB: ";
	}
	
	if ( isSB ){
		depth = 10.;
		r0 = 3.2;
		a = 0.375;
	}
	else{
		depth = 8.;
		r0 = 2.8;
		a = 0.;
	}
	
	lda = d( hb.donor, hb.acceptor );
	E = safeLJ( lda, r0, a, depth ); //
		
	if ( chatty ) cerr << "E(dist)= " << E << " ;";
	//then angular check if not salt bridge and got base
	if (isSB){
		hb.E = E;
		if ( chatty ) cerr << " SB done." << endl;
		return; //we are done.
	}
	
	//angular check on D-H-A angle
	lha = d( hb.h, hb.acceptor );
	ldh = d( hb.donor, hb.h ); //bond length!
	cosT = ( (ldh*ldh) + (lha*lha) - (lda*lda) ) / ( 2* ldh * lha );
	//by construction this must be negative...
	E *= cosT * cosT; // quadratic penalty reaching *0 at ninety degrees
	if ( chatty ) cerr << " cosT= " << cosT << ", E(T)= " << E << "; ";
	
	//now, if base exists, angular check on H-A-B angle
	//noting that this should ideally be the sp3 or sp2 bond angle!
	//do using difference of cosines.
	if ( hb.base < 0 ){
		hb.E = E;
		if ( chatty ) cerr << "No base." << endl;
		return; //done, no base checking as no base!
	}
	
	if ( sp2.at( hb.acceptor ) ){
		cos0 = -0.5; // -1/2 for sp2 bond angle of 120 degrees
	}
	else{
		cos0 = -0.333333333; // -1/3 for sp3 bond angle 0f 109.5 degrees
	}
	
	//got lha already, need lab, lhb
	lab = d( hb.acceptor, hb.base ); //bond length
	lhb = d( hb.h, hb.base );
	cosP = ( (lab*lab) + (lha*lha) - (lhb*lhb) ) / ( 2*lha*lab );
	if ( cosP > 0.1 ){
		E = 0; //too tight angle at base! around eighty degrees, bad
	}
	else{
		diff = cosP - cos0 ;//difference of cosine from ideal
		E *= ( 1 - ( diff * diff ) ); //quadratic penalty on cosine difference
	}
	hb.E = E;
	if ( chatty ) cerr << "cos0= " << cos0 << ", cosP= " << cosP << ", E(P)= " << E << endl;
	return; //all done!
	
} // evaluate an hbond energy by simple geometry
*/


void PDBThing::fillPolarH(){
	//fills an array for faster hbond searching! by seeking polar Hs from list
     polarH.clear();
     for ( int i = 0 ; i < atom.size(); i++ ){
         if ( atom.at(i).isPolar && atom.at(i).elH() ){
              polarH.push_back( i );
         }
     }
     return;
}

void PDBThing::detectGuani(){
	for ( int i = 0 ; i < atom.size(); i++ ){
		if ( !atom.at(i).elC() ) continue; //only look at carbon
		if ( cov.at(i).size() != 3 ) continue; //only sp2
		bool isG = true; //set to false if any neighbour is NOT n
		for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			int k = cov.at(i).at(j); //my bonded partner atom;
			if ( !atom.at(k).elN() ){
				isG = false; //not guani
				break; //one is enough
			}
		}
		if ( !isG ) continue; //not guani
		atom.at(i).isGuaniC = true; //I am the centre of a guanidinium group
		atom.at(i).inGuaniWith = i ;//know thyself
		for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			int k = cov.at(i).at(j); //my bonded partner atom;
			atom.at(k).isGuaniN = true; //N of CN3
			atom.at(k).inGuaniWith = i; //know parent
		}
		cerr << "Detected GUANIDINIUM carbon " << i+1 << " " << atom.at(i).atom2id() << endl;	
	}
} //scans structure labelling guanidinium groups (important for arginine salt bridges).


void PDBThing::detectImi(){
	//this one is best searched by RESIDUE not ATOM
	//because we seek two non-mainchain nitrogens
	//to allow nonstandards and HETs, we do NOT assume that only histidine matters!

	for ( int r = 0 ; r < res.size(); r++ ){
		//r is the RESIDUE in the pdbthing
		int nN = 0;
		vector< int > N; //array of nitrogens if I find any
		
		for ( int j = 0 ; j < res.at(r).atom.size() ; j++ ){
			int i = res.at(r).atom.at(j); // using the residue's atom list member
			if ( !atom.at(i).elN() ) continue; //skip all non-nitrogens!
			if ( atom.at(i).isMain ) continue; //skip not sidechain
			//if we got this far we are a sidechain nitrogen
			if ( atom.at(i).isGuaniN ) continue; //do not re-test guanidinium groups!
			nN++;
			N.push_back( i );
		}
		
		if ( nN < 2 ) continue; //pass to next residue
		cerr << "Found " << nN << " non-guanidinium non-mainchain nitrogen atoms";
		cerr << " in residue " << r+1 << " " << res.at(r).res2id() << endl;
		//first identify each possible pair of nitrogens
		for ( int a = 0 ; a < N.size(); a++ ){
			for ( int b = a+1 ; b < N.size() ; b++ ){
				int n1 = N.at(a);
				int n2 = N.at(b); //these are our two nitrogen atoms
				
				//do they have a common carbon neighbour?
				int c1 = -1;
				bool commonC = foundCommonMember( cov.at(n1), cov.at(n2), c1 );
				if ( !commonC ) break; //abandon searching this pair, not in imi geometry
				
				//okay, now test for five-ring-ism
				//this is the case if n1 has a cov member c2
				//that is carbon
				//and not c1
				//and c2 has a cov member c3
				//that is carbon
				//and not c1
				//which is a cov member for n2
				bool musashi = false; //go rin no sho
				for ( int k = 0 ; k < cov.at(n1).size(); k++ ){
					if ( musashi ) break; //already good
					int c2 = cov.at(n1).at(k);
					if ( !atom.at(c2).elC() ) continue; //not carbon, skip
					if ( c2 == c1 ) continue; //not this guy
					for ( int p = 0 ; p < cov.at(c2).size() ; p++ ){
						if ( musashi ) break; //already good
						int c3 = cov.at(c2).at(p);
						if ( !atom.at(c3).elC() ) continue; //not carbon
						if ( c3 == c1 ) continue; //not this guy (shouldn't happen anyway!)
						if ( isin( c3, cov.at(n2) ) ) musashi = true;
					}
					
				}
				
				cerr << "Identified IMIDAZOLATE ring in residue " << r+1 << " " << res.at(r).res2id() << endl;
				cerr << "Containing " << n1+1 << " " << atom.at(n1).atom2id() << " and ";
				cerr << n2+1 << " " << atom.at(n2).atom2id() << endl;
				atom.at(n1).isImiN = true;
				atom.at(n1).inImiWith = n2;
				atom.at(n2).isImiN = true;
				atom.at(n2).inImiWith = n1;				
			}
		}
		
		
		
		
	}
	
} //scans structure for imidazolate groups (important for histidines).


void PDBThing::detectCOO(){
	for ( int i = 0 ; i < atom.size(); i++ ){
		if ( !atom.at(i).elC() ) continue; //only look at carbon
		if ( cov.at(i).size() != 3 ) continue; //only sp2
		int nO = 0; //count =O bonded neighbours
		
		for ( int j = 0 ; j < cov.at(i).size() ; j++ ){
			int k = cov.at(i).at(j); //my bonded partner atom;
			if ( !atom.at(k).elO() ) continue; //pass over non-oxygens
			if ( cov.at(k).size() == 1 ) nO++; //got a suitable partner
		}
		
		if ( nO < 2 ) continue; //not COO-
		atom.at(i).isCofCOO = true; //I am Spartacus
		
		cerr << "Detected COO- carbon " << i+1 << " " << atom.at(i).atom2id() << endl;	
	}
} //scans structure for COO groups (not triggered by COOH!) for salt bridge acceptors.




string PDBThing::verboseHbondLine( Hbond &hb){
       stringstream record;
       record << hb.hbond2line() << endl;
       record << "DONOR    " << atom.at( hb.donor ).atom2id() << endl;
       record << "HYDROGEN " << atom.at( hb.h ).atom2id() << endl;
       record << "ACCEPTOR " << atom.at( hb.acceptor ).atom2id() << endl;
       if ( hb.base < 0 ){
            record << "BASE     NONE" << endl;
       }
       else{
            record << "BASE     " << atom.at( hb.base ).atom2id() << endl;
       }
       record << "ENERGY   " << hb.E;
       return record.str();
}

string PDBThing::terseHbondLine( Hbond &hb){
       stringstream record;
       record << "H " << hb.h+1 << " " << atom.at( hb.h ).atom2id() << " :: ";
       record << "A " << hb.acceptor+1 << " " << atom.at( hb.acceptor ).atom2id() << " :: ";
       record << "ENERGY " << hb.E << " E(d): " << hb.Ed << " Angular: " << hb.angular;
       return record.str();
}

string PDBThing::fullHbondLine( Hbond &hb){
       stringstream record;
		//identify type
		
		record << hb.h+1 << " " << hb.acceptor+1 << " ";
		record << hb.E << " ";
		
		if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ){
			record << ":HB:charged: ";
		}
		else{
			if ( sp2.at( hb.donor ) ){
				if ( sp2.at( hb.acceptor ) ){
					record << ":HB:sp2-sp2: ";
				}
				else{
					record << ":HB:sp2-sp3: ";
				}
			}
			else{
				if ( sp2.at( hb.acceptor ) ){
					record << ":HB:sp3-sp2: ";
				}
				else{
					record << ":HB:sp3-sp: ";
				}				
			}
		}
		record << atom.at( hb.h ).atom2id() << " :: ";
		record << atom.at( hb.acceptor ).atom2id();
		double dha = d( hb.h, hb.acceptor );
		double dda = d( hb.donor, hb.acceptor );
		record << " dHA: " << dha;
		record << " dDA: " << dda;
		
        return record.str();
}



bool PDBThing::covBonded( int i, int j ){
     if ( isin( j, cov.at(i) ) ) return true;
     //trust construction for consistency!
     //if ( isin( i, cov.at(j) ) ) return true;
     return false;
}





double LJ( double r, double r0, double a, double depth ){
	if ( r < 1E-20 ) return 0.; //just for sanity!
	double ratio = r0/(r+a);
	double r4 = ratio*ratio*ratio*ratio;
	return depth*( 5*r4*r4*r4 - 6*r4*r4*ratio*ratio );
}

double safeLJ( double r, double r0, double a, double depth ){
	if ( ( r+a) < r0 ) return -depth; //on the repulsive side of the curve
	return LJ( r, r0, a, depth );
}

double cubicWell( double r, double r1, double r2, double depth ){
	if ( r < r1 ) return -depth; //flat bottom
	if ( r > r2 ) return 0.; //outside completely
	double diff = r2 - r;
	double range = r2 - r1;
	return -1.*depth* diff * diff * diff / ( range * range * range ); // cubic from r2,0 down to r1, -D
		//E = scale * diff^3;
		//therefore -D = scale * range^3
		//therefore scale = - depth / range^3
		//therefore energy = -depth * diff^3 / range^3
}

double cosThetaTriangle( double a, double b, double c ){
	//sides are a b c
	//returns cosine of angle between sides a and b
	return ( a*a + b*b - c*c ) / (2 * a * b );
}

double cosSqExp6( double theta ){
	//this function needs THETA as an arg not cos theta
	double ct = cos( theta );
	double diff = mypi - theta;
	return ct * ct * exp( -1.* diff*diff*diff*diff*diff*diff );
}

Vector PDBThing::normalDH( Hbond &hb ){
	if ( cov.at( hb.donor ).size() < 2 ) return nullvec; //nothing to work with!
	if ( !sp2.at( hb.donor ) ) return nullvec; //not sp2
	vector< Vector > v; //vectors to work with;
	for ( int i = 0 ; i < cov.at( hb.donor ).size() ; i++ ){
		int k = cov.at( hb.donor ).at( i ); //other atom;
		if ( k == hb.h ) continue; //skip the h!
		v.push_back(  rel( hb.h, k ) ); // h to other vector
	}
	if ( v.size() < 2 ){
		//need another vector: use the H to D bond vector
		v.push_back( rel( hb.h, hb.donor ) );
	}
	//now form and norm a result;
	Vector n = v.at(0).cross( v.at(1) ); // cross product
	n.norm(); //normalise to unit
	return n; //result!
} // return the normal to an sp2 donor group

Vector PDBThing::normalAB( Hbond &hb ){
	if ( cov.at( hb.base ).size() < 2 ) return nullvec; //nothing to work with!
	if ( !sp2.at( hb.base ) ) return nullvec; //not sp2
	vector< Vector > v; //vectors to work with;
	for ( int i = 0 ; i < cov.at( hb.base ).size() ; i++ ){
		int k = cov.at( hb.base ).at( i ); //other atom;
		if ( k == hb.acceptor ) continue; //skip the A!
		v.push_back(  rel( hb.acceptor, k ) ); // h to other vector
	}
	if ( v.size() < 2 ){
		//need another vector: use the A to B bond vector
		v.push_back( rel( hb.acceptor, hb.base ) );
	}
	//now form and norm a result;
	Vector n = v.at(0).cross( v.at(1) ); // cross product
	n.norm(); //normalise to unit
	return n; //result!	
} //return the normal to an sp2 acceptor group

double PDBThing::planeNormalCosine( Hbond &hb ){
	//get back two vectors for the hbond and return a result
	//close to 1 or -1 means groups are coplanar
	if ( !sp2.at( hb.donor ) ) return 1.;
	if ( !sp2.at( hb.base ) ) return 1.;
	
	Vector nd, nb; //acceptor and base group norm vectors
	nd = normalDH( hb );
	nb = normalAB( hb );
	return nd.dot( nb ); //don't get abs value yet - handle in calling function
} //return a nonplanarity cosine for sp2 donor and acceptor groups


double PDBThing::bondPlaneCosine( Hbond &hb ){
	//compare the DH bond vector to the acceptor group normal
	//values near 0 are good
	if ( !sp2.at( hb.base ) ) return 0.; //inappropriate case
	Vector dh = rel( hb.donor, hb.h );
	dh.norm(); //unit vector equivalent
	Vector n = normalAB( hb ); //acceptor group normal
	return dh.dot( n ); //scalar product
} //return outofplane measure for DH bond versus sp2 acceptor group plane

double PDBThing::LJenergy( Hbond &hb ){
	//use LJ( r, r0, a, depth )
	double r = d ( hb.donor, hb.acceptor ); //DA distance
	double depth, r0, a;

	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ){
		depth = 10.;
		r0 = 3.2;
		a = 0.375;
	}
	else{
		depth = 8.;
		r0 = 2.8;
		a = 0.;
	}

	return LJ( r, r0, a, depth );
	
} //returns the FIRST-style LJ energy based on bond D-A distance


double PDBThing::safeLJenergy( Hbond &hb ){
	//use safeLJ( r, r0, a, depth)

	double r = d ( hb.donor, hb.acceptor ); //DA distance
	double depth, r0, a;

	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ){
		depth = 10.;
		r0 = 3.2;
		a = 0.375;
	}
	else{
		depth = 8.;
		r0 = 2.8;
		a = 0.;
	}

	return safeLJ( r, r0, a, depth );
	
}//returns LJ energy without repulsize penalty at short distance


double PDBThing::cosSqTheta( Hbond &hb ){
	//use cosThetaTriangle( a,b,c) for DHA
	double a = d( hb.donor, hb.h );
	double b = d( hb.h, hb.acceptor );
	double c = d( hb.donor, hb.acceptor );
	double t= cosThetaTriangle( a, b, c );
	return t*t;
} // returns cos-squared of the DHA angle


double PDBThing::FIRSTprefactor( Hbond &hb ){
	//use cosThetaTriangle then acos and exp...
	double a = d( hb.donor, hb.h );
	double b = d( hb.h, hb.acceptor );
	double c = d( hb.donor, hb.acceptor );
	double t= cosThetaTriangle( a, b, c );
	double tt = acos( t );
	return cosSqExp6( tt ); 
} //returns cos squared DHA angle times exp ( - ( pi-theta)^6 ) as in FIRST, bafflingly


double PDBThing::cosSqPhiPhi0( Hbond &hb, double phi0 ){
	//use cosThetaTriangle for HAB then acos and cos diff
	//safety: what if no base?
	if ( hb.base < 0 ) return 1.;
	
	double radphi0 = phi0 * degtorad;
	double a = d( hb.h, hb.acceptor );
	double b = d( hb.acceptor, hb.base );
	double c = d( hb.h, hb.base );
	double t= cosThetaTriangle( a, b, c );
	double tt = acos( t );
	
	t = tt - radphi0; // reusing t :)
	tt = cos( t ); //reusing tt :)
	return tt*tt; //cosine squared
	
} // returns penalty on HAB angle compared to phi0 ideal in degrees

double PDBThing::cosSqPhiSp2( Hbond &hb ){
	//first obtain the cos2 of the phi angle itself
	if ( hb.base < 0 ) return 1.; //safety trap
	double a = d( hb.h, hb.acceptor );
	double b = d( hb.acceptor, hb.base );
	double c = d( hb.h, hb.base );
	double cp = cosThetaTriangle(a, b, c); //this is cos phi
	double cp2 = cp*cp;
	if  ( cp2 > 0.75 ) return cp2; //done for angles above 150 degrees
	//if here, add a little bonus because 120 degrees is also favourable
	double bonus = cosSqPhiPhi0( hb, 120. ); // cos phi squared for 120 degree angle
	bonus -= 0.75; //subtract off value at 150 degrees
	return cp2 + bonus; //little fillip
}


double PDBThing::cosSqPsi( Hbond &hb ){
	//construct normal vectors on the two sp2 groups, then use dot and return cos2 equivalent
	if ( hb.base < 0 ) return 1.; //trap for bad case!
	double t = planeNormalCosine( hb );
	return t*t; //nice and simple
	//near 1 is good.
	//note can return near 0 if cannot get a normal to either group!
	//watch out...
} // returns Drieding-style planarity penalty as in FIRST


double PDBThing::cosSqOmega( Hbond &hb ){
	//construct base-group normal and DH bond vector 
	//final result is GOOD (near 1 ) if these are ORTHOGONAL
	//my function returns 0 for orthogonal...
	if ( hb.base < 0 ) return 1.; //safety trap
	double t = bondPlaneCosine( hb );
	return 1. - (t*t);
} //returns my bond/planarity penalty on DH bond vector compared to base plane

void PDBThing::energyFIRST( Hbond &hb ){
	//call LJ and appropriate angular factor
	
	//identify hb type
	bool isSB = false;
	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ) isSB = true;
	
	double E = LJenergy( hb ); //function contains its own SB trap too
	double angular = 1.0;
	
	if( isSB ){
		hb.Ed = E;
		hb.angular = 1.;
		hb.E = E;
	}
	else{
		//vary by donor, acceptor type
		if ( sp2.at( hb.donor ) ){
			if ( sp2.at( hb.acceptor ) ){
				//sp2-sp2 case
				angular = FIRSTsp2sp2( hb );
			}
			else{
				//sp2-sp3 case
				angular = FIRSTsp2sp3( hb );
			}
		}
		else{
			if ( sp2.at( hb.acceptor ) ){
				//sp3-sp2 case
				angular = FIRSTsp3sp2( hb );
			}
			else{
				//sp3-sp3 case
				angular= FIRSTsp3sp3( hb );
			}			
		}
		hb.Ed = E;
		hb.angular = angular;
		hb.E = E * angular;		 
	}
	//hb now contains its assigned energy as appropriate
	
}

void PDBThing::safetyFIRST( Hbond &hb ){
	//call safeLJ and appropriate angular factor
	
	//identify hb type
	bool isSB = false;
	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ) isSB = true;
	
	double E = safeLJenergy( hb ); //function contains its own SB trap too
	double angular = 1.0;
	
	if( isSB ){
		hb.Ed = E;
		hb.angular = 1.;
		hb.E = E;
	}
	else{
		//vary by donor, acceptor type
		if ( sp2.at( hb.donor ) ){
			if ( sp2.at( hb.acceptor ) ){
				//sp2-sp2 case
				angular = safetyFIRSTsp2sp2( hb );
			}
			else{
				//sp2-sp3 case
				angular = safetyFIRSTsp2sp3( hb );
			}
		}
		else{
			if ( sp2.at( hb.acceptor ) ){
				//sp3-sp2 case
				angular = safetyFIRSTsp3sp2( hb );
			}
			else{
				//sp3-sp3 case
				angular= safetyFIRSTsp3sp3( hb );
			}			
		}
		hb.Ed = E;
		hb.angular = angular;
		hb.E = E * angular;		 
	}
	//hb now contains its assigned energy as appropriate


} // FIRST with safeLJ

void PDBThing::energySECOND( Hbond &hb ){
	//call safeLJ and appropriate angular factor

	//identify hb type
	bool isSB = false;
	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ) isSB = true;
	
	double E = safeLJenergy( hb ); //function contains its own SB trap too
	double angular = 1.0;
	
	if( isSB ){
		hb.Ed = E;
		hb.angular = 1.;
		hb.E = E;
	}
	else{
		//vary by donor, acceptor type
		if ( sp2.at( hb.donor ) ){
			if ( sp2.at( hb.acceptor ) ){
				//sp2-sp2 case
				angular = my_sp2sp2( hb );
			}
			else{
				//sp2-sp3 case
				angular = my_sp2sp3( hb );
			}
		}
		else{
			if ( sp2.at( hb.acceptor ) ){
				//sp3-sp2 case
				angular = my_sp3sp2( hb );
			}
			else{
				//sp3-sp3 case
				angular= my_sp3sp3( hb );
			}			
		}
		hb.Ed = E;
		hb.angular = angular;
		hb.E = E * angular;		 
	}
	//hb now contains its assigned energy as appropriate

	
} //returns my modified energy function using safeLJ and angle handling


void PDBThing::energyDREIDING( Hbond &hb ){
	//call LJ and cos^4 angular term
	
	//identify hb type
	bool isSB = false;
	if ( atom.at( hb.donor ).isChargedDonor && atom.at( hb.acceptor ).isChargedAcceptor ) isSB = true;
	
	double E = LJenergy( hb ); //function contains its own SB trap too
	double angular = 1.0;
	
	if( isSB ){
		hb.Ed = E;
		hb.angular = 1.;
		hb.E = E;
	}
	else{
		double f = cosSqTheta( hb );
		angular = f*f;
		hb.Ed = E;
		hb.angular = angular;
		hb.E = E * angular;		 
	}
	//hb now contains its assigned energy as appropriate
	
}



double PDBThing::FIRSTsp3sp3( Hbond &hb ){
	//prefactor times phi penalty relative to tetrahedral angle
	return FIRSTprefactor( hb ) * cosSqPhiPhi0( hb, 109.5 );
} //FIRST case 3-3 angular factor 

double PDBThing::FIRSTsp2sp3( Hbond &hb ){
	//is square of FIRST prefactor (so cos^4 times expfactor^2
	double f = FIRSTprefactor( hb );
	return f*f;
} //FIRST case 2-3 angular factor

double PDBThing::FIRSTsp3sp2( Hbond &hb ){
	//prefactor times phi penalty relative to linear
	return FIRSTprefactor( hb ) * cosSqPhiPhi0( hb, 180.0 );
} //FIRST case 3-2 angular factor  

double PDBThing::FIRSTsp2sp2( Hbond &hb ){
	//prefactor times better out of phi and psi penalties relative to linear/coplanar
	double p = FIRSTprefactor( hb );
	double ph = cosSqPhiPhi0( hb, 180.0 ); //favours linear hbond
	double ps = cosSqPsi( hb ); //returns planarity factor
	return p * max( ph, ps );
} //FIRST case 2-2 angular factor

double PDBThing::safetyFIRSTsp3sp3( Hbond &hb ){
	//prefactor times phi penalty relative to tetrahedral angle
	return cosSqTheta( hb ) * cosSqPhiPhi0( hb, 109.5 );
} //FIRST case 3-3 angular factor 

double PDBThing::safetyFIRSTsp2sp3( Hbond &hb ){
	//is square of FIRST prefactor (so cos^4 times expfactor^2
	double f = cosSqTheta( hb );
	return f*f;
} //FIRST case 2-3 angular factor

double PDBThing::safetyFIRSTsp3sp2( Hbond &hb ){
	//prefactor times phi penalty relative to linear
	return cosSqTheta( hb ) * cosSqPhiPhi0( hb, 180.0 );
} //FIRST case 3-2 angular factor  

double PDBThing::safetyFIRSTsp2sp2( Hbond &hb ){
	//prefactor times better out of phi and psi penalties relative to linear/coplanar
	double p = cosSqTheta( hb );
	double ph = cosSqPhiPhi0( hb, 180.0 ); //favours linear hbond
	double ps = cosSqPsi( hb ); //returns planarity factor
	return p * max( ph, ps );
} //FIRST case 2-2 angular factor

double PDBThing::my_sp3sp3( Hbond &hb ){
	//return cos-squared theta times cos-squared of phi penalty relative to tetrahedral angle
	return cosSqTheta( hb ) * cosSqPhiPhi0( hb, 109.5 );
} //my case 3-3 angular factor 

double PDBThing::my_sp2sp3( Hbond &hb ){
	//return identical to sp3sp3 case by my logic!
	return cosSqTheta( hb ) * cosSqPhiPhi0( hb, 109.5 );
} //my case 2-3 angular factor

double PDBThing::my_sp3sp2( Hbond &hb ){
	//return cos-squared theta times the better of phi penalties relative to linear or trigonal angle
	return cosSqTheta( hb ) * cosSqPhiSp2( hb );
} //my case 3-2 angular factor  

double PDBThing::my_sp2sp2( Hbond &hb ){
	//return cos-squared theta times better of:
	//phi penalty with sp2 120 degree correction
	//bond planarity penalty
	double ph = cosSqPhiSp2( hb );
	double om = cosSqOmega( hb );
	return cosSqTheta( hb ) * max( ph, om );
} //my case 2-2 angular factor  

Vector PDBThing::residueCentralPoint( int i ){
	//the list of atoms is res.at(i).atom list
	//skip hydrogens elH();
	//return average position of heavy atoms
	//for use as a node in ENM
	Vector result; //create null vector
	int nHeavy = 0;
	for ( int j = 0 ; j < res.at(i).atom.size() ; j++ ){
		int k = res.at(i).atom.at(j); //this atom in master list.
		if (  atom.at(k).elH() ) continue; //skip hydrogens
		nHeavy++;
		result += atom.at(k).pos ;
	}
	if ( nHeavy == 0 ) return result; //danger - null vector, no heavy atoms! - shouldn't happen!
	result /= nHeavy;
	return result; //average position of heavy atoms in this residue.
}

Vector PDBThing::residueAlphaPoint( int i ){
	//the list of atoms is res.at(i).atom list
	//return CA position
	//for use as a node in ENM
	Vector result; //create null vector
	if ( !res.at(i).isStandard ) return result; //danger - null vector for nonStandard - check before calling!
	for ( int j = 0 ; j < res.at(i).atom.size() ; j++ ){
		int k = res.at(i).atom.at(j); //this atom in master list.
		if ( atom.at(k).isAlpha ){
			result = atom.at(k).pos ;
			return result; //CA position in this residue.
		}
	}
	return result;//null vector returned if I fail badly....
}

vector< Vector > PDBThing::nodesForENM(){
	// return one node per residue
	//return CA for standard
	//return resCP for other
	//use isStandard boolean in residue object
	vector< Vector > result;
	Vector current;
	
	for ( int i = 0 ; i < res.size() ; i++ ){
		current = nullvec;
		if ( res.at(i).isStandard ){
			current = residueAlphaPoint( i );
			//cerr << "Residue " << i+1 << " gives node position " << current << " from ";
			//cerr << atom.at(  res.at(i).atom.at(0)  ).atom2line() << endl;
		}
		else{
			current = residueCentralPoint( i );
			//cerr << "Residue " << i+1 << " gives node position " << current << endl;
		}
		result.push_back( current );
	}
	return result;
}

vector< Vector > PDBThing::allXYZs(){
	//return pos for every atom
	vector< Vector > result;
	result.resize( atom.size() );
	for ( int i = 0 ; i < atom.size() ; i++ ){
		result.at(i) = atom.at(i).pos;
	}
	return result;
}

vector< Vector > PDBThing::fileXYZs( int i ){
	vector< Vector > result;
	result.resize( file.at(i).atom.size() );

	for ( int j = 0 ; j < file.at(i).atom.size() ; j++ ){
		int k = file.at(i).atom.at(j);
		result.at(j) = atom.at(k).pos;
	}
	return result;
}

vector< Vector > PDBThing::modelXYZs( int i ){
	vector< Vector > result;
	result.resize( model.at(i).atom.size() );

	for ( int j = 0 ; j < model.at(i).atom.size() ; j++ ){
		int k = model.at(i).atom.at(j);
		result.at(j) = atom.at(k).pos;
	}
	return result;
}

vector< Vector > PDBThing::chainXYZs( int i ){
	vector< Vector > result;
	result.resize( chain.at(i).atom.size() );

	for ( int j = 0 ; j < chain.at(i).atom.size() ; j++ ){
		int k = chain.at(i).atom.at(j);
		result.at(j) = atom.at(k).pos;
	}
	return result;
}

vector< Vector > PDBThing::residueXYZs( int i ){
	vector< Vector > result;
	result.resize( res.at(i).atom.size() );

	for ( int j = 0 ; j < res.at(i).atom.size() ; j++ ){
		int k = res.at(i).atom.at(j);
		result.at(j) = atom.at(k).pos;
	}
	return result;
}


bool PDBThing::writePseudoPDBforENM( string filename ){
	//write all the nodes as CA atoms with residue information
	//even if they are fictitious het group nodes
	
    //open the file for writing;
    ofstream outf;
    outf.open( filename.c_str(), ios::out );
    if ( outf.fail() ){
         cerr << "ERROR: cannot open file " << filename << endl;
         return false;
	}
	
	PDBAtom wop;

	vector< Vector > nodes = nodesForENM();
	for ( int i = 0 ; i < res.size(); i++ ){
		int j =res.at(i).atom.at(0); //first atom in the residue has residue information
		wop = atom.at(j); //copy the atom to get all the relevant info
		wop.recname = "ATOM  "; //overwrites HETATM
		wop.name = " CA "; //pseudo CA
		wop.element = " C";
		wop.serial = i+1 ; //just to avoid jumping about? Consider deletion
		wop.pos = nodes.at(i); //node position
		outf << wop.atom2line() << endl;
	}
	
	outf << "END   " << endl;
    outf.close();
    return true;
}

bool PDBThing::writeXYZforENM( string filename ){
	//write all the nodes as XYZ positions
	
    //open the file for writing;
    ofstream outf;
    outf.open( filename.c_str(), ios::out );
    if ( outf.fail() ){
         cerr << "ERROR: cannot open file " << filename << endl;
         return false;
	}

	vector< Vector > nodes = nodesForENM();
	for ( int i = 0 ; i < res.size(); i++ ){
		outf << nodes.at(i)  << endl;
	}

    outf.close();
    return true;
}

bool PDBThing::writePseudoPDBforENM( string filename, vector< Vector > &nodes ){
	//write all the nodes as CA atoms with residue information
	//even if they are fictitious het group nodes
	
    //open the file for writing;
    ofstream outf;
    outf.open( filename.c_str(), ios::out );
    if ( outf.fail() ){
         cerr << "ERROR: cannot open file " << filename << endl;
         return false;
	}
	
	if ( nodes.size() != res.size() ){
		cerr << "ERROR: mismatch between nodes and residues." << endl;
		return false;
	}
	
	PDBAtom wop;

	//vector< Vector > nodes = nodesForENM();
	for ( int i = 0 ; i < res.size(); i++ ){
		int j =res.at(i).atom.at(0); //first atom in the residue has residue information
		cerr << "Starting with info of atom " << j+1 << " : " <<  atom.at( j ).atom2line() << endl;
		wop = atom.at(j); //copy the atom to get all the relevant info
		wop.recname = "ATOM  "; //overwrites HETATM
		wop.name = " CA "; //pseudo CA
		wop.element = " C";
		wop.serial = i+1 ; //just to avoid jumping about? Consider deletion
		wop.pos = nodes.at(i); //node position
		cerr << "Writing info for node atom: " << wop.atom2line() << endl;
		outf << wop.atom2line() << endl;
	}
	
	outf << "END   " << endl;
    outf.close();
    return true;
}

bool PDBThing::writeXYZforENM( string filename, vector< Vector > &nodes ){
	//write all the nodes as XYZ positions
	
    //open the file for writing;
    ofstream outf;
    outf.open( filename.c_str(), ios::out );
    if ( outf.fail() ){
         cerr << "ERROR: cannot open file " << filename << endl;
         return false;
	}

	for ( int i = 0 ; i < res.size(); i++ ){
		outf << nodes.at(i)  << endl;
	}

    outf.close();
    return true;
}

vector< vector< int > > PDBThing::smallRC(){
	cerr << "Building local rigid clusters." << endl;
	//step i; build candidate clusters centred on each atom with neighbours
	vector< vector< int > > candidate;
	for ( int i = 0 ; i < atom.size() ; i++ ){
		vector< int > candid;
		candid.push_back( i );
		mergeAintoB( cov.at(i), candid );
		candidate.push_back( candid );
	}
	//so in this version "candidate" is exactly the size of "atom"
	//step ii; unify based on bars
	//by looking at the bars entry between candidate heads
	vector< vector< int > > result; //
	int nc = candidate.size(); //number of initial clusters to examine
	vector< int > bondto ;//unify with this cluster if not self
	bool gotmates = false;
	bondto.resize( nc );
	
	for ( int i = 0 ; i < nc ; i++ ){
		bondto.at(i) = i;
	}
	
	bool changing = true; //track progress or stasis
	int rounds = 0;
	
	while ( changing ){
		changing = false;
		rounds++;
		cerr << "Garibaldi performing unification round: " << rounds << endl;
		for ( int i = 0; i < nc ; i++ ){
			//i is the index and also the first entry!
			for ( int j = 0 ; j < bars.at(i).size() ; j++ ){
				//need the atom and the bar count
				int k = cov.at(i).at(j); //atom
				int b = bars.at(i).at(j); //bars
				//do not unify if bars not 6:
				if ( b < 6 ) continue; //skipping
				
				//so if we're here, we are unifying these clusters
				gotmates = true;
				int italy = bondto.at(i);
				if ( bondto.at(k) < italy ) italy = bondto.at(k);
				if ( italy != bondto.at(i) || italy != bondto.at(k) ) changing = true;
				bondto.at(i) = italy;
				bondto.at(k) = italy;
				//bondtos propagate like burning algorithm
			}
		}
		
	}

	if ( !gotmates ){
		cerr << "Garibaldi found nothing to do: washing shirts." << endl;
		result = candidate; return result; // don't do anything
	}
	//first pass;all clusters not bonded to self unify with their bondto partner.
	for ( int i = 0; i < nc; i++ ){
		if ( bondto.at(i) != i ){
			int borg = bondto.at(i); //feed candidate(i) into candidate( borg );
			mergeAintoB( candidate.at(i), candidate.at( borg ) );
		}
	}	
	
	//second pass; only clusters bonded to self are pushed into result
	for ( int i = 0; i < nc; i++ ){
		if ( bondto.at(i) != i ){
			continue;
		} // don't push these;
		result.push_back (  candidate.at(i) ); // borg clusters pushed into result
	}
	cerr << "Garibaldi completed risorgimento. Grazie, Giuseppe." << endl;	

	//debug output:
	if ( false ){
		cerr << "Rigid cluster listing: " << endl;
		for ( int i = 0 ; i < result.size() ; i++ ){
			cerr << "RC" << i+1 <<":";
			for ( int j = 0 ; j < result.at(i).size() ; j++){
				int k = result.at(i).at(j);
				cerr << " " << k+1;
			}
			cerr << endl;
		}
	}
	
	return result;
}

bool PDBThing::writeRC( string filename, vector< vector< int > > &RC ){
    //open the file for writing;
    ofstream outf;
    outf.open( filename.c_str(), ios::out );
    if ( outf.fail() ){
         cerr << "ERROR: cannot open file " << filename << endl;
         return false;
	}
	outf << RC.size() << endl; //"Number of rigid bodies" on first line
		for ( int i = 0 ; i < RC.size() ; i++ ){
			//outf << "RC" << i+1 <<":";
			outf << RC.at(i).size(); //?check with Tom if this is his convention
			for ( int j = 0 ; j < RC.at(i).size() ; j++){
				int k = RC.at(i).at(j);
				outf << " " << k+1;
			}
			outf << endl;
		}
    outf.close();
    return true;
}

bool PDBThing::writeFrame( string rider ){
	bool goodwrite = true;
	for ( int i = 0 ; i < file.size() ; i ++ ){
		cerr << "Working on file " << i+1;
		string name = file.at(i).name ;
		cerr << " with own name " << name << endl;
		stringstream onamestream;
		onamestream << name << "." << rider << ".pdb";
		string oname = onamestream.str();
		goodwrite = writePDBFile( i, oname );
		if (!goodwrite){
			cerr << "Disastrous failure while trying to write to " << oname << endl;
			return false;
		}
	}
	return true;
}