trunk/tcProps/scdCorr.c

#define _FILE_OFFSET_BITS 64

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#include "params.h"
#include "tcProps.h"
#include "readWrite.h"
#include "scdCorr.h"

double scdParamTot[9];

void accumScdFrame( struct atomCoord* atoms );
void scdSet( int index, struct atomCoord* atoms, double *Sxx, 
             double *Syy, int nLipid, int a, int b, int c );

void calcScdCorr( double startTime, struct atomCoord* atoms, 
                  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;
  
  double Hmat[3][3];

  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;

  for(i=0;i<9;i++)
    scdParamTot[i] = 0.0;

  for(i=startFrame; i<nFrames; i++){

    percentComplete = 
      (int)( 100.0 * (double)framesFinished / (double) corrFrames );
    
    fprintf( stdout,
             "\rScd corr %3d%% complete.",
             percentComplete );
    fflush( stdout );

    readFrame( i, atoms, Hmat );
    
    accumScdFrame( atoms );

    framesFinished++;
  }

  for(i=0;i<9;i++)
    scdParamTot[i] /= (double)corrFrames;

  sprintf( outName, "%s.scdCorr", outPrefix );
  outFile = fopen( outName, "w" );

  fprintf( outFile,
           "#C_num\tScd\n" );
  for(i=1;i<8;i++)
    fprintf( outFile,
             "%d\t%G\n",
             i+1,
             scdParamTot[i] );
  fflush(outFile);
  fclose(outFile);

  percentComplete = 
    (int)( 100.0 * (double)framesFinished / (double) corrFrames );
  
  fprintf( stdout,
           "\rScd corr %3d%% complete.\n"
           "done.\n",
           percentComplete );
  fflush( stdout );
}

void accumScdFrame( struct atomCoord* atoms ){

  // 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;

  int index;
  int i,j,k;
  
  double scdParam[9];
  double Sxx[9];
  double Syy[9];

  for(i=0;i<9;i++){
    Sxx[i] = 0.0;
    Syy[i] = 0.0;
  }
    
  
  for(i=0;i<nLipids;i++){

    
    // Sxx average first

    // chain 1
    
    scdSet( 1, atoms, Sxx, Syy, i, 2, 3, 4 );
    scdSet( 2, atoms, Sxx, Syy, i, 3, 4, 5 );
    scdSet( 3, atoms, Sxx, Syy, i, 4, 5, 6 );
    scdSet( 4, atoms, Sxx, Syy, i, 5, 6, 7 );
    scdSet( 5, atoms, Sxx, Syy, i, 6, 7, 8 );
    scdSet( 6, atoms, Sxx, Syy, i, 7, 8, 9 );
    scdSet( 7, atoms, Sxx, Syy, i, 8, 9, 10 );
    
    // chain 2

    scdSet( 1, atoms, Sxx, Syy, i, 2, 11, 12 );
    scdSet( 2, atoms, Sxx, Syy, i, 11, 12, 13 );
    scdSet( 3, atoms, Sxx, Syy, i, 12, 13, 14 );
    scdSet( 4, atoms, Sxx, Syy, i, 13, 14, 15 );
    scdSet( 5, atoms, Sxx, Syy, i, 14, 15, 16 );
    scdSet( 6, atoms, Sxx, Syy, i, 15, 16, 17 );
    scdSet( 7, atoms, Sxx, Syy, i, 16, 17, 18 );
  }

  // take the average of the frame

  Sxx[0] /= 2*nLipids;
  for(i=1;i<9;i++)
    Sxx[i] /= 4*nLipids;

  Syy[0] /= 2*nLipids;
  for(i=1;i<9;i++)
    Syy[i] /= 4*nLipids;

  for(i=0;i<9;i++)
    scdParam[i] = 2.0*Sxx[i]/3.0 + Syy[i]/3.0;

  // accumulate into the running average
  
  for(i=0;i<9;i++)
    scdParamTot[i] += scdParam[i];
  
 
}

void scdSet( int index, struct atomCoord* atoms, double *Sxx, 
             double *Syy, int nLipid, int a, int b, int c ){

  int i,k;

  double dist,dSqr;
  double molX[3];
  double molY[3];
  double molZ[3];
  double temp[3];

  k = nLipid*NL_ATOMS;
  a += k;
  b += k;
  c += k;

  // calc the molecular z axis

  for(i=0;i<3;i++)
    molZ[i] = atoms[c].pos[i] - atoms[a].pos[i];

  dSqr = 0.0;
  for(i=0;i<3;i++)
    dSqr += molZ[i] * molZ[i];

  dist = sqrt(dSqr);
  for(i=0;i<3;i++)
    molZ[i] /= dist;

  // calc the temp vector

  for(i=0;i<3;i++)
    temp[i] = atoms[a].pos[i] - atoms[b].pos[i];

  dSqr = 0.0;
  for(i=0;i<3;i++)
    dSqr += temp[i] * temp[i];

  dist = sqrt(dSqr);
  for(i=0;i<3;i++)
    temp[i] /= dist;

  
  // calc the molecular x axis ( x = temp cross z )

  molX[0] = temp[1]*molZ[2] - molZ[1]*temp[2];
  molX[1] = temp[2]*molZ[0] - molZ[2]*temp[0];
  molX[2] = temp[0]*molZ[1] - molZ[0]*temp[1];

  // calc the y axis ( y = z cross x )
  
  molY[0] = molZ[1]*molX[2] - molX[1]*molZ[2];
  molY[1] = molZ[2]*molX[0] - molX[2]*molZ[0];
  molY[2] = molZ[0]*molX[1] - molX[0]*molZ[1];
  
  // calculate the Sxx and Syy params
  
  Sxx[index] += 3 * molX[2] * molX[2] - 1;
  Syy[index] += 3 * molY[2] * molY[2] - 1;
}
Revision:	1062
Committed:	Fri Feb 20 21:20:37 2004 UTC (21 years, 8 months ago) by mmeineke
Content type:	text/plain
File size:	4909 byte(s)
Log Message:	debugging the director code
#	User	Rev	Content
1	mmeineke	1060	#define _FILE_OFFSET_BITS 64
2
3			#include <stdio.h>
4			#include <stdlib.h>
5			#include <math.h>
6
7
8			#include "params.h"
9			#include "tcProps.h"
10			#include "readWrite.h"
11			#include "scdCorr.h"
12
13			double scdParamTot[9];
14
15			void accumScdFrame( struct atomCoord* atoms );
16			void scdSet( int index, struct atomCoord* atoms, double *Sxx,
17			double *Syy, int nLipid, int a, int b, int c );
18
19			void calcScdCorr( double startTime, struct atomCoord* atoms,
20			char* outPrefix ){
21
22			// list of 'a priori' constants
23
24			const int nLipAtoms = NL_ATOMS;
25			const int nBonds = NBONDS;
26			const int nLipids = NLIPIDS;
27			const int nSSD = NSSD;
28			const int nAtoms = nLipAtoms * nLipids + nSSD;
29
30			// variables
31
32			char outName[500];
33			FILE* outFile;
34			int i, j, k;
35			int startFrame, corrFrames, framesFinished;
36			int startFound, percentComplete;
37
38			double Hmat[3][3];
39
40			framesFinished = 0;
41
42			startFound = 0;
43			startFrame = -1;
44			while( !startFound ){
45
46			startFrame++;
47
48			if(startFrame >= nFrames){
49
50			fprintf( stderr,
51			"Start Time, %G, was not found in the dump file.\n",
52			startTime );
53			exit(0);
54			}
55	mmeineke	1062
56			if(startTime <= frameTimes[startFrame])
57			startFound = 1;
58
59
60	mmeineke	1060	}
61
62			corrFrames = nFrames - startFrame;
63
64			for(i=0;i<9;i++)
65			scdParamTot[i] = 0.0;
66
67			for(i=startFrame; i<nFrames; i++){
68
69			percentComplete =
70			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
71
72			fprintf( stdout,
73			"\rScd corr %3d%% complete.",
74			percentComplete );
75			fflush( stdout );
76
77			readFrame( i, atoms, Hmat );
78
79			accumScdFrame( atoms );
80
81			framesFinished++;
82			}
83
84			for(i=0;i<9;i++)
85			scdParamTot[i] /= (double)corrFrames;
86
87			sprintf( outName, "%s.scdCorr", outPrefix );
88			outFile = fopen( outName, "w" );
89
90			fprintf( outFile,
91			"#C_num\tScd\n" );
92			for(i=1;i<8;i++)
93			fprintf( outFile,
94			"%d\t%G\n",
95			i+1,
96			scdParamTot[i] );
97			fflush(outFile);
98			fclose(outFile);
99
100			percentComplete =
101			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
102
103			fprintf( stdout,
104			"\rScd corr %3d%% complete.\n"
105			"done.\n",
106			percentComplete );
107			fflush( stdout );
108			}
109
110			void accumScdFrame( struct atomCoord* atoms ){
111
112			// list of 'a priori' constants
113
114			const int nLipAtoms = NL_ATOMS;
115			const int nBonds = NBONDS;
116			const int nLipids = NLIPIDS;
117			const int nSSD = NSSD;
118			const int nAtoms = nLipAtoms * nLipids + nSSD;
119
120			int index;
121			int i,j,k;
122
123			double scdParam[9];
124			double Sxx[9];
125			double Syy[9];
126
127			for(i=0;i<9;i++){
128			Sxx[i] = 0.0;
129			Syy[i] = 0.0;
130			}
131
132
133			for(i=0;i<nLipids;i++){
134
135
136			// Sxx average first
137
138			// chain 1
139
140			scdSet( 1, atoms, Sxx, Syy, i, 2, 3, 4 );
141			scdSet( 2, atoms, Sxx, Syy, i, 3, 4, 5 );
142			scdSet( 3, atoms, Sxx, Syy, i, 4, 5, 6 );
143			scdSet( 4, atoms, Sxx, Syy, i, 5, 6, 7 );
144			scdSet( 5, atoms, Sxx, Syy, i, 6, 7, 8 );
145			scdSet( 6, atoms, Sxx, Syy, i, 7, 8, 9 );
146			scdSet( 7, atoms, Sxx, Syy, i, 8, 9, 10 );
147
148			// chain 2
149
150			scdSet( 1, atoms, Sxx, Syy, i, 2, 11, 12 );
151			scdSet( 2, atoms, Sxx, Syy, i, 11, 12, 13 );
152			scdSet( 3, atoms, Sxx, Syy, i, 12, 13, 14 );
153			scdSet( 4, atoms, Sxx, Syy, i, 13, 14, 15 );
154			scdSet( 5, atoms, Sxx, Syy, i, 14, 15, 16 );
155			scdSet( 6, atoms, Sxx, Syy, i, 15, 16, 17 );
156			scdSet( 7, atoms, Sxx, Syy, i, 16, 17, 18 );
157			}
158
159			// take the average of the frame
160
161			Sxx[0] /= 2*nLipids;
162			for(i=1;i<9;i++)
163			Sxx[i] /= 4*nLipids;
164
165			Syy[0] /= 2*nLipids;
166			for(i=1;i<9;i++)
167			Syy[i] /= 4*nLipids;
168
169			for(i=0;i<9;i++)
170			scdParam[i] = 2.0*Sxx[i]/3.0 + Syy[i]/3.0;
171
172			// accumulate into the running average
173
174			for(i=0;i<9;i++)
175			scdParamTot[i] += scdParam[i];
176
177
178			}
179
180			void scdSet( int index, struct atomCoord* atoms, double *Sxx,
181			double *Syy, int nLipid, int a, int b, int c ){
182
183			int i,k;
184
185			double dist,dSqr;
186			double molX[3];
187			double molY[3];
188			double molZ[3];
189			double temp[3];
190
191			k = nLipid*NL_ATOMS;
192			a += k;
193			b += k;
194			c += k;
195
196			// calc the molecular z axis
197
198			for(i=0;i<3;i++)
199			molZ[i] = atoms[c].pos[i] - atoms[a].pos[i];
200
201			dSqr = 0.0;
202			for(i=0;i<3;i++)
203			dSqr += molZ[i] * molZ[i];
204
205			dist = sqrt(dSqr);
206			for(i=0;i<3;i++)
207			molZ[i] /= dist;
208
209			// calc the temp vector
210
211			for(i=0;i<3;i++)
212			temp[i] = atoms[a].pos[i] - atoms[b].pos[i];
213
214			dSqr = 0.0;
215			for(i=0;i<3;i++)
216			dSqr += temp[i] * temp[i];
217
218			dist = sqrt(dSqr);
219			for(i=0;i<3;i++)
220			temp[i] /= dist;
221
222
223			// calc the molecular x axis ( x = temp cross z )
224
225			molX[0] = temp[1]molZ[2] - molZ[1]temp[2];
226			molX[1] = temp[2]molZ[0] - molZ[2]temp[0];
227			molX[2] = temp[0]molZ[1] - molZ[0]temp[1];
228
229			// calc the y axis ( y = z cross x )
230
231			molY[0] = molZ[1]molX[2] - molX[1]molZ[2];
232			molY[1] = molZ[2]molX[0] - molX[2]molZ[0];
233			molY[2] = molZ[0]molX[1] - molX[0]molZ[1];
234
235			// calculate the Sxx and Syy params
236
237			Sxx[index] += 3 * molX[2] * molX[2] - 1;
238			Syy[index] += 3 * molY[2] * molY[2] - 1;
239			}