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 (21 years, 8 months ago) by mmeineke
Content type:	text/plain
File size:	7095 byte(s)
Log Message:	added gofz code
#	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
33	// 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
45	FILE* outFile;
46
47	int i, j, k;
48	int startFrame, corrFrames, framesFinished;
49	int startFound, percentComplete;
50	int index;
51	int nCounts;
52	int bin;
53
54	double Hmat[3][3];
55	double diffTime;
56	double rmsdDt;
57	double timeOut, outVal;
58
59	struct atomCoord* atoms;
60
61	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
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
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	myFrames[index].time = frameTimes[i];
173	readFrame(i, atoms, myFrames[index].Hmat );
174	}
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