--- trunk/tcProps/rmsd.c	2004/02/27 20:39:17	1072
+++ trunk/tcProps/rmsd.c	2004/02/28 16:45:57	1073
@@ -17,10 +17,19 @@ strcut frameStruct* myFrames;
   double Hmat[3][3];
 };
 
-strcut frameStruct* myFrames;
+struct binStruct{
+  double rmsd;
+  int nValues;
+};
 
-void rmsd(double startTime, char* outPrefix ){
+struct binStruct rmsdCorr[RMSDBINS];
 
+struct frameStruct* myFrames;
+
+void calcRMSDpair( int bin, int i, int j, enum atomNames rType );
+
+void rmsd( enum atomNames rType, double startTime, char* outPrefix ){
+
   // list of 'a priori' constants
 
   const int nLipAtoms = NL_ATOMS;
@@ -32,12 +41,20 @@ void rmsd(double startTime, char* outPrefix ){
   // variables
 
   char outName[500];
+  
   FILE* outFile;
+  
   int i, j, k;
   int startFrame, corrFrames, framesFinished;
   int startFound, percentComplete;
-  
+  int index;
+  int nCounts;
+  int bin;
+
   double Hmat[3][3];
+  double diffTime;
+  double rmsdDt;
+  double timeOut, outVal;
 
   framesFinished = 0;
   
@@ -64,3 +81,191 @@ void rmsd(double startTime, char* outPrefix ){
   corrFrames = nFrames - startFrame;
   myFrames = (struct frameStruct*)calloc(corrFrames, 
 					 sizeof(struct frameStruct));
+
+  index=0;
+  for(i=startFrame;i<nFrames;i++){
+    
+    if( index >= corrFrames ){
+      fprintf( stderr,
+	       "Error, number of frames, %d, exceeds corrFrames, %d.\n",
+	       index,
+	       corrFrames );
+    }
+   
+    myFrames[index].time = frameTimes[i];
+    readFrame(i, myFrames[index].atoms, myFrames[index].Hmat );
+  }
+
+  // initialize the counts and the correlation
+
+  nCounts = 0;
+  for(i=0;i<(corrFrames-1);i++)
+    for(j=i+1;j<corrFrames;j++)
+      nCounts++;
+ 
+  for(i=0;i<RMSDBINS;i++){
+    rmsdCorr[i].rmsd = 0.0;
+    rmsdCorr[i].nValues = 0;
+  }
+
+  rmsdDt = (myFrames[corrFrames-1].time - myFrames[0].time) / RMSDBINS;
+   
+  // perform the calculation
+
+  index = 0;
+  for(i=0;i<(corrFrames-1);i++){
+    for(j=i+1;j<corrFrames;j++){
+
+      diffTime = myFrames[j].time - myFrames[i].time;
+      if( diffTime > 0.0 ){
+	
+	percentComplete = 
+	  (int)( 100.0 * (double)index / (double) nCounts );
+	
+	fprintf( stdout,
+		 "\rRMSD corr %3d%% complete.",
+		 percentComplete );
+	fflush( stdout );
+	
+	bin = (int)(diffTime / rmsdDt);
+	if( bin < RMSDBINS)
+	  calcRMSDpair( bin, i, j, rType );
+      }
+
+      index++;
+    }
+  }
+
+  // print out the correlation
+
+  sprintf( outName, "%s.rmsd", outPrefix );
+  outFile = fopen( outName, "w" );
+
+  fprintf( outFile,
+	   "#time\trmsd\n");
+
+  for(i=0;i<RMSDBINS;i++){
+    
+    if(rmsdCorr[i].nValues){
+      
+      timeOut = (2.0*(double)i+1.0)*rmsdDt*0.5;
+      outVal = rmsdCorr[i].rmsd / (double)rmsdCorr[i].nValues;
+
+      fprintf( outFile,
+	       "%6G\t%6G\n",
+	       timeOut, outVal );
+    }
+  }
+  
+  fflush(outFile);
+  fclose(outFile);
+
+  percentComplete = 
+    (int)( 100.0 * (double)index / (double) nCounts );
+  
+  fprintf( stdout,
+	   "\rRMSD corr %3d%% complete.",
+	   percentComplete );
+  fflush( stdout );
+
+  free( myFrames );
+  myFrames = NULL;
+}
+
+
+void calcRMSDpair( int bin, int frameI, int frameJ, enum atomNames rType ){
+
+  struct atomCoord* atomsI;
+  struct atomCoord* atomsJ;
+  
+  int i,j,k,l,m,n;
+  int nAccums;
+
+  double d[3];
+  double dSqr;
+  double comI[3];
+  double comJ[3];
+  double totMass;
+  double accum;
+  
+  
+  
+  // list of 'a priori' constants
+
+  const int nLipAtoms = NL_ATOMS;
+  const int nBonds    = NBONDS;
+  const int nLipids   = NLIPIDS;
+  const int nSSD      = NSSD;
+  const int nAtoms    = nLipAtoms * nLipids + nSSD;
+  
+  
+  atomsI = myFrames[frameI].atoms;
+  atomsJ = myFrames[frameJ].atoms;
+
+  if( rType != COM ){
+    
+    nAccums = 0;
+    for(i=0;i<nAtoms;i++){
+      
+      if (atomsI[i].type == rType){
+	
+	for(j=0;j<3;j++)
+	  d[j] = atomsJ[i].pos[j] - atomsI[i].pos[j];
+	
+	dSqr = 0.0;
+	for(j=0;j<3;j++)
+	  dSqr += d[j] * d[j];
+	
+	accum += sqrt(dSqr);
+	nAccums++;
+      }    
+    }
+  }
+  else{
+    
+    nAccums = 0;
+    for(i=0;i<nLipids;i++){
+      
+      k = i*nLipAtoms;
+      
+      // com calc
+    
+      totMass = 0.0;
+      for(m=0;m<3;m++){
+	comI[m] = 0.0;
+	comJ[m] = 0.0;
+      }
+      for(j=0;j<nLipAtoms;j++){
+	l=k+j;
+      
+	for(m=0;m<3;m++){
+	  comI[m] += atomsI[l].mass * atomsI[l].pos[m];
+	  comJ[m] += atomsJ[l].mass * atomsJ[l].pos[m];
+	  
+	  totMass += atomsI[l].mass;
+	}
+	for(m=0;m<3;m++){
+	  comI[m] /= totMass;
+	  comJ[m] /= totMass;
+	}
+      }
+      
+      for(j=0;j<3;j++)
+	d[j] = comJ[j] - comI[j];
+      
+      dSqr = 0.0;
+      for(j=0;j<3;j++)
+	dSqr += d[j] * d[j];
+      
+      accum += sqrt(dSqr);
+      nAccums++;
+    }
+  }
+
+  accum /= (double)nAccums;
+  rmsdCorr[bin].rmsd += accum;
+  rmsdCorr[bin].nValues++;
+}
+  
+  
+