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, 3 months ago) by mmeineke
Content type:	text/plain
File size:	7159 byte(s)
Log Message:	who knows?, but it works for the most part
#	Content
1	#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	struct binStruct{
21	double rmsd;
22	int nValues;
23	};
24
25	struct binStruct rmsdCorr[RMSDBINS];
26
27	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	char* theName){
33
34	// 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
46	FILE* outFile;
47
48	int i, j, k;
49	int startFrame, corrFrames, framesFinished;
50	int startFound, percentComplete;
51	int index;
52	int nCounts;
53	int bin;
54
55	double Hmat[3][3];
56	double diffTime;
57	double rmsdDt;
58	double timeOut, outVal;
59
60	struct atomCoord* atoms;
61
62	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
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
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	myFrames[index].time = frameTimes[i];
174	readFrame(i, atoms, myFrames[index].Hmat );
175
176	index++;
177	}
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
201	if( diffTime > 0.0 ){
202
203	percentComplete =
204	(int)( 100.0 * (double)index / (double) nCounts );
205
206	fprintf( stdout,
207	"\rMSD corr %3d%% complete.",
208	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	sprintf( outName, "%s-%s.msd", outPrefix, theName );
223	outFile = fopen( outName, "w" );
224
225	fprintf( outFile,
226	"#time\tmsd\n");
227
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	"\rMSD corr %3d%% complete.\n",
249	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	accum = 0.0;
287
288	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	accum += dSqr;
303	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	accum += dSqr;
344	nAccums++;
345	}
346	}
347
348	accum /= (double)nAccums;
349	rmsdCorr[bin].rmsd += accum;
350	rmsdCorr[bin].nValues++;
351	}
352
353
354