trunk/madProps/chrisProps.cpp

// Bloated Analysis Program
// Shamefully written by Chris Fennell
// Updated: 4/17/02
// Rewritten: 5/15/02

#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <string>
#include <cstring>
#include <vector>

using namespace std;

class Props {
public:
  Props(string fileName, const char *splitter, vector<int> &pT, int delR);
  ~Props();
  void DoRMSDCorrelation();
  void DoVelCorrelation();
  void DoGofR();
  void DoDipoleCorr();
  void DoOrientCorr();
  void DoMagFileCorr();
  void DoCosineCorr();
private:
  string filer;
  const char *delimiter,*file;
  char *token;
  char inLine[1000],inValue[200], allList[200], velList[200];
  char smallList[200], largeList[200],gofrList[200],sepFiler[200];
  char atomNeon[3],atomX[3],atomO[3],atomH[3],atomSSD[4];
  int nAtoms, i, j, k, l, simTime, simStep, counter, coorCount, count;
  int blackHole, m, start, end, atomer, countSm, countBg, fineGrain;
  int bin, option, outTyper, counter1, counter2, counter3, counter4, counter5;
  //  int *pAtomIdent;
  vector<int> pAtomIdent;
  double result, rX1, rY1, rZ1, rX2, rY2, rZ2, vX1, vY1, vZ1, vX2, vY2, vZ2;
  double uX1, uY1, uZ1, uX2, uY2, uZ2, velDot, velDotNaught, dotDivide;
  double mSdist, sum, average, xsum, ysum, xxsum, xysum, boxLength;
  double densityStar, rijDist, rXij, rYij, rZij, deltaR, constant;
  double rLower, rUpper, nIdeal, rho, rhoA, rhoB, rhoO, rhoH, rhoX;
  double constantA, constantB, constantX, constantO, constantH, nIdealA;
  double nIdealB, nIdealO, nIdealH, valueOuts, earlyTime, lateTime;
  double earlyAvg, lateAvg, totalAvg, timeDiffer, crossTime;
  double boxx, boxy, boxz, dipoleXo, dipoleYo, dipoleZo;
  double dipoleXt, dipoleYt, dipoleZt, dipODotdipO, dipTDotdipO;
};

Props::Props(string fileName, const char *splitter, vector<int> &pT, int delR){
  filer = fileName;
  file = filer.c_str();
  delimiter = splitter;
  pAtomIdent = pT;
  fineGrain = delR;
  

  strcpy(atomNeon,"Ne");
  strcpy(atomX,"X");
  strcpy(atomO,"O");
  strcpy(atomH,"H");
  strcpy(atomSSD,"SSD");

  ifstream inputer(file);

  inputer.getline(inLine,999,'\n');
  token = strtok(inLine, delimiter);
  strcpy(inValue,token);
  nAtoms = atoi(inValue);

  counter = 0; counter1 = 0; counter2 = 0; counter3 = 0; counter4 = 0; 
  counter5 = 0;
 
  // Identify the particles we're dealing with
  inputer.getline(inLine,999,'\n');
  for (i = 0; i < nAtoms; i++) {
    inputer.getline(inLine,999,'\n');
    token = strtok(inLine, delimiter);
    strcpy(inValue,token);
    if (strcmp(atomNeon,inValue) == 0) {
      pAtomIdent.push_back(0);
      counter++;
    }
    else if (strcmp(atomX,inValue) == 0) {
      pAtomIdent.push_back(2);
      counter2++;
    }
    else if (strcmp(atomO,inValue) == 0) {
      pAtomIdent.push_back(3);
      counter3++;
    }
    else if (strcmp(atomH,inValue) == 0) {
      pAtomIdent.push_back(4);
      counter4++;
    }
    else if (strcmp(atomSSD, inValue) == 0) {
      pAtomIdent.push_back(5);
      counter5++;
    }
    else {
      pAtomIdent.push_back(1);
      counter1++;
    }
  }

  inputer.clear(); 
}

Props::~Props() {
}


void Props::DoRMSDCorrelation() {
  strcpy(allList,file);
  strcat(allList,"rcorr");

  counter = nAtoms;

  if (counter5 > 0)
    counter = nAtoms;
  else {
    nAtoms = nAtoms/4; //A quick hack for SSD water (4 atoms in each SSD)
    counter = nAtoms;
  }
  
  vector<double> time;
  vector<double> pxSpot;
  vector<double> pySpot;
  vector<double> pzSpot;
  vector<double> avg;
  
  cout << "Reading in Trajectory...\n";
  
  // Fill the vectors with the needed information from the .dump or .xyz file
  ifstream inputer(file);
  coorCount = 0;

  if (counter5 > 0) {
    inputer.getline(inLine,999,'\n');
    while (!inputer.eof()) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      strcpy(inValue,token);
      time.push_back(atof(inValue));
      for (i = 0; i < nAtoms; i++) {
        inputer.getline(inLine,999,'\n');
        token = strtok(inLine, delimiter);
        token = strtok(NULL, delimiter);
        strcpy(inValue,token);
        pxSpot.push_back(atof(inValue));
        token = strtok(NULL, delimiter);
        strcpy(inValue,token);
        pySpot.push_back(atof(inValue));
        token = strtok(NULL, delimiter);
        strcpy(inValue,token);
        pzSpot.push_back(atof(inValue));
      }
      coorCount++;
      inputer.getline(inLine,999,'\n');
    }
  }
  else {
    inputer.getline(inLine,999,'\n');
    while (!inputer.eof()) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      strcpy(inValue,token);
      time.push_back(atof(inValue));
      for (i = 0; i < nAtoms; i++) {
        inputer.getline(inLine,999,'\n');
        token = strtok(inLine, delimiter);
        token = strtok(NULL, delimiter);
        strcpy(inValue,token);
        pxSpot.push_back(atof(inValue));
        token = strtok(NULL, delimiter);
        strcpy(inValue,token);
        pySpot.push_back(atof(inValue));
        token = strtok(NULL, delimiter);
        strcpy(inValue,token);
        pzSpot.push_back(atof(inValue));
        inputer.getline(inLine,999,'\n');
        inputer.getline(inLine,999,'\n');
        inputer.getline(inLine,999,'\n');
      }
      coorCount++;
      inputer.getline(inLine,999,'\n');
    }
  }
  cout << "Performing RMSD Calculations...\n";
  
  // The main calculation loop...
  for (j = 0; j < coorCount-1; j++) {
    sum = 0;
    blackHole = 0;
    for (i = 0; (i + j) < coorCount; i++) {
      for (k = 0; k < counter; k++) {
        rX1 = pxSpot[k + (i * counter)];
        rY1 = pySpot[k + (i * counter)];
        rZ1 = pzSpot[k + (i * counter)];
        rX2 = pxSpot[k + (i * counter) + (j * counter)];
        rY2 = pySpot[k + (i * counter) + (j * counter)];
        rZ2 = pzSpot[k + (i * counter) + (j * counter)];
        
        mSdist = (pow((rX2-rX1),2) + pow((rY2-rY1),2)
                  + pow((rZ2-rZ1),2));
        
        sum += mSdist;
      }
      blackHole++;
    }
    average = sum / counter;
    average /= blackHole;
    avg.push_back(average);
  }
  
  // Write out the results...
  ofstream deltaOut(allList);
  deltaOut << "#dt\tRMSD\n";
  for (i = 0; i < coorCount-1; i++)
    deltaOut << time[i]-time[0] << "\t" << avg[i] << "\n";
  
  inputer.clear();
}


void Props::DoVelCorrelation() {
  strcpy(velList,file);
  strcat(velList,"vcorr");

  counter = nAtoms;

  vector<double> time;
  vector<double> pxVel;
  vector<double> pyVel;
  vector<double> pzVel;
  vector<double> avg;

  ifstream inputer(file);

  if (counter5 == 0) {
    cout << "Note: VCorr should only be run on a .dump file\n";
    return;
  }

  //Load all velocities into the appropriate vectors declared above.
  cout << "Reading in Trajectory...\n";
  coorCount = 0;
  inputer.getline(inLine,999,'\n');
  while (!inputer.eof()) {
    inputer.getline(inLine,999,'\n');
    token = strtok(inLine, delimiter);
    strcpy(inValue,token);
    time.push_back(atof(inValue));
    for (i = 0; i < nAtoms; i++) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pxVel.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pyVel.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pzVel.push_back(atof(inValue));
     }
    coorCount++;
    inputer.getline(inLine,999,'\n');
  }
    
  // Do Velocity Correlation Calculations...
  cout << "Performing VCorr Calculations...\n";
  for (j = 0; j < coorCount-1; j++) {
    sum = 0;
    blackHole = 0;
    for (i = 0; (i + j) < coorCount-1; i++) {
      for (k = 0; k < nAtoms; k++) {
        vX1 = pxVel[(k + (i * nAtoms))];
        vY1 = pyVel[(k + (i * nAtoms))];
        vZ1 = pzVel[(k + (i * nAtoms))];
        vX2 = pxVel[k + (i * nAtoms) + (j * nAtoms)];
        vY2 = pyVel[k + (i * nAtoms) + (j * nAtoms)];
        vZ2 = pzVel[k + (i * nAtoms) + (j * nAtoms)];
        velDot = vX1*vX2 + vY1*vY2 + vZ1*vZ2;
        velDotNaught = vX1*vX1 + vY1*vY1 + vZ1*vZ1;
        dotDivide = velDot / velDotNaught;
        
        sum += dotDivide;
      }
      blackHole++;
    }
    average = sum / nAtoms;
    average /= blackHole;
    avg.push_back(average);
  }

  // Write out the results...
  ofstream deltaOut(velList);
  deltaOut << "#dt\tVCorr\n";
  for (i=0; i < coorCount-1; i++)
    deltaOut << (time[i]-time[0]) << "\t" << avg[i] << "\n";  
  inputer.clear();
}


void Props::DoGofR() {
  strcpy(gofrList,file);
  strcat(gofrList,"gofr");

  count = nAtoms/4;

  vector<double> pxSpot;
  vector<double> pySpot;
  vector<double> pzSpot;
  //  vector<int> pAtomIdent;
  int *histogram = new int[fineGrain];
  int *OOhistogram = new int[fineGrain];
  int *OHhistogram = new int[fineGrain];
  int *HHhistogram = new int[fineGrain];
  double *prValues= new double[fineGrain];
  double *OOGofR = new double[fineGrain];
  double *HHGofR = new double[fineGrain];
  double *OHGofR = new double[fineGrain];
  double *GofR = new double[fineGrain];

  ifstream inputer(file);
  
  //Load all positions into the appropriate vectors declared above.
  coorCount = 0;
  cout << "Reading in Trajectory...\n";
  inputer.getline(inLine,999,'\n');
  inputer.getline(inLine,999,'\n');
  token = strtok(inLine,delimiter);
  token = strtok(NULL,delimiter);
  strcpy(inValue,token);
  boxx = atof(inValue);
  token = strtok(NULL,delimiter);
  strcpy(inValue,token);
  boxy = atof(inValue);
  token = strtok(NULL,delimiter);
  strcpy(inValue,token);
  boxz = atof(inValue);
  while (!inputer.eof()) {
    for (i = 0; i < nAtoms; i++) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pxSpot.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pySpot.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pzSpot.push_back(atof(inValue));
    }
    coorCount++;
    inputer.getline(inLine,999,'\n');
    inputer.getline(inLine,999,'\n');
  }

  boxLength = boxy;
  deltaR = (boxLength/2)/fineGrain;  
  rho = (nAtoms / pow(boxLength, 3));
  constant = (4.0 * M_PI * rho) / 3.0; 
  rhoO = ((nAtoms/4) / pow(boxLength, 3));
  rhoH = ((nAtoms/2) / pow(boxLength, 3));
  rhoX = ((nAtoms) / pow(boxLength, 3));
  constantX = (4.0 * M_PI * rhoX) / 3.0;
  constantO = (4.0 * M_PI * rhoO) / 3.0;
  constantH = (4.0 * M_PI * rhoH) / 3.0;
  
  //Zero out all of the g(r)'s and histogram bins
  for (i = 0; i < fineGrain; i++) {
    prValues[i] = 0;
    OOGofR[i] = 0;
    HHGofR[i] = 0;
    OHGofR[i] = 0;
    GofR[i] = 0;
    histogram[i] = 0;
    OOhistogram[i] = 0;
    HHhistogram[i] = 0;
    OHhistogram[i] = 0;
  }
  
  cout << "Performing G(r) Calculations...\n";
  
  // The main calculation loops...
  for (j = 0; j < coorCount; j++) { 
    for (i = 0; i < nAtoms-1; i++) {
      rX1 = pxSpot[(i + (j * nAtoms))];
      rY1 = pySpot[(i + (j * nAtoms))];
      rZ1 = pzSpot[(i + (j * nAtoms))];
      for (k = 1; (k + i) < nAtoms; k++) {
        rX2 = pxSpot[(k + i + (j * nAtoms))];
        rY2 = pySpot[(k + i + (j * nAtoms))];
        rZ2 = pzSpot[(k + i + (j * nAtoms))];
        
        rXij = (rX1-rX2);
        rYij = (rY1-rY2);
        rZij = (rZ1-rZ2);
        
        //Do minimum image filter...
        rXij -= boxLength*copysign(1.0,rXij)*floor(fabs(rXij/boxLength)+0.5);
        rYij -= boxLength*copysign(1.0,rYij)*floor(fabs(rYij/boxLength)+0.5);
        rZij -= boxLength*copysign(1.0,rZij)*floor(fabs(rZij/boxLength)+0.5);
        rijDist = sqrt(rXij*rXij + rYij*rYij + rZij*rZij);

        //Add to appropriate histogram bin...
        bin = (int)(rijDist / deltaR);
        if (bin < fineGrain) {
          if (pAtomIdent[i] == 3 && pAtomIdent[(k+i)] == 3) {
            OOhistogram[bin] = OOhistogram[bin]+2;
          }
          else if(pAtomIdent[i] == 4 && pAtomIdent[(k+i)] == 4) {
            HHhistogram[bin] = HHhistogram[bin]+2;
          }
          else if(pAtomIdent[i] == 3 && pAtomIdent[(k+i)] == 4) {
            OHhistogram[bin] = OHhistogram[bin]+2;
          }
          else if(pAtomIdent[i] == 4 && pAtomIdent[(k+i)] == 3) {
            OHhistogram[bin] = OHhistogram[bin]+2;
          }
          else if (counter5 > 0) {
            histogram[bin] = histogram[bin]+2;
          }
        }       
      }
    }
    
    // Calculate the G(r) values...
    for (i = 0; i < fineGrain; i++) {
      rLower = i * deltaR;
      rUpper = rLower + deltaR;
      nIdeal = constant * (pow(rUpper, 3) - pow(rLower, 3));
      nIdealO = constantO * (pow(rUpper, 3) - pow(rLower, 3));
      nIdealH = constantH * (pow(rUpper, 3) - pow(rLower, 3));
      
      GofR[i] = histogram[i] / (coorCount * nAtoms * nIdeal);
      OOGofR[i] = OOhistogram[i] / (coorCount * (nAtoms/4) * nIdealO);
      HHGofR[i] = HHhistogram[i] / (coorCount * (nAtoms/2) * nIdealH);
      //Check if this one works...
      OHGofR[i] = OHhistogram[i] / (coorCount * (nAtoms/4) * nIdealO * 4);
      
      prValues[i] = rLower + (deltaR/2);
    }
  }  

  // Write out the results...
  strcpy(sepFiler, gofrList);
  ofstream deltaOut(gofrList);
  if (counter5 > 0) {
    deltaOut << "#rDist\tSSDG(r)\n";
    for (i=0; i<fineGrain; i++)
      deltaOut << prValues[i] << "\t" << GofR[i] << "\n";
  }
  else {
    deltaOut << "#rDist\tOOG(r)\tHHG(r)\tOHG(r)\n";
    for (i=0; i < fineGrain; i++)
      deltaOut << prValues[i] << "\t" << OOGofR[i] << "   \t" 
               << HHGofR[i] << "   \t" << OHGofR[i] << "\n";
  }
  inputer.clear();
  delete[] OOhistogram; delete[] OHhistogram; delete[] HHhistogram; 
  delete[] prValues; delete[] GofR; delete[] OOGofR; 
  delete[] OHGofR; delete[] HHGofR;
}


void Props::DoDipoleCorr() {
  strcpy(allList,file);
  strcat(allList,"dcorr");

  counter = nAtoms;

  if (counter5 > 0) {
    cout << "Not yet working for this file type.  Use an .xyz file.\n";
    return;
  }
  else {
    nAtoms = nAtoms/4; //A quick hack for SSD water (4 atoms in each SSD)
    counter = nAtoms;
  }

  vector<double> time;
  vector<double> pUnitX;
  vector<double> pUnitY;
  vector<double> pUnitZ;
  double dot;
  double zeroval;
  int nloops;
  
  ifstream inputer(file);

  // Load all positions and unit vectors into the appropriate 
  // vectors declared above.
  coorCount = 0;
  cout << "Reading in Trajectory...\n";
  inputer.getline(inLine,999,'\n');
  while (!inputer.eof()) {
    inputer.getline(inLine,999,'\n');
    token = strtok(inLine, delimiter);
    strcpy(inValue,token);
    time.push_back(atof(inValue));
    for (i = 0; i < nAtoms; i++) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitX.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitY.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitZ.push_back(atof(inValue));
      inputer.getline(inLine,999,'\n');
      inputer.getline(inLine,999,'\n');
      inputer.getline(inLine,999,'\n');
    }
    coorCount++;
    inputer.getline(inLine,999,'\n');
  }

  double *mx = new double[coorCount];
  double *my = new double[coorCount];
  double *mz = new double[coorCount];
  double *avg = new double[coorCount];
  
  cout << "Performing DCorr Calculations...\n";

  // Compute total magnetizations
  for (j = 0; j < coorCount; j++) {
    mx[j] = 0;
    my[j] = 0;
    mz[j] = 0;
    for (k = 0; k < counter; k++) {
      mx[j] += pUnitX[k + (j*counter)];
      my[j] += pUnitY[k + (j*counter)];
      mz[j] += pUnitZ[k + (j*counter)];
    }
  }

  // Perform time correlations
  for (j = 0; j < coorCount-1; j++) {
    nloops = 0;
    dot = 0.0;
    
    for (i = 0 ; (i+j)<coorCount-1; i++) {
      
      dot += mx[i]*mx[i+j] + my[i]*my[i+j] + mz[i]*mz[i+j];
      nloops++;
    }
    
    avg[j] = dot / (double)nloops;
    
    if (j == 0) {
      zeroval = avg[j];
    }
    
     avg[j] = avg[j] / zeroval;
  }

  // Write out the results...
  ofstream deltaOut(allList);
  deltaOut << "#dt\tDipole Corr\n";
  for (i = 0; i < time.size()-1; i++) 
    deltaOut << (time[i]-time[0]) << "\t" << avg[i] << "\n";

  inputer.clear();
  delete[] mx; delete[] my; delete[] mz; delete[] avg;
}

void Props::DoOrientCorr() {
  strcpy(allList,file);
  strcat(allList,"mucorr");

  counter = nAtoms;

  if (counter5 > 0) {
    cout << "Not yet working for this file type.  Use an .xyz file.\n";
    return;
  }
  else {
    nAtoms = nAtoms/4; //A quick hack for SSD water (4 atoms in each SSD)
    counter = nAtoms;
  }
  
  vector<double> time;
  vector<double> pUnitX;
  vector<double> pUnitY;
  vector<double> pUnitZ;
  double dipDot;
  double dipDotNaught;
  double legendreNum;

  ifstream inputer(file);
  
  // Load all positions and unit vectors into the appropriate 
  // vectors declared above.
  coorCount = 0;
  cout << "Reading in Trajectory...\n";
  inputer.getline(inLine,999,'\n');
  while (!inputer.eof()) {
    inputer.getline(inLine,999,'\n');
    token = strtok(inLine, delimiter);
    strcpy(inValue,token);
    time.push_back(atof(inValue));
    for (i = 0; i < counter; i++) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitX.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitY.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitZ.push_back(atof(inValue));
      inputer.getline(inLine,999,'\n'); // Jump over O and H atom lines
      inputer.getline(inLine,999,'\n'); // "
      inputer.getline(inLine,999,'\n'); // "
    }
    coorCount++;
    inputer.getline(inLine,999,'\n');
  }

  double *pC1 = new double[coorCount];
  double *pC2 = new double[coorCount];
  double *dipDotSum = new double[coorCount];
  double *dipDotSum2 = new double[coorCount];

  cout << "Performing muCorr Calculations...\n";
  // The main calculation loops...
  for (j = 0; j < coorCount-1; j++) {
    dipDotSum[j] = 0;
    dipDotSum2[j] = 0;
    blackHole = 0;
    for (i = 0; (i + j) < coorCount-1; i++) {
      for (k = 0; k < counter; k++) {
        uX1 = pUnitX[k + (i * counter)];
        uY1 = pUnitY[k + (i * counter)];
        uZ1 = pUnitZ[k + (i * counter)];
        uX2 = pUnitX[k + (i * counter) + (j * counter)];
        uY2 = pUnitY[k + (i * counter) + (j * counter)];
        uZ2 = pUnitZ[k + (i * counter) + (j * counter)];
        
        dipDot = uX1*uX2 + uY1*uY2 + uZ1*uZ2;
        legendreNum = 1.5 * dipDot * dipDot - 0.5;
        
        dipDotSum[j] += dipDot;
        dipDotSum2[j] += legendreNum;
      }
      blackHole++;
    }
    
    dipDotSum[j] = dipDotSum[j] / (blackHole*counter);
    dipDotSum2[j] = dipDotSum2[j] / (blackHole*counter);
  }
  
  for (j = 0; j < coorCount-1; j++) {
    pC1[j] = dipDotSum[j] / dipDotSum[0];
    pC2[j] = dipDotSum2[j] / dipDotSum2[0];
  }
  
  
  ofstream deltaOut(allList);
  deltaOut << "#dt\tC1(t)\tC2(t)\n";
  for (i = 0; i < coorCount-1; i++)
    deltaOut << (time[i]-time[0]) << "\t" << pC1[i] << "\t" << pC2[i] << "\n";
  inputer.clear();
  
  delete[] pC1; delete[] pC2;
}


void Props::DoMagFileCorr() {
  strcpy(allList,file);
  strcat(allList,"magcorr");

  vector<double> time;
  vector<double> pUnitX;
  vector<double> pUnitY;
  vector<double> pUnitZ;
  double dot;
  double zeroval;
  int nloops;
  
  ifstream inputer(file);

  // Load all positions and unit vectors into the appropriate 
  // vectors declared above.
  coorCount = 0;
  cout << "Reading in Trajectory...\n";
  inputer.getline(inLine,999,'\n');
  while (!inputer.eof()) {
    token = strtok(inLine, delimiter);
    strcpy(inValue,token);
    time.push_back(atof(inValue));
    token = strtok(NULL, delimiter);
    strcpy(inValue,token);
    pUnitX.push_back(atof(inValue));
    token = strtok(NULL, delimiter);
    strcpy(inValue,token);
    pUnitY.push_back(atof(inValue));
    token = strtok(NULL, delimiter);
    strcpy(inValue,token);
    pUnitZ.push_back(atof(inValue));
    coorCount++;
    inputer.getline(inLine,999,'\n');
  }

  double *avg = new double[time.size()];
  
  cout << "Performing magCorr Calculations...\n";

  // Perform time correlations
  for (j = 0; j < coorCount-1; j++) {
    nloops = 0;
    dot = 0.0;
    
    for (i = 0 ; (i+j)<coorCount-1; i++) {      
      dot += pUnitX[i]*pUnitX[i+j] + pUnitY[i]*pUnitY[i+j] 
        + pUnitZ[i]*pUnitZ[i+j];
      nloops++;
    }
    
    avg[j] = dot / (double)nloops;
    
    if (j == 0) {
      zeroval = avg[j];
    }    
     avg[j] = avg[j] / zeroval;
  }

  // Write out the results...
  ofstream deltaOut(allList);
  deltaOut << "#dt\tMag Corr\n";
  for (i = 0; i < time.size()-1; i++) 
    deltaOut << (time[i]-time[0]) << "\t" << avg[i] << "\n";

  inputer.clear();
  delete[] avg;
}

void Props::DoCosineCorr() {
  strcpy(gofrList,file);
  strcat(gofrList,"coscorr");

  counter = nAtoms/4;

  vector<double> pxSpot;
  vector<double> pySpot;
  vector<double> pzSpot;
  vector<double> pUnitX;
  vector<double> pUnitY;
  vector<double> pUnitZ;
  double dot;
  int *histogram = new int[fineGrain];
  double *prValues= new double[fineGrain];
  double *avgCos = new double[fineGrain];

  if (counter5 > 0) {
    cout << "File type currently unsupported.  Please use a .xyz file.\n";
    return;
  }

  ifstream inputer(file);
  
  //Load all positions into the appropriate vectors declared above.
  coorCount = 0;
  cout << "Reading in Trajectory...\n";
  inputer.getline(inLine,999,'\n');
  inputer.getline(inLine,999,'\n');
  token = strtok(inLine,delimiter);
  token = strtok(NULL,delimiter);
  strcpy(inValue,token);
  boxx = atof(inValue);
  token = strtok(NULL,delimiter);
  strcpy(inValue,token);
  boxy = atof(inValue);
  token = strtok(NULL,delimiter);
  strcpy(inValue,token);
  boxz = atof(inValue);
  while (!inputer.eof()) {
    for (i = 0; i < counter; i++) {
      inputer.getline(inLine,999,'\n');
      token = strtok(inLine, delimiter);
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pxSpot.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pySpot.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pzSpot.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitX.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitY.push_back(atof(inValue));
      token = strtok(NULL, delimiter);
      strcpy(inValue,token);
      pUnitZ.push_back(atof(inValue));
      inputer.getline(inLine,999,'\n');
      inputer.getline(inLine,999,'\n');
      inputer.getline(inLine,999,'\n');
    }
    coorCount++;
    inputer.getline(inLine,999,'\n');
    inputer.getline(inLine,999,'\n');
  }

  boxLength = boxy;
  deltaR = (boxLength/2)/fineGrain;  
  
  //Zero out all of the g(r)'s and histogram bins
  for (i = 0; i < fineGrain; i++) {
    prValues[i] = 0;
    histogram[i] = 0;
    avgCos[i] = 0;
  }
  
  cout << "Performing cosCorr Calculations...\n";
  // Main sorting loops...
  for (j = 0; j < coorCount; j++) { 
    for (i = 0; i < counter-1; i++) {
      rX1 = pxSpot[(i + (j * counter))];
      rY1 = pySpot[(i + (j * counter))];
      rZ1 = pzSpot[(i + (j * counter))];
      uX1 = pUnitX[(i + (j * counter))];
      uY1 = pUnitY[(i + (j * counter))];
      uZ1 = pUnitZ[(i + (j * counter))];
      for (k = 1; (k + i) < counter; k++) {
        rX2 = pxSpot[(k + i + (j * counter))];
        rY2 = pySpot[(k + i + (j * counter))];
        rZ2 = pzSpot[(k + i + (j * counter))];
        
        rXij = (rX1-rX2);
        rYij = (rY1-rY2);
        rZij = (rZ1-rZ2);
        
        //Do minimum image filter...
        rXij -= boxLength*copysign(1.0,rXij)*floor(fabs(rXij/boxLength)+0.5);
        rYij -= boxLength*copysign(1.0,rYij)*floor(fabs(rYij/boxLength)+0.5);
        rZij -= boxLength*copysign(1.0,rZij)*floor(fabs(rZij/boxLength)+0.5);
        rijDist = sqrt(rXij*rXij + rYij*rYij + rZij*rZij);
        
        //Now do the average cosine b/t dipoles
        uX2 = pUnitX[(k + i + (j * counter))];
        uY2 = pUnitY[(k + i + (j * counter))];
        uZ2 = pUnitZ[(k + i + (j * counter))];
        
        dot = uX1*uX2 + uY1*uY2 + uZ1*uZ2;

        //Add to appropriate histogram bin...
        bin = (int)(rijDist / deltaR);
        if (bin < fineGrain) {
            histogram[bin] = histogram[bin] + 1;
            avgCos[bin] = avgCos[bin] + dot;
        }       
      }
    }
  }  

  for (i = 0; i < fineGrain; i++) {
    if (histogram[i] != 0)
      avgCos[i] = avgCos[i] / histogram[i];
  }

  for (i = 0; i < fineGrain; i++) {
    rLower = i * deltaR;
    prValues[i] = rLower + (deltaR/2);
  }
  
  // Write out the results...
  strcpy(sepFiler, gofrList);
  ofstream deltaOut(gofrList);
  deltaOut << "#Radius\tcosCorr\n";
  for (i=0; i<fineGrain; i++)
      deltaOut << prValues[i] << "\t" << avgCos[i] << "\n";
  
  inputer.clear();
  delete[] prValues; delete[] histogram; delete[] avgCos;
}


int main(int argc, char *argv[]) {
  char choice[200], inplace[200], atomChar[10], inLine[1000], inValue[200];
  char temp[4];
  char fracString[4] = "";
  char *token;
  const char *delimit = " \t\n";
  int selection,nParticles,i,grain;
  vector<int> pAtomType;
  string strungName;

  // Command line input format test
  if (argc != 3) {
    cout 
      << "Usage: " << argv[0] << " <file name> <correlation choice>\n"
      << "\n"
      << "Correlation choice options:\n"
      << "\n"
      << "  rcorr - Root Mean-Squared Displacement\n"
      << "  vcorr - Velocity Autocorrelation Function\n"
      << "  gofr - Pair Correlation Function\n"
      << "  dcorr - Collective Dipole Correlation Function\n"
      << "  mucorr - Self Dipolar Orientational Correlation Function (C1 & C2)\n"
      << "  magfile - Collective Dipole Correlation Function\n"
      << "  coscorr - Distance-Dependent Average Cosine Correlation\n";
    return 0;
  }
 
  ifstream prayer(argv[1]);

  // Make sure the file exists
  if (!prayer) { 
    cout << "Unable to open " << argv[1] << " for reading.\n";
    return 0;
  }

  prayer.close();

  // Build a filename string
  strungName = argv[1];
  grain = 500; //detail of g(r) curve (# of bins) - increase for smoother curve

  // Assign the correlation to perform
  if (!strcmp(argv[2], "rcorr"))
    selection = 1;
  else if (!strcmp(argv[2], "vcorr"))
    selection = 2;
  else if (!strcmp(argv[2], "gofr"))
    selection = 3;
  else if (!strcmp(argv[2], "dcorr"))
    selection = 4;
  else if (!strcmp(argv[2], "mucorr"))
    selection = 5;
  else if (!strcmp(argv[2], "magfile"))
    selection = 6;
  else if (!strcmp(argv[2], "coscorr"))
    selection = 7;
  else {
    cout 
      << "Not a valid correlation. Choose one of the following:\n"
      << "\n"
      << "  rcorr - Root Mean-Squared Displacement\n"
      << "  vcorr - Velocity Autocorrelation Function\n"
      << "  gofr - Pair Correlation Function\n"
      << "  dcorr - Collective Dipole Correlation Function\n"
      << "  mucorr - Self Dipolar Orientational Correlation Function (C1 & C2)\n"
      << "  magfile - Collective Dipole Correlation Function\n"
      << "  coscorr - Distance-Dependent Average Cosine Correlation\n";
    return 0;
  }

  // Call up correlation calculations
  
  Props analyzer(strungName,delimit,pAtomType,grain);
  
  switch( selection ){
    
  case 1:
    analyzer.DoRMSDCorrelation();
    break;

  case 2:
    analyzer.DoVelCorrelation();
    break;

  case 3:
    analyzer.DoGofR();
    break;
    
  case 4:
    analyzer.DoDipoleCorr();
    break;

  case 5:
    analyzer.DoOrientCorr();
    break;
    
  case 6:
    analyzer.DoMagFileCorr();
    break;
    
  case 7:
    analyzer.DoCosineCorr();
    break;
   
  default:
    cout << "Error in understanding the selection\n";
  }
  
  return 0;
}
Revision:	39
Committed:	Fri Jul 19 01:37:38 2002 UTC (21 years, 11 months ago) by mmeineke
File size:	27754 byte(s)
Log Message:	This commit was generated by cvs2svn to compensate for changes in r38, which included commits to RCS files with non-trunk default branches.