#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
#include <cmath>
#include <stdlib.h>
#include <cctype>
#include <cassert>
#include <unistd.h> //for chdir
#include <time.h>
#include <sys/time.h>
#include <vector>
#include <algorithm>
#include <functional> //for transform
#include "nr3.h" //for shape from "Numerical Recipes in C" W. H. Press et al
#include "svd.h" //for shape

//NEEDS C++11

#define NUMBINS 40 //R HISTO
#define MAXDIST 40 //R HISTO
#define NUMBINSRG 20 //RG HISTO set to MAXDISTRG (see explanation below under HISTO
#define MAXDISTRG 20 //RG HISTO "standard one used for PEP_X
#define NUMBINSETOE 40 //RG HISTO set to MAXDISTRG (see explanation below under HISTO
#define MAXDISTETOE 40 //RG HISTO "standard one used for PEP_X
//how many samps of input Ens to keep instead of use them all
//#define NUMSAMPSUSED 10
#define WRITEINTERVAL 10 //for pdbs and RG_EXT
#define RANGE  0.0//RN (-range/2,range/2) //if want to "fuzz out" phi/psi angles.
#define DISCUT 2.4  //similar to min dis in MD FF - van der waals >3rd neighbor exclusion

using namespace std;

//program DihedToCC.cpp
//PURPOSE: Convert dihed angles to CCs and evaluate properties.
//SET CCs of atoms A,B,C and dihed angles, bond distances, angles
//IN: Ens*.dat: Into array  Phi[numSamp][numRes] and Psi[numSamp][numRes] 
//(numSamp samples of Phi and Psi diheds for numRes residues generated externally)

//METHOD: init with CCs of A,B,C and: CD bond distances, BCD angles and some BC dihed. Return D CCs
//Then, cp B,C,D CCs to A,B.C CCs, respectively. For read in BC dihed, CD bond dist, BCD angle
//get D CCs. Also write to E the carbonyl oxygen.
//Iterate this to form the chain CCs.
//METHOD generate all atom distances for vdw clashes.
//Use mapAt[m]=at+1 so written out atoms are 1-based for comparison to pdbs
//METHOD then delete pdbs based on clash rejection.
//N.B. skip O1-O2 and O1-N2 from vdw clash because init and read in dihPsi aren't "connected"

//get RG (use same mass value since N,C,CA,O close enough) and EtoE 
//RG = sqrt (sum squares of pos (CC[at][i]) - COM)/sqrt(numAtoms)     
//EtoE dist of CA res first to last

//N.B. bondAngle = (bondAngle-180); That means geometry A-B-C (in a row)
//corresponds to input angle of 180 degrees. Then 90 degrees is B-C rotated up and 270
//degrees B-C rotated down relative to A-B.
//init  psi_0 and cycle over (omega1 phi_1 psi_1) (omega2 phi_2 psi_2)....

//Get "HB" for alpha helix using res1(O)-res4(N) and for PPII res1(O)-res3(N). 
//Using atom 1 etc numbering
//1-4: starts with atoms 4 - 13 (res1-res4) 
//1-3: starts with atoms 4 - 9 (res1-res3) 

//OUTFILES
//OUT: cout has list of dihed angles used for each state and if they exist 
//distances less than disCut.
//OUT: /DATA/*.pdb for each kept state
//N.B. have to mv /DATA to something, else program won't work.
//set range if want to "fuzz out" phi and psi angles
//OUT: *_RG_EXT.txt  RG delta and EtoE by kept sample
//OUT:  *_PofR.txt P(R) pdf of R (distance distribution)
//OUT:  *_PofRG.txt P(RG) pdf of RG distribution
//OUT:  *_PofEtoE.txt P(EtoE) pdf of EtoE distribution
/*
  COUT: End of log
  Shape information - based on expressions in:
  Vitalis, A.; Wang, X. L.; Pappu, R. V., Quantitative Characterization of 
  Intrinsic Disorder in Polyglutamine: Insights from Analysis 
  Based on Polymer Theories. Biophys. J. 2007, 93, 1923-1937
  Dima, R. I.; Thirumalai, D., Asymmetry in the Shapes of Folded and Denatured 
  States of Proteins. J. Phys. Chem. B 2004, 108, 6564-6570
  RGsq: 
  Shapes[0]=RGsq
  delta: asphericity (ranges (0,1)
  Shapes[1] = 1.0-3.0*(eigVals[0]*eigVals[1]+eigVals[1]*eigVals[2]+eigVals[2]*eigVals[0])
  /(Shapes[0]*Shapes[0]);
  S: shape parameter (ranges -1/4, 2) 
  Shapes[2]=27.0*numerator/denom;
  double lambAve=Shapes[0]/3.0;
  double numerator= (eigVals[0]-lambAve)*(eigVals[1]-lambAve)*(eigVals[2]-lambAve);
  double denom = pow(Shapes[0],3);

  percent HB14 (alphaHelix) 
  percent HB13 (PPII)

*/

void refoldNRF(double*, double*, double*, double*, double, double, double, double);
void writeIt(char*, int, int, double*); 
//void fwriteIt(ofstream &, char*, int, int, double*); 
void fwriteIt(ofstream &, char*, int, int, double*, vector<string>); 
double getdist(double*, double*);
vector<double> getPofR(double*,  int, int, double);
void getCovMatrix(double**, int, double**); //in CC and numAtoms; out covMatrix 
vector <double> getEigvals(double**, int);
vector <double> getShapes(vector <double>, int);
vector<double> getHistoRG(double, int, double);
vector<double> getHistoEtoE(double, int, double);
double getJHNH(double); //for NMR
string oneToThree(string); //convert fasta to three letter

main()
{
  string sdirInputs ="/path/analyticEnsOutDirectory";
  string iFileName="Ensemble_nR_9_TCB.dat"; //name in the above dir
  chdir(sdirInputs.c_str());
  system("echo -n '1. Data In Directory is '; pwd");
  //put results in named dir
  string sdirOut ="/path/DihedToCCOutDirectory";
  int numSamps; //read from IFILE
  int numRes; //read from, ILFILE
  int numDoub; //numRes-1
  double range=RANGE; //for RN (-range/2,range/2)
  double disCut=DISCUT;//for vdw

  //for histo P(R)
  int numBins=NUMBINS; 
  double maxDist=MAXDIST;
  double delBin =maxDist/numBins;
  double delBin2=delBin/2.0;
  vector <double> cumeHisto(numBins); //to cume histo for P(R)
  //for histo of P(RG)
  int numBinsRG=NUMBINSRG; 
  double maxDistRG=MAXDISTRG;
  double delBinRG =maxDistRG/numBinsRG;
  double delBin2RG=delBinRG/2.0;
  vector <double> cumeHistoRG(numBinsRG); //to cume histo for p(RG)
  double cumeRG=0.0; //cume RG for write to cout
  int numBinsEtoE=NUMBINSETOE;
  double maxDistEtoE=MAXDISTETOE;
  double delBinEtoE =maxDistEtoE/numBinsEtoE;
  double delBin2EtoE=delBinEtoE/2.0;
  vector <double> cumeHistoEtoE(numBinsEtoE); //to cume histo of EtoE

  int numShapes=3; //how many shape parameters to write
  vector <double> cumeShapes(numShapes,0);

  int numSampsKept=0; //count how many samps kept (no overlap)
  //read Phi Psi from Ensemble generated by realization_dataBase.cpp
  ifstream iFile;
  //string iFileName="Ensemble_nR_9_TCB.dat";
  iFile.open(iFileName.data(),ios::in);
  assert(iFile.is_open());
  string drop;
  string sequenceName;
  string line;
  while (iFile >> std::ws && std::getline(iFile, line)); // skip empty lines and read last line
  numSamps=std::stoi(line); //total number samps in Ensemble for realz
  cout<<"there are "<<numSamps<<" samples in the Esemble"<<endl; 
  //numSamps = NUMSAMPSUSED; // If want to change number samples read<total in input
  cout<<"with "<<numSamps<<" used"<<endl;
  iFile.close();
  iFile.open(iFileName.data(),ios::in);
 
  iFile>>drop>>numRes>>drop>>sequenceName;
  numDoub=numRes-1;
  int lmax=numDoub*3;
  int numAtoms=4*(numDoub+1);//for CC 4 atoms * number of residues
  double **Phi; //Phi[numSamps][numRes]
  Phi = new double*[numSamps];
  for(int ns=0;ns<numSamps;ns++)
    Phi[ns] = new double[numDoub];
  double **Psi;
  Psi = new double *[numSamps];
  for(int ns=0;ns<numSamps;ns++)
    Psi[ns] = new double[numDoub];
  //over residues
  for(int ns=0;ns<numSamps;ns++)
    {
      for(int nd=0;nd<numDoub;nd++)
	iFile>>Phi[ns][nd];
      for(int nd=0;nd<numDoub;nd++)
	iFile>>Psi[ns][nd];
    }
  iFile.close();
  /*
  //debug in case want to check input properly read
  for(int ns=0;ns<numSamps;ns++){
  for(int nd=0;nd<numDoub;nd++)
  cout<<Phi[ns][nd]<<" ";
  cout<<endl;;
  for(int nd=0;nd<numDoub;nd++)
  cout<<Psi[ns][nd]<<" ";
  cout<<endl;
  }
  exit(-1);
  */
  double **CC; //CC[numAtoms][3]
  CC = new double *[numAtoms];
  for(int at=0;at<numAtoms;at++)
    CC[at] = new double[3];
  //set up arrays for keeping desired distances
  int numDist = numAtoms*(numAtoms-1)/2;   
  double *distArray;
  distArray = new double[numDist];
  int *mapAt, *mapAtp;
  mapAt = new int[numDist]; //dim numDist to go with distArray[numDist]
  mapAtp = new int[numDist];
  double *COM; // COM first
  COM = new double[3];
  double *distArrayAll;
  distArrayAll = new double[numDist];
  double *atCC;
  atCC = new double[3];

  //set up HB
  int nHB=numRes-4; //how many 14 HBs
  int *cntHB;
  cntHB = new int[nHB];
  for(int nh=0;nh<nHB;nh++)
    cntHB[nh]=0;
  int nHB13=numRes-3; //how many 13 HBs
  int *cntHB13;
  cntHB13 = new int[nHB13];
  for(int nh=0;nh<nHB13;nh++)
    cntHB13[nh]=0;
  //set up for NMR
  double *jHNH;
  double *jHNHSQ;
  jHNH = new double[numDoub];
  jHNHSQ = new double[numDoub];
  for(int nd=0;nd<numDoub;nd++){ 
    jHNH[nd]=0.0;
    jHNHSQ[nd]=0.0;
  } 

  chdir(sdirOut.c_str()); //put ofiles here 
  system("echo -n '2. Data Out Directory is '; pwd");
  //for write OUTs
  //ofiles 
  string Sens=static_cast<ostringstream *> (&(ostringstream() <<numSamps/1000) )->str();

  string ofname=sequenceName+"_ENS_"+Sens+"K.pdb"; //OUTFILE all in one file
  ofstream ofile(ofname.data());
  if(!ofile.is_open()){
    cout<<ofname<<" not available"<<endl;
    exit(-1);
  }
  ofname=sequenceName+"_ENS_"+Sens+"K_KEEP.pdb"; //OUTFILE all in one file
  ofstream ofileKP(ofname.data());
  if(!ofileKP.is_open()){
    cout<<ofname<<" not available"<<endl;
    exit(-1);
  }
  string ofnameRG=sequenceName+"_ENS_"+Sens+"K_RG_EXT.txt";
  ofstream ofRG(ofnameRG.data());
  if(!ofRG.is_open()){
    cout<<ofnameRG<<" not available"<<endl;
    exit(-1);
  }
  ofRG<<"ns\t"<<"RG\t"<<"delta\t"<<"EtoE\t"<<endl;
  string ofnamePofR=sequenceName+"_ENS_"+Sens+"K_PofR.txt";
  ofstream ofilePofR(ofnamePofR.data());
  if(!ofilePofR.is_open()){
    cout<<ofnamePofR<<" not available"<<endl;
    exit(-1);
  }
  ofilePofR<<"R"<<"\t"<<"P(R)"<<endl;
  ofilePofR<<0.00<<"\t"<<0.00<<endl; //always

  string ofnamePofRG=sequenceName+"_ENS_"+Sens+"K_PofRG.txt";
  ofstream ofilePofRG(ofnamePofRG.data());
  if(!ofilePofRG.is_open()){
    cout<<ofnamePofRG<<" not available"<<endl;
    exit(-1);
  }
  string ofnamePofEtoE=sequenceName+"_ENS_"+Sens+"K_PofEtoE.txt";
  ofstream ofilePofEtoE(ofnamePofEtoE.data());
  if(!ofilePofEtoE.is_open()){
    cout<<ofnamePofEtoE<<" not available"<<endl;
    exit(-1);
  }
  string ofnameNMR=sequenceName+"_ENS_"+Sens+"K_NMR.txt";
  ofstream ofileNMR(ofnameNMR.data());
  if(!ofileNMR.is_open()){
    cout<<ofnameNMR<<" not available"<<endl;
    exit(-1);
  }
  string ofnameHB=sequenceName+"_ENS_"+Sens+"K_HB.txt";
  ofstream ofileHB(ofnameHB.data());
  if(!ofileHB.is_open()){
    cout<<ofnameHB<<" not available"<<endl;
    exit(-1);
  }

  //put in headers
  ofilePofRG<<"RG"<<"\t"<<"P(RG)"<<endl;
  ofilePofEtoE<<"EtoE"<<"\t"<<"P(EtoE)"<<endl;

  //convert fasta to three letter
  vector <string> name;
  char tempCString[2];   //do a fasta string
  string str=sequenceName;
  for(std::string::iterator it = str.begin(); it != str.end(); ++it) {
    tempCString[0]=*it; //convert  to  C-string of length 1
    string answ = oneToThree(tempCString);
    answ=answ.substr(0, 3); //pull out 3 letter code
    name.push_back(answ);
  }

  //constants set phi/psi to alpha helix omega to planar
  //for carbonyl O (CO) dihed use read dihPsi+180 
  double bndOmg=1.46; double bndPhi=1.51; double bndPsi=1.33;
  double angOmg=122; double angPhi=111; double angPsi=116;
  double dihOmg=180; double dihPhiBase=-60; double dihPsi=0;
  double bndCO=1.24; double angCACO=120.5; double dihNCACO=180+dihPsi; //place O
  double dihPhi;
 
  //atom CC vectors
  double *a, *b, *c, *d, *e;
  a = new double[3];
  b = new double[3];
  c = new double[3];
  d = new double[3];
  e = new double[3]; //e for write O

  //CCs stored for eliminate vdw clashes
  for(int at=0;at<numAtoms;at++)
    for(int i=0;i<3;i++)
      CC[at][i]=0.0;
  
  cout<<"numAtoms="<<numAtoms<<endl;

  //for shapes
  double **covMatrix;
  covMatrix = new double*[3];
  for(int i=0;i<3;i++)
    covMatrix[i] = new double[3];
  MatDoub covMatrixTot(3,3); 
  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
      covMatrixTot[i][j]=0.0;
  
  //for Random Number if "fuzz"
  timeval val;
  gettimeofday(&val,NULL);
  srand(val.tv_usec);//get the microseconds for seed

  //loop over numSamps
  //timer start
  clock_t start = clock();
  int writeInterval=WRITEINTERVAL;
  double sumDiag=0.0; //DEBUG

  for(int ns=0;ns<numSamps;ns++) //BEGIN SAMP LOOP 
    {
      //write header to ofile
      if(ns%writeInterval==0)
	ofile<<"MODEL"<<"        "<<setiosflags(ios::left)<<ns+1<<resetiosflags(ios::left)<<endl;
      
      a[0]=-1.0; a[1]=1.0; a[2]=0.0; //N
      b[0]=-0.0; b[1]=0.0; b[2]=0.0; //CA
      c[0]=1.5;  c[1]=0.0; c[2]=0.0; //C

      //first time set psi
      refoldNRF(a, b, c, d, bndPsi, angPsi,dihPsi, 1.0); //set psi 
      int countAt=1;
      int countRes=1;
      char atN[]="N";
      if(ns%writeInterval==0) fwriteIt(ofile,atN,countAt,countRes,a,name);
      for (int i=0;i<3;i++)
	CC[countAt-1][i]=a[i];
      countAt++;
 
      char atCA[]="CA";
      if(ns%writeInterval==0) fwriteIt(ofile,atCA,countAt,countRes,b,name);
      for (int i=0;i<3;i++)
	CC[countAt-1][i]=b[i];
      countAt++;

      char atC[]="C";
      if(ns%writeInterval==0) fwriteIt(ofile,atC,countAt,countRes,c,name);
      for (int i=0;i<3;i++)
	CC[countAt-1][i]=c[i];
      countAt++;

      refoldNRF(a, b, c, e, bndCO, angCACO,dihNCACO, 1.0); //set O first time
      char atO[]="O";
  
      e[2]=0.0; //Kludge for getting 0.00  instead of -1.09e-11 
      if(ns%writeInterval==0) fwriteIt(ofile,atO,countAt,countRes,e,name);
      for (int i=0;i<3;i++)
	CC[countAt-1][i]=e[i];
 
      countAt++;
      countRes++;

      //for l%3==0 omega write N, l%3==1 phi write CA, l%3==2 psi write C and O
      for (int l=0;l<lmax;l++)
	{
	  for (int i=0;i<3;i++)
	    {
	      a[i]=b[i];
	      b[i]=c[i];
	      c[i]=d[i];
	    }
	  //get value of each dih
	  int lby3=l/3; //pick out every third l for Phid[ns][nr] //nr residue
	  if(l%3==0)
	    {
	      dihPhi = Phi[ns][lby3];
	      double randNum = ((double)rand())/(RAND_MAX); //(0,1)
	      //shift to (-range,range) with range=RANGE/2
	      randNum *=range;
	      randNum -=range/2;
	      dihPhi += randNum;
	      dihPsi = Psi[ns][lby3];
	      randNum = ((double)rand())/(RAND_MAX); //(0,1)
	      randNum *=range;
	      randNum -=range/2;
	      dihPsi += randNum;
	    }
	  //get CC next atom
	  if(l%3==0)
	    {
	      refoldNRF(a, b, c, d, bndOmg, angOmg, dihOmg, 1.0);
	      if(ns%writeInterval==0) fwriteIt(ofile,atN,countAt,countRes,c,name);
	      for (int i=0;i<3;i++)
		CC[countAt-1][i]=c[i];
	    }
	  else if(l%3==1)
	    {
	      refoldNRF(a, b, c, d, bndPhi, angPhi,dihPhi, 1.0);
	      if(ns%writeInterval==0) fwriteIt(ofile,atCA,countAt,countRes,c,name);
	      for (int i=0;i<3;i++)
		CC[countAt-1][i]=c[i];
	    }
	  else if(l%3==2)
	    {
	      refoldNRF(a, b, c, d, bndPsi, angPsi,dihPsi, 1.0);
	      if(ns%writeInterval==0) fwriteIt(ofile,atC,countAt,countRes,c,name);
	      for (int i=0;i<3;i++)
		CC[countAt-1][i]=c[i];
	      countAt++; //"fourth" atom O
	      refoldNRF(a, b, c, e, bndCO, angCACO,180+dihPsi, 1.0); //evaluate O CCs
	      if(ns%writeInterval==0) fwriteIt(ofile,atO,countAt,countRes,e,name);
	      for (int i=0;i<3;i++)
		CC[countAt-1][i]=e[i];
	    }
  
	  countRes=countAt/4;
	  countRes++;
	  countAt++;
	} //end l loop

      //write termination to ofile
      if(ns%writeInterval==0) ofile<<"TER"<<endl;
      if(ns%writeInterval==0) ofile<<"END"<<endl;
  
      //get distArray   
      for(int m=0;m<numDist;m++)
	distArray[m]=-1;
      for(int m=0;m<numDist;m++)
	{
	  mapAt[m]=-1;
	  mapAtp[m]=-1;
	}
      //get all more than 3rd neighbor distances
      int  m=0;
      double dist=0.0;
      for(int at=0;at<numAtoms-1;at++)
	for (int atp=at+1;atp<numAtoms;atp++)
	  if(atp-at>3)
	    { 
	      mapAt[m]=at+1;
	      mapAtp[m]=atp+1;
	      dist = getdist(CC[at],CC[atp]); 
	      distArray[m]=dist;
	      m++;
	    }
      //write out short distances with their atom pairs 
      m=0;
      int cntBadDist=0;
      // cout<<"for samp = "<<ns<<" distances less than disCut ="<<disCut<< " are"<<endl; //temp don't write
      int atBegin=0; 
      for(int at=atBegin;at<numAtoms-1;at++)
	for (int atp=at+1;atp<numAtoms;atp++)
	  {
	    if(distArray[m]>0 && distArray[m]<disCut)
	      if( mapAt[m]!=4 && mapAtp[m]!=8){ //exclude an initial first O to O distance issue
		//cout<<m+1<<"  "<<mapAt[m]<<" "<<mapAtp[m]<<" "<<distArray[m]<<endl; //temp don't write
		cntBadDist++;
	      }
	    m++;
	  }
      if(cntBadDist!=0)
	//cout<<"there are "<<cntBadDist<<" bad distances "<<endl; //temp don't write
	cout<<"for samp = "<<ns<<" there are "<<cntBadDist<<" bad distances "<<endl; 

      //get RG and EtoE (CA) and p(R) 
      if(cntBadDist==0)
	{
	  //incr numSampsKept at end
	  for(int i=0;i<3;i++)
	    COM[i]=0.0;
	  for(int i=0;i<3;i++)
	    {
	      for(int at=0;at<numAtoms;at++)
		COM[i] +=  CC[at][i];
	      COM[i] /= numAtoms;
	    }
	  double RGSq=0.0;
	  for(int at=0;at<numAtoms;at++)
	    for(int i=0;i<3;i++)
	      {
		RGSq += (CC[at][i]-COM[i])*(CC[at][i]-COM[i]);
	      }
	  RGSq /= numAtoms;
	  //cume RG
	  cumeRG += RGSq; 
	  
	  vector <double> eigs;
	  eigs = getEigvals(CC, numAtoms);
	  vector <double> Shapes(numShapes,0);
	  Shapes =  getShapes(eigs,numShapes);
	  //cume the shapes
	  transform(cumeShapes.begin(),cumeShapes.end(),Shapes.begin(),cumeShapes.begin(), plus<double>());

	  double EtoE=0.0;
	  int firstCA=1; //0-based indexing
	  int lastCA=4*(numDoub+1)-1; //0-based indexing
	  EtoE = getdist(CC[firstCA],CC[lastCA]); //CAs of first and last CA
	  
	  //cumulate HistoREtoEover samps
	  vector <double> HistoEtoE = getHistoEtoE(EtoE, numBinsEtoE, maxDistEtoE);
	  for(int bn=0;bn<numBinsEtoE;bn++)
	    cumeHistoEtoE[bn] += HistoEtoE[bn];

	  if(ns%writeInterval==0) 
	    ofRG<<ns<<"\t"<<sqrt(RGSq)<<"\t"<<Shapes[1]<<"\t"<<EtoE<<endl;

	  //get 14 HBs
	  int interval=13; //between firstAt and lastAt 
	  int firstAt=-1; //0 based index
	  int period=4;
	  int lastAt=0; //constructed
	  double distHB=0.0;
	  int firstRes=0; //for write purposes
	  double HBCut=3.5;
	  for (int nh=0;nh<nHB;nh++){
	    firstAt += period;
	    lastAt = firstAt + interval;
	    distHB=getdist(CC[firstAt],CC[lastAt]);
	    firstRes +=1;
	    if(distHB<HBCut) 
	      cntHB[nh]++;
	  }
	  //get 13 HBs
	  interval=9; //between firstAt and lastAt 
	  firstAt=-1; //0 based index
	  period=4;
	  lastAt=0; //constructed
	  distHB=0.0;
	  double HBCut13=3.5;
	  for (int nh=0;nh<nHB13;nh++){
	    firstAt += period;
	    lastAt = firstAt + interval;
	    distHB=getdist(CC[firstAt],CC[lastAt]);
	    firstRes +=1;
	    if(distHB<HBCut13) 
	      cntHB13[nh]++;
	  }

	  //cumulate histo for  PofR over kept samps	  
	  for(int m=0;m<numDist;m++)
	    distArrayAll[m]=-1;
	  dist=0.0;
	  int mc=0;
	  for(int at=0;at<numAtoms-1;at++)
	    for (int atp=at+1;atp<numAtoms;atp++)
	      { 
		dist = getdist(CC[at],CC[atp]); 
		distArrayAll[mc]=dist;
		mc++;
	      }
	  vector <double> Histo =  getPofR(distArrayAll,numDist, numBins, maxDist);
	  //cumulate Histo over samps
	  for(int bn=0;bn<numBins;bn++)
	    cumeHisto[bn] += Histo[bn];

	  //cumulate HistoRG over kept samps
	  vector <double> HistoRG = getHistoRG(RGSq, numBinsRG, maxDistRG);
	  for(int bn=0;bn<numBinsRG;bn++)
	    cumeHistoRG[bn] += HistoRG[bn];

	  //NMR section
	  //get karplus or BAX 
	  double temp=0.0;
	  double ang=0.0;
	  for(int nd=0;nd<numDoub;nd++){
	    ang=Phi[ns][nd];
	    jHNH[nd] += getJHNH(ang);
	    jHNHSQ[nd] += getJHNH(ang)*getJHNH(ang);
	  }

	  numSampsKept++; //for norm of cumulated things

	  if(ns%writeInterval==0) 
	    ofileKP<<"MODEL"<<"        "<<setiosflags(ios::left)
		   <<ns+1<<resetiosflags(ios::left)<<endl;
	  int cntAt=1;
	  //write pdbs to ofileKP
	  for(int nr=0;nr<numDoub+1;nr++) //+1 because of added init residue
	    {
	      for(int i=0;i<3;i++)
		atCC[i]=CC[cntAt-1][i];
	      if(ns%writeInterval==0) fwriteIt(ofileKP,atN,cntAt,nr+1,atCC,name);
	      cntAt++;
	      for(int i=0;i<3;i++)
		atCC[i]=CC[cntAt-1][i];
	      if(ns%writeInterval==0) fwriteIt(ofileKP,atCA,cntAt,nr+1,atCC,name);
	      cntAt++;
	      for(int i=0;i<3;i++)
		atCC[i]=CC[cntAt-1][i];
	      if(ns%writeInterval==0) fwriteIt(ofileKP,atC,cntAt,nr+1,atCC,name);
	      cntAt++;
	      for(int i=0;i<3;i++)
		atCC[i]=CC[cntAt-1][i];
	      if(ns%writeInterval==0) fwriteIt(ofileKP,atO,cntAt,nr+1,atCC,name);
	      cntAt++;
	    }
	  if(ns%writeInterval==0) ofileKP<<"TER"<<endl;
	  if(ns%writeInterval==0) ofileKP<<"ENDMDL"<<endl;

	} //end !bad 
    } //END SAMP LOOP
  
  cout<<endl<<"after numsamps"<<endl;
  //write cumeHistos
  for(int bn=0;bn<numBins;bn++)
    ofilePofR<<delBin2+delBin*bn<<"\t"<<cumeHisto[bn]/numSampsKept<<endl; 
  //get RG from PofR (RGsq=<R^2*P(R)>)
  double RG_from_PofR=0.0;
  double normHisto=0.0;
  for(int bn=0;bn<numBins;bn++){
    RG_from_PofR += cumeHisto[bn]*(delBin2+delBin*bn)*(delBin2+delBin*bn);
    normHisto += cumeHisto[bn];
  }
  cout<<"RG : RG_from_PorR"<<endl;
  cout<< sqrt(cumeRG/double(numSampsKept))<<" : "<<sqrt(RG_from_PofR/normHisto/2.0)<<endl;
  //write cumeHistoRG from cumeHistoRG
  for(int bn=0;bn<numBinsRG;bn++){
    ofilePofRG<<delBin2RG+delBinRG*bn<<"\t"<<cumeHistoRG[bn]/double(numSampsKept)<<endl;
  }
  for(int bn=0;bn<numBinsEtoE;bn++){
    ofilePofEtoE<<delBin2EtoE+delBinEtoE*bn<<"\t"<<cumeHistoEtoE[bn]/double(numSampsKept)<<endl;
  }

  //write RGsq delta and S
  cout<<"cumeShapes Averaged"<<endl;
  cout<<"RGsq delta S"<<endl;
  for(const auto& e:  cumeShapes)
    cout<<e/double(numSampsKept)<<" ";
  cout<<endl;
 
 //write HB to ofileHB
//HB14
  int intervalRes=4; //for write purposes
  int firstRes=0;
 ofileHB<<"percent\t14_HB"<<endl;
  for (int nh=0;nh<nHB;nh++){
    firstRes +=1;
    ofileHB<<firstRes<<"_"<<firstRes+intervalRes<<"\t"
	<<double(cntHB[nh])/double(numSampsKept)<<endl;
  }
  //HB13
  intervalRes=3; //for write purposes
  firstRes=0;
  ofileHB<<"percent\t13_HB"<<endl;
  for (int nh=0;nh<nHB13;nh++){
    firstRes +=1;
    ofileHB<<firstRes<<"_"<<firstRes+intervalRes<<"\t"
	<<double(cntHB13[nh])/double(numSampsKept)<<endl;
  }

  //write Jcoups to ofileNMR
  ofileNMR<<"nd\taveJHNH\tstdJHNH"<<endl;
  for(int nd=0;nd<numDoub;nd++){
    double A=jHNHSQ[nd]/numSampsKept;
    double B=(jHNH[nd]/numSampsKept)*(jHNH[nd]/numSampsKept);
    double var = A-B;
    if(A-B>1e-10)
      ofileNMR<<nd+1<<"\t"<<setprecision(2)<<jHNH[nd]/numSampsKept
	      <<"\t"<<sqrt(var)<<endl;
    else
      ofileNMR<<nd+1<<"\t"<</*setprecision(2)<<*/jHNH[nd]/numSampsKept
	      <<"\t"<<0.0<<endl; 
  }

  ofile.close();
  ofRG.close();
  ofileKP.close();
  ofilePofR.close();
  ofilePofRG.close();
  ofilePofEtoE.close();
  ofileHB.close();
  ofileNMR.close();
 
  clock_t end = clock();
  clock_t clockTicksTaken = end - start;
  double timeInSeconds = clockTicksTaken / (double) CLOCKS_PER_SEC;
  cout<<endl<<"For "<<numSampsKept<<" samps kept out of "<<numSamps<<
    " samps, took "<<timeInSeconds<<" seconds"<<endl;
  delete [] a;
  delete [] b;
  delete [] c;
  delete [] d;
  delete [] e;

  delete [] distArray;
  delete [] mapAt;
  delete [] mapAtp;
  delete [] COM;
  delete [] atCC;
  delete [] cntHB;

  for(int at=0;at<numAtoms;at++)
    delete [] CC[at];
  delete [] CC;
}

void writeIt(char* atom, int atomNum, int resNum, double* CC)
{
  cout<<"ATOM  "<<setw(5)<<atomNum<<"  "<<setw(4)<<setiosflags(ios::left)
      <<atom<<setw(3)<<"ALA"<<"  "<<setw(4)<<resetiosflags(ios::left)<<resNum
      <<"    "<<setw(8);
  for (int i=0;i<3;i++)
    cout<<CC[i]<<setprecision(3)<<fixed<<setw(8);
  cout<<setw(6)<<"1.00"<<endl;
}

//version with residue names
void fwriteIt(ofstream & ofile, char* atom, int atomNum, int resNum, double* CC, vector <string> name)
{
  ofile<<"ATOM  "<<setw(5)<<atomNum<<"  "<<setw(4)<<setiosflags(ios::left)
       <<atom<<setw(3)<<name[resNum-1]<<"  "<<setw(4)<<resetiosflags(ios::left)<<resNum
       <<"    "<<setw(8);
  for (int i=0;i<3;i++)
    ofile<<CC[i]<<setprecision(3)<<fixed<<setw(8);
  ofile<<setw(6)<<"1.00"<<" "<<10.00<<endl;
}

void refoldNRF(double *a, double *b, double *c, double* d, double bondLength, double bondAngle, double torsionAngle,double bondLengthn1)

{
  //IN: atom a,b,c CCs; bond CD (bondLength) , bond angle BCD (bondAngle) 
  //, torsion ABCD (torsionAngle), bondLengthn1 not used in this method.
  //OUT: CCs of atom d. 
#define rad .01745329252

  int ii;
  double D, length;
  float a_bx, a_by, a_bz;
  float cTD, cPsTD,sbA,sPsTD;
  float Mxx, Mxy, Mxz, Myx, Myy, Myz, Mzx, Mzy, Mzz;
  
  bondAngle = (bondAngle-180); 
  
  bondAngle*=rad;
  torsionAngle*=rad;

  sbA=sin(bondAngle);
  cTD =   bondLength * cos(bondAngle);

  a_bx = (b[0] - a[0]);
  a_by = (b[1] - a[1]);
  a_bz = (b[2] - a[2]);
  
  Mxx = (c[0] - b[0]) ;
  Mxy = (c[1] - b[1]) ;
  Mxz = (c[2] - b[2]) ;
  D = sqrt((Mxx)*(Mxx) + (Mxy)*(Mxy) + (Mxz)*(Mxz));
  Mxx/=D;
  Mxy/=D;
  Mxz/=D;

  cPsTD = bondLength * cos(torsionAngle) * sbA;
  sPsTD = bondLength * sin(torsionAngle) * sbA;

  Mzx = Mxy*a_bz - Mxz*a_by;
  Mzy = Mxz*a_bx - Mxx*a_bz;
  Mzz = Mxx*a_by - Mxy*a_bx;
  length = sqrt((Mzx * Mzx) + (Mzy * Mzy) + (Mzz * Mzz));
  Mzx = (Mzx / length);
  Mzy = (Mzy / length);
  Mzz = (Mzz / length);

  Myx = Mzy*Mxz - Mzz*Mxy;
  Myy = Mzz*Mxx - Mzx*Mxz;
  Myz = Mzx*Mxy - Mzy*Mxx;
  length = sqrt((Myx * Myx) + (Myy * Myy) + (Myz * Myz));
  Myx = (Myx / length);
  Myy = (Myy / length);
  Myz = (Myz / length);

  d[0] = Mxx*cTD + Myx*cPsTD + Mzx*sPsTD + c[0];
  d[1] = Mxy*cTD + Myy*cPsTD + Mzy*sPsTD + c[1];
  d[2] = Mxz*cTD + Myz*cPsTD + Mzz*sPsTD + c[2];
  
}

double getdist(double *a, double *b)
{
  double dist=0.0;
  for(int i=0;i<3;i++)
    {
      dist += (a[i]-b[i])*(a[i]-b[i]);
    }
  return dist = sqrt(dist);
}

vector<double> getPofR(double* distArray, int numDist, int numBins, double maxDist){
  //histo stuff
  double minDist=0.0;
  double delBin =maxDist/numBins;
  double delBin2=delBin/2.0;
  int hisBin; //local bn variable
  double value;
  int hisCount=0;
  vector <double> Histo(numBins);
  for(int bn=0;bn<numBins;bn++)
    Histo[bn]=0.0;
  for(int m=0;m<numDist;m++){
    value = distArray[m];
    hisBin = int (floor((value-minDist)/delBin) ); 
    if(hisBin>=0 && hisBin<numBins)
      {
	Histo[hisBin] += 1.0;
	hisCount++;
      }
  }
  for(int bn=0;bn<numBins;bn++)
    Histo[bn] /= hisCount;
  return Histo;
}

void getCovMatrix(double** CC,  int numAtoms , double** covMatrix){ //in CC and numAtoms; out covMatrix 
  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
      covMatrix[i][j]=0.0;
  //get aves first
  VecDoub ave(3);
  for(int i=0;i<3;i++)
    ave[i]=0.0;
  for(int i=0;i<3;i++){
    for(int at=0;at<numAtoms;at++)
      ave[i] += CC[at][i];
    ave[i]  /= numAtoms;
  }
  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
      {
	for(int at=0;at<numAtoms;at++)
	  covMatrix[i][j] += (CC[at][i]-ave[i])*(CC[at][j]-ave[j]);
	covMatrix[i][j] /= numAtoms;
      }
}

vector <double> getEigvals(double **CC,int numAtoms){
  vector <double> eigVals;    
  //VecDoub eigVals(3);
  MatDoub covMatrix(3,3); 
  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
      covMatrix[i][j]=0.0;
  //get aves first
  VecDoub ave(3);
  for(int i=0;i<3;i++)
    ave[i]=0.0;
  for(int i=0;i<3;i++){
    for(int at=0;at<numAtoms;at++)
      ave[i] += CC[at][i];
    ave[i]  /= numAtoms;
  }
  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
      {
	for(int at=0;at<numAtoms;at++)
	  covMatrix[i][j] += (CC[at][i]-ave[i])*(CC[at][j]-ave[j]);
	covMatrix[i][j] /= numAtoms;
      }
  SVD ob(covMatrix);
  for(int i=0;i<3;i++)
    eigVals.push_back(ob.w[i]); //convert array to vector
  return eigVals;
}

vector <double> getShapes(vector <double> eigVals,int numShapes){
  vector <double> Shapes(numShapes,0);
  double temp=0.0;
  for(int i=0;i<3;i++)
    temp += eigVals[i];
  Shapes[0]=temp; //RGsq
  //delta: asphericity
  Shapes[1] = 1.0-3.0*(eigVals[0]*eigVals[1]+eigVals[1]*eigVals[2]+eigVals[2]*eigVals[0])
    /(Shapes[0]*Shapes[0]);
  // S: shape parameter
  double lambAve=Shapes[0]/3.0;
  double numerator= (eigVals[0]-lambAve)*(eigVals[1]-lambAve)*(eigVals[2]-lambAve);
  double denom = pow(Shapes[0],3);
  Shapes[2]=27.0*numerator/denom;
  return Shapes;
}

vector<double> getHistoRG(double RGsq, int nB, double maxD){
  //histo stuff
  double minDist=0.0;
  double delBin =maxD/nB;
  double delBin2=delBin/2.0;
  int hisBin; //local bn variable
  double value;
  int hisCount=0;
  vector <double> Histo(nB);
  for(int bn=0;bn<nB;bn++)
    Histo[bn]=0.0;
  value = sqrt(RGsq); //we are binning RG values 
  hisBin = int (floor((value-minDist)/delBin) ); 
  if(hisBin>=0 && hisBin<nB)
    {
      Histo[hisBin] += 1.0;
      hisCount++;
    }

  for(int bn=0;bn<nB;bn++)
    Histo[bn] /= hisCount;
  return Histo;
}

vector<double> getHistoEtoE(double EtoE, int nB, double maxD){
  //histo stuff
  double minDist=0.0;
  double delBin =maxD/nB;
  double delBin2=delBin/2.0;
  int hisBin; //local bn variable
  double value;
  int hisCount=0;
  vector <double> Histo(nB);
  for(int bn=0;bn<nB;bn++)
    Histo[bn]=0.0;
  value = EtoE; 
  hisBin = int (floor((value-minDist)/delBin) ); 
  if(hisBin>=0 && hisBin<nB)
    {
      Histo[hisBin] += 1.0;
      hisCount++;
    }

  for(int bn=0;bn<nB;bn++)
    Histo[bn] /= hisCount;
  return Histo;
}

/*
//original Karplus as in Rauscher
double getJHNH(double phi){
double jhn=6.51*cos(phi-60)*cos(phi-60) -1.76*cos(phi-60)+1.60;
return jhn;
}
*/
//With BAX correction
double getJHNH(double phi){
  double degToRad = 0.0174533;
  double jhn=6.98*cos( (phi-60)*degToRad )*cos( (phi-60)*degToRad ) -1.38*cos( (phi-60)*degToRad )+1.72;
  return jhn;
}

//convert fasta to three letter code
string oneToThree(string res_type)
{
  if(res_type.compare(0,1,"A")==0)
    {
      return res_type.replace(0,3,"ALA");
    }
  else if(res_type.compare(0,1,"R")==0)
    {
      return res_type.replace(0,3,"ARG");
    }
  else if(res_type.compare(0,1,"N")==0)
    {
      return  res_type.replace(0,3,"ASN");
    }
  else if(res_type.compare(0,1,"D")==0)
    {
      return  res_type.replace(0,3,"ASP");
    }
  else if(res_type.compare(0,1,"C")==0)
    {
      return  res_type.replace(0,3,"CYS");
    }
  else if(res_type.compare(0,1,"Q")==0)
    {
      return  res_type.replace(0,3,"GLN");
    }
  else if(res_type.compare(0,1,"E")==0)
    {
      return  res_type.replace(0,3,"GLU");
    }
  else if(res_type.compare(0,1,"G")==0)
    {
      return  res_type.replace(0,3,"GLY");
    }
  else if(res_type.compare(0,1,"H")==0)
    {
      return  res_type.replace(0,3,"HIS");
    }
  else if(res_type.compare(0,1,"I")==0)
    {
      return  res_type.replace(0,3,"ILE");
    }
  else if(res_type.compare(0,1,"L")==0)
    {
      return  res_type.replace(0,3,"LEU");
    }
  else if(res_type.compare(0,1,"K")==0)
    {
      return  res_type.replace(0,3,"LYS");
    }
  else if(res_type.compare(0,1,"M")==0)
    {
      return  res_type.replace(0,3,"MET");
    }
  else if(res_type.compare(0,1,"F")==0)
    {
      return  res_type.replace(0,3,"PHE");
    }
  else if(res_type.compare(0,1,"P")==0)
    {
      return  res_type.replace(0,3,"PRO");
    }
  else if(res_type.compare(0,1,"S")==0)
    {
      return  res_type.replace(0,3,"SER");
    }
  else if(res_type.compare(0,1,"T")==0)
    {
      return  res_type.replace(0,3,"THR");
    }
  else if(res_type.compare(0,1,"W")==0)
    {
      return  res_type.replace(0,3,"TRP");
    }
  else if(res_type.compare(0,1,"Y")==0)
    {
      return  res_type.replace(0,3,"TYR");
    }
  else if(res_type.compare(0,1,"V")==0)
    {
      return  res_type.replace(0,3,"VAL");
    }

}