trunk/tcProps/rmsd.c

#define _FILE_OFFSET_BITS 64

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


#include "params.h"
#include "tcProps.h"
#include "readWrite.h"
#include "rmsd.h"

struct frameStruct{
  struct atomCoord atoms[NL_ATOMS*NLIPIDS+NSSD];
  double time;
  double Hmat[3][3];
};

struct binStruct{
  double rmsd;
  int nValues;
};

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;
  const int nBonds    = NBONDS;
  const int nLipids   = NLIPIDS;
  const int nSSD      = NSSD;
  const int nAtoms    = nLipAtoms * nLipids + nSSD;

  // 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;
  
  startFound = 0;
  startFrame = -1;
  while( !startFound ){

    startFrame++;

    if(startFrame >= nFrames){
      
      fprintf( stderr,
               "Start Time, %G, was not found in the dump file.\n",
               startTime );
      exit(0);
    }
    
    if(startTime <= frameTimes[startFrame])
      startFound = 1;


  }

  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++;
}
  
  
Revision:	1073
Committed:	Sat Feb 28 16:45:57 2004 UTC (20 years, 6 months ago) by mmeineke
Content type:	text/plain
File size:	4982 byte(s)
Log Message:	finished the rmsd coding. hasn't been debugged
#	User	Rev	Content
1	mmeineke	1072	#define _FILE_OFFSET_BITS 64
2
3			#include <stdio.h>
4			#include <stdlib.h>
5			#include <string.h>
6			#include <math.h>
7
8
9			#include "params.h"
10			#include "tcProps.h"
11			#include "readWrite.h"
12			#include "rmsd.h"
13
14			struct frameStruct{
15			struct atomCoord atoms[NL_ATOMS*NLIPIDS+NSSD];
16			double time;
17			double Hmat[3][3];
18			};
19
20	mmeineke	1073	struct binStruct{
21			double rmsd;
22			int nValues;
23			};
24	mmeineke	1072
25	mmeineke	1073	struct binStruct rmsdCorr[RMSDBINS];
26	mmeineke	1072
27	mmeineke	1073	struct frameStruct* myFrames;
28
29			void calcRMSDpair( int bin, int i, int j, enum atomNames rType );
30
31			void rmsd( enum atomNames rType, double startTime, char* outPrefix ){
32
33	mmeineke	1072	// list of 'a priori' constants
34
35			const int nLipAtoms = NL_ATOMS;
36			const int nBonds = NBONDS;
37			const int nLipids = NLIPIDS;
38			const int nSSD = NSSD;
39			const int nAtoms = nLipAtoms * nLipids + nSSD;
40
41			// variables
42
43			char outName[500];
44	mmeineke	1073
45	mmeineke	1072	FILE* outFile;
46	mmeineke	1073
47	mmeineke	1072	int i, j, k;
48			int startFrame, corrFrames, framesFinished;
49			int startFound, percentComplete;
50	mmeineke	1073	int index;
51			int nCounts;
52			int bin;
53
54	mmeineke	1072	double Hmat[3][3];
55	mmeineke	1073	double diffTime;
56			double rmsdDt;
57			double timeOut, outVal;
58	mmeineke	1072
59			framesFinished = 0;
60
61			startFound = 0;
62			startFrame = -1;
63			while( !startFound ){
64
65			startFrame++;
66
67			if(startFrame >= nFrames){
68
69			fprintf( stderr,
70			"Start Time, %G, was not found in the dump file.\n",
71			startTime );
72			exit(0);
73			}
74
75			if(startTime <= frameTimes[startFrame])
76			startFound = 1;
77
78
79			}
80
81			corrFrames = nFrames - startFrame;
82			myFrames = (struct frameStruct*)calloc(corrFrames,
83			sizeof(struct frameStruct));
84	mmeineke	1073
85			index=0;
86			for(i=startFrame;i<nFrames;i++){
87
88			if( index >= corrFrames ){
89			fprintf( stderr,
90			"Error, number of frames, %d, exceeds corrFrames, %d.\n",
91			index,
92			corrFrames );
93			}
94
95			myFrames[index].time = frameTimes[i];
96			readFrame(i, myFrames[index].atoms, myFrames[index].Hmat );
97			}
98
99			// initialize the counts and the correlation
100
101			nCounts = 0;
102			for(i=0;i<(corrFrames-1);i++)
103			for(j=i+1;j<corrFrames;j++)
104			nCounts++;
105
106			for(i=0;i<RMSDBINS;i++){
107			rmsdCorr[i].rmsd = 0.0;
108			rmsdCorr[i].nValues = 0;
109			}
110
111			rmsdDt = (myFrames[corrFrames-1].time - myFrames[0].time) / RMSDBINS;
112
113			// perform the calculation
114
115			index = 0;
116			for(i=0;i<(corrFrames-1);i++){
117			for(j=i+1;j<corrFrames;j++){
118
119			diffTime = myFrames[j].time - myFrames[i].time;
120			if( diffTime > 0.0 ){
121
122			percentComplete =
123			(int)( 100.0 * (double)index / (double) nCounts );
124
125			fprintf( stdout,
126			"\rRMSD corr %3d%% complete.",
127			percentComplete );
128			fflush( stdout );
129
130			bin = (int)(diffTime / rmsdDt);
131			if( bin < RMSDBINS)
132			calcRMSDpair( bin, i, j, rType );
133			}
134
135			index++;
136			}
137			}
138
139			// print out the correlation
140
141			sprintf( outName, "%s.rmsd", outPrefix );
142			outFile = fopen( outName, "w" );
143
144			fprintf( outFile,
145			"#time\trmsd\n");
146
147			for(i=0;i<RMSDBINS;i++){
148
149			if(rmsdCorr[i].nValues){
150
151			timeOut = (2.0(double)i+1.0)rmsdDt*0.5;
152			outVal = rmsdCorr[i].rmsd / (double)rmsdCorr[i].nValues;
153
154			fprintf( outFile,
155			"%6G\t%6G\n",
156			timeOut, outVal );
157			}
158			}
159
160			fflush(outFile);
161			fclose(outFile);
162
163			percentComplete =
164			(int)( 100.0 * (double)index / (double) nCounts );
165
166			fprintf( stdout,
167			"\rRMSD corr %3d%% complete.",
168			percentComplete );
169			fflush( stdout );
170
171			free( myFrames );
172			myFrames = NULL;
173			}
174
175
176			void calcRMSDpair( int bin, int frameI, int frameJ, enum atomNames rType ){
177
178			struct atomCoord* atomsI;
179			struct atomCoord* atomsJ;
180
181			int i,j,k,l,m,n;
182			int nAccums;
183
184			double d[3];
185			double dSqr;
186			double comI[3];
187			double comJ[3];
188			double totMass;
189			double accum;
190
191
192
193			// list of 'a priori' constants
194
195			const int nLipAtoms = NL_ATOMS;
196			const int nBonds = NBONDS;
197			const int nLipids = NLIPIDS;
198			const int nSSD = NSSD;
199			const int nAtoms = nLipAtoms * nLipids + nSSD;
200
201
202			atomsI = myFrames[frameI].atoms;
203			atomsJ = myFrames[frameJ].atoms;
204
205			if( rType != COM ){
206
207			nAccums = 0;
208			for(i=0;i<nAtoms;i++){
209
210			if (atomsI[i].type == rType){
211
212			for(j=0;j<3;j++)
213			d[j] = atomsJ[i].pos[j] - atomsI[i].pos[j];
214
215			dSqr = 0.0;
216			for(j=0;j<3;j++)
217			dSqr += d[j] * d[j];
218
219			accum += sqrt(dSqr);
220			nAccums++;
221			}
222			}
223			}
224			else{
225
226			nAccums = 0;
227			for(i=0;i<nLipids;i++){
228
229			k = i*nLipAtoms;
230
231			// com calc
232
233			totMass = 0.0;
234			for(m=0;m<3;m++){
235			comI[m] = 0.0;
236			comJ[m] = 0.0;
237			}
238			for(j=0;j<nLipAtoms;j++){
239			l=k+j;
240
241			for(m=0;m<3;m++){
242			comI[m] += atomsI[l].mass * atomsI[l].pos[m];
243			comJ[m] += atomsJ[l].mass * atomsJ[l].pos[m];
244
245			totMass += atomsI[l].mass;
246			}
247			for(m=0;m<3;m++){
248			comI[m] /= totMass;
249			comJ[m] /= totMass;
250			}
251			}
252
253			for(j=0;j<3;j++)
254			d[j] = comJ[j] - comI[j];
255
256			dSqr = 0.0;
257			for(j=0;j<3;j++)
258			dSqr += d[j] * d[j];
259
260			accum += sqrt(dSqr);
261			nAccums++;
262			}
263			}
264
265			accum /= (double)nAccums;
266			rmsdCorr[bin].rmsd += accum;
267			rmsdCorr[bin].nValues++;
268			}
269
270
271