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(startTime <= frameTimes[startFrame])
      startFound = 1;

    if(startFrame >= nFrames){
      
      fprintf( stderr,
               "Start Time, %G, was not found in the dump file.\n",
               startTime );
      exit(0);
    }
  }

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