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,
           char* theName){

  // 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 );
    
    index++;
  }

  // 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,
                 "\rMSD 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-%s.msd", outPrefix, theName );
  outFile = fopen( outName, "w" );

  fprintf( outFile,
           "#time\tmsd\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,
           "\rMSD corr %3d%% complete.\n",
           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;

  accum = 0.0;

  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 += 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 += dSqr;
      nAccums++;
    }
  }

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