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;

  struct atomCoord* atoms;

  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 );
    }
     
    // 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++;
}
  
  
Revision:	1080
Committed:	Wed Mar 3 15:16:15 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	text/plain
File size:	7095 byte(s)
Log Message:	added gofz code
#	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	mmeineke	1080	struct atomCoord* atoms;
60
61	mmeineke	1072	framesFinished = 0;
62
63			startFound = 0;
64			startFrame = -1;
65			while( !startFound ){
66
67			startFrame++;
68
69			if(startFrame >= nFrames){
70
71			fprintf( stderr,
72			"Start Time, %G, was not found in the dump file.\n",
73			startTime );
74			exit(0);
75			}
76
77			if(startTime <= frameTimes[startFrame])
78			startFound = 1;
79
80
81			}
82
83			corrFrames = nFrames - startFrame;
84			myFrames = (struct frameStruct*)calloc(corrFrames,
85			sizeof(struct frameStruct));
86	mmeineke	1073
87			index=0;
88			for(i=startFrame;i<nFrames;i++){
89
90			if( index >= corrFrames ){
91			fprintf( stderr,
92			"Error, number of frames, %d, exceeds corrFrames, %d.\n",
93			index,
94			corrFrames );
95			}
96	mmeineke	1080
97			// initialize the arrays
98
99			atoms = myFrames[index].atoms;
100			for(j=0;j<nLipids;j++){
101
102			atoms[nLipAtoms*j+0].type = HEAD;
103			atoms[nLipAtoms*j+0].mass = 72;
104			atoms[nLipAtoms*j+0].u[0] = 0.0;
105			atoms[nLipAtoms*j+0].u[1] = 0.0;
106			atoms[nLipAtoms*j+0].u[2] = 1.0;
107
108
109
110			atoms[nLipAtoms*j+1].type = CH2;
111			atoms[nLipAtoms*j+1].mass = 14.03;
112
113			atoms[nLipAtoms*j+2].type = CH;
114			atoms[nLipAtoms*j+2].mass = 13.02;
115
116			atoms[nLipAtoms*j+3].type = CH2;
117			atoms[nLipAtoms*j+3].mass = 14.03;
118
119			atoms[nLipAtoms*j+4].type = CH2;
120			atoms[nLipAtoms*j+4].mass = 14.03;
121
122			atoms[nLipAtoms*j+5].type = CH2;
123			atoms[nLipAtoms*j+5].mass = 14.03;
124
125			atoms[nLipAtoms*j+6].type = CH2;
126			atoms[nLipAtoms*j+6].mass = 14.03;
127
128			atoms[nLipAtoms*j+7].type = CH2;
129			atoms[nLipAtoms*j+7].mass = 14.03;
130
131			atoms[nLipAtoms*j+8].type = CH2;
132			atoms[nLipAtoms*j+8].mass = 14.03;
133
134			atoms[nLipAtoms*j+9].type = CH2;
135			atoms[nLipAtoms*j+9].mass = 14.03;
136
137			atoms[nLipAtoms*j+10].type = CH3;
138			atoms[nLipAtoms*j+10].mass = 15.04;
139
140			atoms[nLipAtoms*j+11].type = CH2;
141			atoms[nLipAtoms*j+11].mass = 14.03;
142
143			atoms[nLipAtoms*j+12].type = CH2;
144			atoms[nLipAtoms*j+12].mass = 14.03;
145
146			atoms[nLipAtoms*j+13].type = CH2;
147			atoms[nLipAtoms*j+13].mass = 14.03;
148
149			atoms[nLipAtoms*j+14].type = CH2;
150			atoms[nLipAtoms*j+14].mass = 14.03;
151
152			atoms[nLipAtoms*j+15].type = CH2;
153			atoms[nLipAtoms*j+15].mass = 14.03;
154
155			atoms[nLipAtoms*j+16].type = CH2;
156			atoms[nLipAtoms*j+16].mass = 14.03;
157
158			atoms[nLipAtoms*j+17].type = CH2;
159			atoms[nLipAtoms*j+17].mass = 14.03;
160
161			atoms[nLipAtoms*j+18].type = CH3;
162			atoms[nLipAtoms*j+18].mass = 15.04;
163			}
164
165			for(j=(nLipAtoms*nLipids);j<nAtoms;j++){
166			atoms[j].type = SSD;
167			atoms[j].mass = 18.03;
168			atoms[j].u[0] = 0.0;
169			atoms[j].u[1] = 0.0;
170			atoms[j].u[2] = 1.0;
171			}
172	mmeineke	1073	myFrames[index].time = frameTimes[i];
173	mmeineke	1080	readFrame(i, atoms, myFrames[index].Hmat );
174	mmeineke	1073	}
175
176			// initialize the counts and the correlation
177
178			nCounts = 0;
179			for(i=0;i<(corrFrames-1);i++)
180			for(j=i+1;j<corrFrames;j++)
181			nCounts++;
182
183			for(i=0;i<RMSDBINS;i++){
184			rmsdCorr[i].rmsd = 0.0;
185			rmsdCorr[i].nValues = 0;
186			}
187
188			rmsdDt = (myFrames[corrFrames-1].time - myFrames[0].time) / RMSDBINS;
189
190			// perform the calculation
191
192			index = 0;
193			for(i=0;i<(corrFrames-1);i++){
194			for(j=i+1;j<corrFrames;j++){
195
196			diffTime = myFrames[j].time - myFrames[i].time;
197			if( diffTime > 0.0 ){
198
199			percentComplete =
200			(int)( 100.0 * (double)index / (double) nCounts );
201
202			fprintf( stdout,
203			"\rRMSD corr %3d%% complete.",
204			percentComplete );
205			fflush( stdout );
206
207			bin = (int)(diffTime / rmsdDt);
208			if( bin < RMSDBINS)
209			calcRMSDpair( bin, i, j, rType );
210			}
211
212			index++;
213			}
214			}
215
216			// print out the correlation
217
218			sprintf( outName, "%s.rmsd", outPrefix );
219			outFile = fopen( outName, "w" );
220
221			fprintf( outFile,
222			"#time\trmsd\n");
223
224			for(i=0;i<RMSDBINS;i++){
225
226			if(rmsdCorr[i].nValues){
227
228			timeOut = (2.0(double)i+1.0)rmsdDt*0.5;
229			outVal = rmsdCorr[i].rmsd / (double)rmsdCorr[i].nValues;
230
231			fprintf( outFile,
232			"%6G\t%6G\n",
233			timeOut, outVal );
234			}
235			}
236
237			fflush(outFile);
238			fclose(outFile);
239
240			percentComplete =
241			(int)( 100.0 * (double)index / (double) nCounts );
242
243			fprintf( stdout,
244			"\rRMSD corr %3d%% complete.",
245			percentComplete );
246			fflush( stdout );
247
248			free( myFrames );
249			myFrames = NULL;
250			}
251
252
253			void calcRMSDpair( int bin, int frameI, int frameJ, enum atomNames rType ){
254
255			struct atomCoord* atomsI;
256			struct atomCoord* atomsJ;
257
258			int i,j,k,l,m,n;
259			int nAccums;
260
261			double d[3];
262			double dSqr;
263			double comI[3];
264			double comJ[3];
265			double totMass;
266			double accum;
267
268
269
270			// list of 'a priori' constants
271
272			const int nLipAtoms = NL_ATOMS;
273			const int nBonds = NBONDS;
274			const int nLipids = NLIPIDS;
275			const int nSSD = NSSD;
276			const int nAtoms = nLipAtoms * nLipids + nSSD;
277
278
279			atomsI = myFrames[frameI].atoms;
280			atomsJ = myFrames[frameJ].atoms;
281
282			if( rType != COM ){
283
284			nAccums = 0;
285			for(i=0;i<nAtoms;i++){
286
287			if (atomsI[i].type == rType){
288
289			for(j=0;j<3;j++)
290			d[j] = atomsJ[i].pos[j] - atomsI[i].pos[j];
291
292			dSqr = 0.0;
293			for(j=0;j<3;j++)
294			dSqr += d[j] * d[j];
295
296			accum += sqrt(dSqr);
297			nAccums++;
298			}
299			}
300			}
301			else{
302
303			nAccums = 0;
304			for(i=0;i<nLipids;i++){
305
306			k = i*nLipAtoms;
307
308			// com calc
309
310			totMass = 0.0;
311			for(m=0;m<3;m++){
312			comI[m] = 0.0;
313			comJ[m] = 0.0;
314			}
315			for(j=0;j<nLipAtoms;j++){
316			l=k+j;
317
318			for(m=0;m<3;m++){
319			comI[m] += atomsI[l].mass * atomsI[l].pos[m];
320			comJ[m] += atomsJ[l].mass * atomsJ[l].pos[m];
321
322			totMass += atomsI[l].mass;
323			}
324			for(m=0;m<3;m++){
325			comI[m] /= totMass;
326			comJ[m] /= totMass;
327			}
328			}
329
330			for(j=0;j<3;j++)
331			d[j] = comJ[j] - comI[j];
332
333			dSqr = 0.0;
334			for(j=0;j<3;j++)
335			dSqr += d[j] * d[j];
336
337			accum += sqrt(dSqr);
338			nAccums++;
339			}
340			}
341
342			accum /= (double)nAccums;
343			rmsdCorr[bin].rmsd += accum;
344			rmsdCorr[bin].nValues++;
345			}
346
347
348