ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/tcProps/rmsd.c
(Generate patch)

Comparing trunk/tcProps/rmsd.c (file contents):
Revision 1072 by mmeineke, Fri Feb 27 20:39:17 2004 UTC vs.
Revision 1073 by mmeineke, Sat Feb 28 16:45:57 2004 UTC

# Line 17 | Line 17 | strcut frameStruct* myFrames;
17    double Hmat[3][3];
18   };
19  
20 < strcut frameStruct* myFrames;
20 > struct binStruct{
21 >  double rmsd;
22 >  int nValues;
23 > };
24  
25 < void rmsd(double startTime, char* outPrefix ){
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;
# Line 32 | Line 41 | void rmsd(double startTime, char* outPrefix ){
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 <  
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    framesFinished = 0;
60    
# Line 64 | Line 81 | void rmsd(double startTime, char* outPrefix ){
81    corrFrames = nFrames - startFrame;
82    myFrames = (struct frameStruct*)calloc(corrFrames,
83                                           sizeof(struct frameStruct));
84 +
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 +      

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines