trunk/tcProps/directorHead.c

#define _FILE_OFFSET_BITS 64

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


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


struct directStr{
  double uTop[3];
  double uBottom[3];
  double orderTop;
  double orderBottom;
  double time;
};

struct directStr* directorHead;


void accumDHFrame( int index, struct atomCoord *atoms );

void calcDirHeadCorr(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;
  directorHead = (struct directStr*)calloc(corrFrames, 
                                           sizeof(struct directStr));
  

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

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

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

    framesFinished++;
  }

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


//   fflush(outFile);
//   fclose(outFile);

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


void accumDHFrame( int index, 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;
  const double oneThird = 1.0 / 3.0;

  int i,j,k;
  int nTop;
  int nBot;
  int lWork;
  int nfilled;
  int ndiag;
  int ifail;

  double oTop[3][3];   
  double oBottom[3][3];
  double evals[3];
  double work[9];
  double *u;
  double max;
  int which;

  char job, uplo;

  job = 'V';
  uplo = 'U';
  nfilled = 3;
  ndiag = 3;
  lWork = 9;

  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      oTop[i][j] = 0.0;   
      oBottom[i][j] = 0.0;
    }
  }

  nTop = 0;
  nBot = 0;
  for(i=0;i<nLipids;i++){
    
    k = i*nLipAtoms;
    u = atoms[k].u;
    if(atoms[k].pos[2] > 0){
      
      oTop[0][0] += u[0] * u[0] - oneThird;
      oTop[0][1] += u[0] * u[1];
      oTop[0][2] += u[0] * u[2];

      oTop[1][0] += u[1] * u[0];
      oTop[1][1] += u[1] * u[1] - oneThird;
      oTop[1][2] += u[1] * u[2];

      oTop[2][0] += u[2] * u[0];
      oTop[2][1] += u[2] * u[1];
      oTop[2][2] += u[2] * u[2] - oneThird;
    
      nTop++;
    }
    else{
      
      oBottom[0][0] += u[0] * u[0] - oneThird;
      oBottom[0][1] += u[0] * u[1];
      oBottom[0][2] += u[0] * u[2];

      oBottom[1][0] += u[1] * u[0];
      oBottom[1][1] += u[1] * u[1] - oneThird;
      oBottom[1][2] += u[1] * u[2];

      oBottom[2][0] += u[2] * u[0];
      oBottom[2][1] += u[2] * u[1];
      oBottom[2][2] += u[2] * u[2] - oneThird;
    
      nBot++;
    }
  }

  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      oTop[i][j] /= (double)nTop;   
      oBottom[i][j] /= (double)nBot;
    }
  }

  ifail = 0;

  dsyev(&job, &uplo, &nfilled, oTop, &ndiag, evals, work, &lWork, &ifail);

  if (ifail) {
    fprintf(stderr, "dsyev screwed something up!\n");
    exit(0);
  }
  
  max = 0.0;
  for (i=0; i<3;i++) {
    if (fabs(evals[i]) > max) {
      which = i;
      max = fabs(evals[i]);
    }
  }

  for (i = 0; i < 3; i++) {
    directorHead[index].uTop[i] = oTop[which][i];
  }

  directorHead[index].orderTop = 1.5 * max;

  ifail = 0;

  dsyev(&job, &uplo, &nfilled, oBottom, &ndiag, evals, work, &lWork, &ifail);

  if (ifail) {
    fprintf(stderr, "dsyev screwed something up!\n");
    exit(0);
  }
  
  max = 0.0;
  for (i=0; i<3;i++) {
    if (fabs(evals[i]) > max) {
      which = i;
      max = fabs(evals[i]);
    }
  }

  for (i = 0; i < 3; i++) {
    directorHead[index].uBottom[i] = oBottom[which][i];
  }

  directorHead[index].orderBottom = 1.5 * max;
  
  directorHead[index].time = frameTimes[index];

  fprintf(stderr,
          "frame[%d] => orderTop = %6G; < %6G, %6G, %6G >\n"
          "          orderBottom = %6G; < %6G, %6G, %6G >\n\n",
          index,
          directorHead[index].orderTop,
          directorHead[index].uTop[0],
          directorHead[index].uTop[1],
          directorHead[index].uTop[2],
          directorHead[index].orderBottom,
          directorHead[index].uBottom[0],
          directorHead[index].uBottom[1],
          directorHead[index].uBottom[2] );
}

Revision:	1060
Committed:	Thu Feb 19 21:10:06 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	text/plain
File size:	5182 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			#include <mkl_lapack64.h>
7
8
9			#include "params.h"
10			#include "tcProps.h"
11			#include "readWrite.h"
12			#include "directorHead.h"
13
14
15
16			struct directStr{
17			double uTop[3];
18			double uBottom[3];
19			double orderTop;
20			double orderBottom;
21			double time;
22			};
23
24			struct directStr* directorHead;
25
26
27			void accumDHFrame( int index, struct atomCoord *atoms );
28
29			void calcDirHeadCorr(double startTime, struct atomCoord* atoms,
30			char* outPrefix ){
31
32			// list of 'a priori' constants
33
34			const int nLipAtoms = NL_ATOMS;
35			const int nBonds = NBONDS;
36			const int nLipids = NLIPIDS;
37			const int nSSD = NSSD;
38			const int nAtoms = nLipAtoms * nLipids + nSSD;
39
40			// variables
41
42			char outName[500];
43			FILE* outFile;
44			int i, j, k;
45			int startFrame, corrFrames, framesFinished;
46			int startFound, percentComplete;
47
48			double Hmat[3][3];
49
50			framesFinished = 0;
51
52			startFound = 0;
53			startFrame = -1;
54			while( !startFound ){
55
56			startFrame++;
57
58			if(startTime <= frameTimes[startFrame])
59			startFound = 1;
60
61			if(startFrame >= nFrames){
62
63			fprintf( stderr,
64			"Start Time, %G, was not found in the dump file.\n",
65			startTime );
66			exit(0);
67			}
68			}
69
70			corrFrames = nFrames - startFrame;
71			directorHead = (struct directStr*)calloc(corrFrames,
72			sizeof(struct directStr));
73
74
75			for(i=startFrame; i<nFrames; i++){
76
77			percentComplete =
78			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
79
80			// fprintf( stdout,
81			// "\rDirector head corr %3d%% complete.",
82			// percentComplete );
83			// fflush( stdout );
84
85			readFrame( i, atoms, Hmat );
86
87			accumDHFrame( i, atoms );
88
89			framesFinished++;
90			}
91
92			// sprintf( outName, "%s.dirHead", outPrefix );
93			// outFile = fopen( outName, "w" );
94
95
96
97
98			// fflush(outFile);
99			// fclose(outFile);
100
101			percentComplete =
102			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
103
104			fprintf( stdout,
105			"\rDirector head corr %3d%% complete.",
106			"done.\n",
107			percentComplete );
108			fflush( stdout );
109			}
110
111
112			void accumDHFrame( int index, struct atomCoord *atoms ){
113
114			// list of 'a priori' constants
115
116			const int nLipAtoms = NL_ATOMS;
117			const int nBonds = NBONDS;
118			const int nLipids = NLIPIDS;
119			const int nSSD = NSSD;
120			const int nAtoms = nLipAtoms * nLipids + nSSD;
121			const double oneThird = 1.0 / 3.0;
122
123			int i,j,k;
124			int nTop;
125			int nBot;
126			int lWork;
127			int nfilled;
128			int ndiag;
129			int ifail;
130
131			double oTop[3][3];
132			double oBottom[3][3];
133			double evals[3];
134			double work[9];
135			double *u;
136			double max;
137			int which;
138
139			char job, uplo;
140
141			job = 'V';
142			uplo = 'U';
143			nfilled = 3;
144			ndiag = 3;
145			lWork = 9;
146
147			for(i=0;i<3;i++){
148			for(j=0;j<3;j++){
149			oTop[i][j] = 0.0;
150			oBottom[i][j] = 0.0;
151			}
152			}
153
154			nTop = 0;
155			nBot = 0;
156			for(i=0;i<nLipids;i++){
157
158			k = i*nLipAtoms;
159			u = atoms[k].u;
160			if(atoms[k].pos[2] > 0){
161
162			oTop[0][0] += u[0] * u[0] - oneThird;
163			oTop[0][1] += u[0] * u[1];
164			oTop[0][2] += u[0] * u[2];
165
166			oTop[1][0] += u[1] * u[0];
167			oTop[1][1] += u[1] * u[1] - oneThird;
168			oTop[1][2] += u[1] * u[2];
169
170			oTop[2][0] += u[2] * u[0];
171			oTop[2][1] += u[2] * u[1];
172			oTop[2][2] += u[2] * u[2] - oneThird;
173
174			nTop++;
175			}
176			else{
177
178			oBottom[0][0] += u[0] * u[0] - oneThird;
179			oBottom[0][1] += u[0] * u[1];
180			oBottom[0][2] += u[0] * u[2];
181
182			oBottom[1][0] += u[1] * u[0];
183			oBottom[1][1] += u[1] * u[1] - oneThird;
184			oBottom[1][2] += u[1] * u[2];
185
186			oBottom[2][0] += u[2] * u[0];
187			oBottom[2][1] += u[2] * u[1];
188			oBottom[2][2] += u[2] * u[2] - oneThird;
189
190			nBot++;
191			}
192			}
193
194			for(i=0;i<3;i++){
195			for(j=0;j<3;j++){
196			oTop[i][j] /= (double)nTop;
197			oBottom[i][j] /= (double)nBot;
198			}
199			}
200
201			ifail = 0;
202
203			dsyev(&job, &uplo, &nfilled, oTop, &ndiag, evals, work, &lWork, &ifail);
204
205			if (ifail) {
206			fprintf(stderr, "dsyev screwed something up!\n");
207			exit(0);
208			}
209
210			max = 0.0;
211			for (i=0; i<3;i++) {
212			if (fabs(evals[i]) > max) {
213			which = i;
214			max = fabs(evals[i]);
215			}
216			}
217
218			for (i = 0; i < 3; i++) {
219			directorHead[index].uTop[i] = oTop[which][i];
220			}
221
222			directorHead[index].orderTop = 1.5 * max;
223
224			ifail = 0;
225
226			dsyev(&job, &uplo, &nfilled, oBottom, &ndiag, evals, work, &lWork, &ifail);
227
228			if (ifail) {
229			fprintf(stderr, "dsyev screwed something up!\n");
230			exit(0);
231			}
232
233			max = 0.0;
234			for (i=0; i<3;i++) {
235			if (fabs(evals[i]) > max) {
236			which = i;
237			max = fabs(evals[i]);
238			}
239			}
240
241			for (i = 0; i < 3; i++) {
242			directorHead[index].uBottom[i] = oBottom[which][i];
243			}
244
245			directorHead[index].orderBottom = 1.5 * max;
246
247			directorHead[index].time = frameTimes[index];
248
249			fprintf(stderr,
250			"frame[%d] => orderTop = %6G; < %6G, %6G, %6G >\n"
251			" orderBottom = %6G; < %6G, %6G, %6G >\n\n",
252			index,
253			directorHead[index].orderTop,
254			directorHead[index].uTop[0],
255			directorHead[index].uTop[1],
256			directorHead[index].uTop[2],
257			directorHead[index].orderBottom,
258			directorHead[index].uBottom[0],
259			directorHead[index].uBottom[1],
260			directorHead[index].uBottom[2] );
261			}
262