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(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;
  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-startFrame, atoms );

    framesFinished++;
  }

  sprintf( outName, "%s.dirHeadTop", outPrefix );
  outFile = fopen( outName, "w" );
  fprintf( outFile,
           "#time\torderParam\tx\ty\tz\n");

  for(i=0;i<corrFrames;i++){
    fprintf( outFile,
             "%6G\t%6G\t%6G\t%6G\t%6G\n",
             directorHead[i].time,
             directorHead[i].orderTop,
             directorHead[i].uTop[0],
             directorHead[i].uTop[1],
             directorHead[i].uTop[2]);
  }
  fflush(outFile);
  fclose(outFile);

  sprintf( outName, "%s.dirHeadBottom", outPrefix );
  outFile = fopen( outName, "w" );
  fprintf( outFile,
           "#time\torderParam\tx\ty\tz\n");

  for(i=0;i<corrFrames;i++){
    fprintf( outFile,
             "%6G\t%6G\t%6G\t%6G\t%6G\n",
             directorHead[i].time,
             directorHead[i].orderBottom,
             directorHead[i].uBottom[0],
             directorHead[i].uBottom[1],
             directorHead[i].uBottom[2]);
  }
  fflush(outFile);
  fclose(outFile);


  percentComplete = 
    (int)( 100.0 * (double)framesFinished / (double) corrFrames );
  
  fprintf( stdout,
           "\rDirector head corr %3d%% complete.\n"
           "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,l,m;
  int nTop;
  int nBot;
  int lWork;
  int nfilled;
  int ndiag;
  int ifail;

  double oTop[3][3];   
  double oBottom[3][3];
  double matWrap[9];
  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;
    }
  }

  // matWrap to fortran convention
  
  for(j=0;j<3;j++)
    for(l=0;l<3;l++)
      matWrap[l+3*j] = oTop[l][j];

  ifail = 0;

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

  if (ifail) {
    fprintf(stderr, "dsyev screwed something up!\n");
    exit(0);
  }

  // matWrap from fortran convention
  
  for(j=0;j<3;j++)
    for(l=0;l<3;l++)
      oTop[l][j] = matWrap[l+3*j];
  
  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[i][which];
  }

  directorHead[index].orderTop = 1.5 * max;

  // matWrap to fortran convention
  
  for(j=0;j<3;j++)
    for(l=0;l<3;l++)
      matWrap[l+3*j] = oBottom[l][j];

  ifail = 0;
  dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);

  if (ifail) {
    fprintf(stderr, "dsyev screwed something up!\n");
    exit(0);
  }
  
  // matWrap from fortran convention
  
  for(j=0;j<3;j++)
    for(l=0;l<3;l++)
      oBottom[l][j] = matWrap[l+3*j];

  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[i][which];
  }

  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:	1067
Committed:	Tue Feb 24 17:13:06 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	text/plain
File size:	6439 byte(s)
Log Message:	worked out the diagonalization bugs in head and whole
#	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(startFrame >= nFrames){
59
60			fprintf( stderr,
61			"Start Time, %G, was not found in the dump file.\n",
62			startTime );
63			exit(0);
64			}
65	mmeineke	1062
66			if(startTime <= frameTimes[startFrame])
67			startFound = 1;
68
69
70	mmeineke	1060	}
71
72			corrFrames = nFrames - startFrame;
73			directorHead = (struct directStr*)calloc(corrFrames,
74			sizeof(struct directStr));
75
76			for(i=startFrame; i<nFrames; i++){
77
78			percentComplete =
79			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
80
81	mmeineke	1067	fprintf( stdout,
82			"\rDirector head corr %3d%% complete.",
83			percentComplete );
84			fflush( stdout );
85	mmeineke	1060
86			readFrame( i, atoms, Hmat );
87
88	mmeineke	1062	accumDHFrame( i-startFrame, atoms );
89	mmeineke	1060
90			framesFinished++;
91			}
92
93	mmeineke	1067	sprintf( outName, "%s.dirHeadTop", outPrefix );
94			outFile = fopen( outName, "w" );
95			fprintf( outFile,
96			"#time\torderParam\tx\ty\tz\n");
97	mmeineke	1060
98	mmeineke	1067	for(i=0;i<corrFrames;i++){
99			fprintf( outFile,
100			"%6G\t%6G\t%6G\t%6G\t%6G\n",
101			directorHead[i].time,
102			directorHead[i].orderTop,
103			directorHead[i].uTop[0],
104			directorHead[i].uTop[1],
105			directorHead[i].uTop[2]);
106			}
107			fflush(outFile);
108			fclose(outFile);
109	mmeineke	1060
110	mmeineke	1067	sprintf( outName, "%s.dirHeadBottom", outPrefix );
111			outFile = fopen( outName, "w" );
112			fprintf( outFile,
113			"#time\torderParam\tx\ty\tz\n");
114	mmeineke	1060
115	mmeineke	1067	for(i=0;i<corrFrames;i++){
116			fprintf( outFile,
117			"%6G\t%6G\t%6G\t%6G\t%6G\n",
118			directorHead[i].time,
119			directorHead[i].orderBottom,
120			directorHead[i].uBottom[0],
121			directorHead[i].uBottom[1],
122			directorHead[i].uBottom[2]);
123			}
124			fflush(outFile);
125			fclose(outFile);
126	mmeineke	1060
127	mmeineke	1067
128	mmeineke	1060	percentComplete =
129			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
130
131			fprintf( stdout,
132	mmeineke	1062	"\rDirector head corr %3d%% complete.\n"
133	mmeineke	1060	"done.\n",
134			percentComplete );
135			fflush( stdout );
136			}
137
138
139			void accumDHFrame( int index, struct atomCoord *atoms ){
140
141			// list of 'a priori' constants
142
143			const int nLipAtoms = NL_ATOMS;
144			const int nBonds = NBONDS;
145			const int nLipids = NLIPIDS;
146			const int nSSD = NSSD;
147			const int nAtoms = nLipAtoms * nLipids + nSSD;
148			const double oneThird = 1.0 / 3.0;
149
150	mmeineke	1067	int i,j,k,l,m;
151	mmeineke	1060	int nTop;
152			int nBot;
153			int lWork;
154			int nfilled;
155			int ndiag;
156			int ifail;
157
158			double oTop[3][3];
159			double oBottom[3][3];
160	mmeineke	1067	double matWrap[9];
161	mmeineke	1060	double evals[3];
162			double work[9];
163			double *u;
164			double max;
165			int which;
166
167			char job, uplo;
168
169			job = 'V';
170			uplo = 'U';
171			nfilled = 3;
172			ndiag = 3;
173			lWork = 9;
174
175			for(i=0;i<3;i++){
176			for(j=0;j<3;j++){
177			oTop[i][j] = 0.0;
178			oBottom[i][j] = 0.0;
179			}
180			}
181
182			nTop = 0;
183			nBot = 0;
184			for(i=0;i<nLipids;i++){
185
186			k = i*nLipAtoms;
187			u = atoms[k].u;
188			if(atoms[k].pos[2] > 0){
189
190			oTop[0][0] += u[0] * u[0] - oneThird;
191			oTop[0][1] += u[0] * u[1];
192			oTop[0][2] += u[0] * u[2];
193
194			oTop[1][0] += u[1] * u[0];
195			oTop[1][1] += u[1] * u[1] - oneThird;
196			oTop[1][2] += u[1] * u[2];
197
198			oTop[2][0] += u[2] * u[0];
199			oTop[2][1] += u[2] * u[1];
200			oTop[2][2] += u[2] * u[2] - oneThird;
201
202			nTop++;
203			}
204			else{
205
206			oBottom[0][0] += u[0] * u[0] - oneThird;
207			oBottom[0][1] += u[0] * u[1];
208			oBottom[0][2] += u[0] * u[2];
209
210			oBottom[1][0] += u[1] * u[0];
211			oBottom[1][1] += u[1] * u[1] - oneThird;
212			oBottom[1][2] += u[1] * u[2];
213
214			oBottom[2][0] += u[2] * u[0];
215			oBottom[2][1] += u[2] * u[1];
216			oBottom[2][2] += u[2] * u[2] - oneThird;
217
218			nBot++;
219			}
220			}
221
222			for(i=0;i<3;i++){
223			for(j=0;j<3;j++){
224			oTop[i][j] /= (double)nTop;
225			oBottom[i][j] /= (double)nBot;
226			}
227			}
228
229	mmeineke	1067	// matWrap to fortran convention
230
231			for(j=0;j<3;j++)
232			for(l=0;l<3;l++)
233			matWrap[l+3*j] = oTop[l][j];
234
235	mmeineke	1060	ifail = 0;
236
237	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
238	mmeineke	1060
239			if (ifail) {
240			fprintf(stderr, "dsyev screwed something up!\n");
241			exit(0);
242			}
243	mmeineke	1067
244			// matWrap from fortran convention
245	mmeineke	1060
246	mmeineke	1067	for(j=0;j<3;j++)
247			for(l=0;l<3;l++)
248			oTop[l][j] = matWrap[l+3*j];
249
250	mmeineke	1060	max = 0.0;
251			for (i=0; i<3;i++) {
252			if (fabs(evals[i]) > max) {
253			which = i;
254			max = fabs(evals[i]);
255			}
256			}
257
258			for (i = 0; i < 3; i++) {
259	mmeineke	1067	directorHead[index].uTop[i] = oTop[i][which];
260	mmeineke	1060	}
261
262			directorHead[index].orderTop = 1.5 * max;
263
264	mmeineke	1067	// matWrap to fortran convention
265
266			for(j=0;j<3;j++)
267			for(l=0;l<3;l++)
268			matWrap[l+3*j] = oBottom[l][j];
269
270	mmeineke	1060	ifail = 0;
271	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
272	mmeineke	1060
273			if (ifail) {
274			fprintf(stderr, "dsyev screwed something up!\n");
275			exit(0);
276			}
277
278	mmeineke	1067	// matWrap from fortran convention
279
280			for(j=0;j<3;j++)
281			for(l=0;l<3;l++)
282			oBottom[l][j] = matWrap[l+3*j];
283
284	mmeineke	1060	max = 0.0;
285			for (i=0; i<3;i++) {
286			if (fabs(evals[i]) > max) {
287			which = i;
288			max = fabs(evals[i]);
289			}
290			}
291
292			for (i = 0; i < 3; i++) {
293	mmeineke	1067	directorHead[index].uBottom[i] = oBottom[i][which];
294	mmeineke	1060	}
295
296			directorHead[index].orderBottom = 1.5 * max;
297
298			directorHead[index].time = frameTimes[index];
299
300	mmeineke	1067	// fprintf(stderr,
301			// "frame[%d] => orderTop = %6G; < %6G, %6G, %6G >\n"
302			// " orderBottom = %6G; < %6G, %6G, %6G >\n\n",
303			// index,
304			// directorHead[index].orderTop,
305			// directorHead[index].uTop[0],
306			// directorHead[index].uTop[1],
307			// directorHead[index].uTop[2],
308			// directorHead[index].orderBottom,
309			// directorHead[index].uBottom[0],
310			// directorHead[index].uBottom[1],
311			// directorHead[index].uBottom[2] );
312	mmeineke	1060	}
313