[ViewVC] Diff of: group/trunk/tcProps/rmsd.c

Comparing trunk/tcProps/rmsd.c (file contents):
Revision 1072 by mmeineke, Fri Feb 27 20:39:17 2004 UTC vs.
Revision 1080 by mmeineke, Wed Mar 3 15:16:15 2004 UTC

+  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;
+  // 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;
-+
+  struct atomCoord* atoms;
-+
+  framesFinished = 0;
+  startFound = 0;
+  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 );
-+
+    }
-+
-+
+    // initialize the arrays
-+
-+
+    atoms = myFrames[index].atoms;
-+
+    for(j=0;j<nLipids;j++){
-+
-+
+      atoms[nLipAtoms*j+0].type = HEAD;
-+
+      atoms[nLipAtoms*j+0].mass = 72;
-+
+      atoms[nLipAtoms*j+0].u[0] = 0.0;
-+
+      atoms[nLipAtoms*j+0].u[1] = 0.0;
-+
+      atoms[nLipAtoms*j+0].u[2] = 1.0;
-+
-+
-+
-+
+      atoms[nLipAtoms*j+1].type = CH2;
-+
+      atoms[nLipAtoms*j+1].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+2].type = CH;
-+
+      atoms[nLipAtoms*j+2].mass = 13.02;
-+
-+
+      atoms[nLipAtoms*j+3].type = CH2;
-+
+      atoms[nLipAtoms*j+3].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+4].type = CH2;
-+
+      atoms[nLipAtoms*j+4].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+5].type = CH2;
-+
+      atoms[nLipAtoms*j+5].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+6].type = CH2;
-+
+      atoms[nLipAtoms*j+6].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+7].type = CH2;
-+
+      atoms[nLipAtoms*j+7].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+8].type = CH2;
-+
+      atoms[nLipAtoms*j+8].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+9].type = CH2;
-+
+      atoms[nLipAtoms*j+9].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+10].type = CH3;
-+
+      atoms[nLipAtoms*j+10].mass = 15.04;
-+
-+
+      atoms[nLipAtoms*j+11].type = CH2;
-+
+      atoms[nLipAtoms*j+11].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+12].type = CH2;
-+
+      atoms[nLipAtoms*j+12].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+13].type = CH2;
-+
+      atoms[nLipAtoms*j+13].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+14].type = CH2;
-+
+      atoms[nLipAtoms*j+14].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+15].type = CH2;
-+
+      atoms[nLipAtoms*j+15].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+16].type = CH2;
-+
+      atoms[nLipAtoms*j+16].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+17].type = CH2;
-+
+      atoms[nLipAtoms*j+17].mass = 14.03;
-+
-+
+      atoms[nLipAtoms*j+18].type = CH3;
-+
+      atoms[nLipAtoms*j+18].mass = 15.04;
-+
+    }
-+
-+
+    for(j=(nLipAtoms*nLipids);j<nAtoms;j++){
-+
+      atoms[j].type = SSD;
-+
+      atoms[j].mass = 18.03;
-+
+      atoms[j].u[0] = 0.0;
-+
+      atoms[j].u[1] = 0.0;
-+
+      atoms[j].u[2] = 1.0;
-+
+    }
-+
+    myFrames[index].time = frameTimes[i];
-+
+    readFrame(i, 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++;
-+
+}
-+
-+
-+

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing trunk/tcProps/rmsd.c (file contents): Revision 1072 by mmeineke, Fri Feb 27 20:39:17 2004 UTC vs. Revision 1080 by mmeineke, Wed Mar 3 15:16:15 2004 UTC

Diff Legend

Comparing trunk/tcProps/rmsd.c (file contents):
Revision 1072 by mmeineke, Fri Feb 27 20:39:17 2004 UTC vs.
Revision 1080 by mmeineke, Wed Mar 3 15:16:15 2004 UTC