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.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.\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;
  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:	1062
Committed:	Fri Feb 20 21:20:37 2004 UTC (21 years, 8 months ago) by mmeineke
Content type:	text/plain
File size:	5192 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			#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			// fprintf( stdout,
82			// "\rDirector head corr %3d%% complete.",
83			// percentComplete );
84			// fflush( stdout );
85
86			readFrame( i, atoms, Hmat );
87
88	mmeineke	1062	accumDHFrame( i-startFrame, atoms );
89	mmeineke	1060
90			framesFinished++;
91			}
92
93			// sprintf( outName, "%s.dirHead", outPrefix );
94			// outFile = fopen( outName, "w" );
95
96
97
98
99			// fflush(outFile);
100			// fclose(outFile);
101
102			percentComplete =
103			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
104
105			fprintf( stdout,
106	mmeineke	1062	"\rDirector head corr %3d%% complete.\n"
107	mmeineke	1060	"done.\n",
108			percentComplete );
109			fflush( stdout );
110			}
111
112
113			void accumDHFrame( int index, struct atomCoord *atoms ){
114
115			// list of 'a priori' constants
116
117			const int nLipAtoms = NL_ATOMS;
118			const int nBonds = NBONDS;
119			const int nLipids = NLIPIDS;
120			const int nSSD = NSSD;
121			const int nAtoms = nLipAtoms * nLipids + nSSD;
122			const double oneThird = 1.0 / 3.0;
123
124			int i,j,k;
125			int nTop;
126			int nBot;
127			int lWork;
128			int nfilled;
129			int ndiag;
130			int ifail;
131
132			double oTop[3][3];
133			double oBottom[3][3];
134			double evals[3];
135			double work[9];
136			double *u;
137			double max;
138			int which;
139
140			char job, uplo;
141
142			job = 'V';
143			uplo = 'U';
144			nfilled = 3;
145			ndiag = 3;
146			lWork = 9;
147
148			for(i=0;i<3;i++){
149			for(j=0;j<3;j++){
150			oTop[i][j] = 0.0;
151			oBottom[i][j] = 0.0;
152			}
153			}
154
155			nTop = 0;
156			nBot = 0;
157			for(i=0;i<nLipids;i++){
158
159			k = i*nLipAtoms;
160			u = atoms[k].u;
161			if(atoms[k].pos[2] > 0){
162
163			oTop[0][0] += u[0] * u[0] - oneThird;
164			oTop[0][1] += u[0] * u[1];
165			oTop[0][2] += u[0] * u[2];
166
167			oTop[1][0] += u[1] * u[0];
168			oTop[1][1] += u[1] * u[1] - oneThird;
169			oTop[1][2] += u[1] * u[2];
170
171			oTop[2][0] += u[2] * u[0];
172			oTop[2][1] += u[2] * u[1];
173			oTop[2][2] += u[2] * u[2] - oneThird;
174
175			nTop++;
176			}
177			else{
178
179			oBottom[0][0] += u[0] * u[0] - oneThird;
180			oBottom[0][1] += u[0] * u[1];
181			oBottom[0][2] += u[0] * u[2];
182
183			oBottom[1][0] += u[1] * u[0];
184			oBottom[1][1] += u[1] * u[1] - oneThird;
185			oBottom[1][2] += u[1] * u[2];
186
187			oBottom[2][0] += u[2] * u[0];
188			oBottom[2][1] += u[2] * u[1];
189			oBottom[2][2] += u[2] * u[2] - oneThird;
190
191			nBot++;
192			}
193			}
194
195			for(i=0;i<3;i++){
196			for(j=0;j<3;j++){
197			oTop[i][j] /= (double)nTop;
198			oBottom[i][j] /= (double)nBot;
199			}
200			}
201
202			ifail = 0;
203
204			dsyev(&job, &uplo, &nfilled, oTop, &ndiag, evals, work, &lWork, &ifail);
205
206			if (ifail) {
207			fprintf(stderr, "dsyev screwed something up!\n");
208			exit(0);
209			}
210
211			max = 0.0;
212			for (i=0; i<3;i++) {
213			if (fabs(evals[i]) > max) {
214			which = i;
215			max = fabs(evals[i]);
216			}
217			}
218
219			for (i = 0; i < 3; i++) {
220			directorHead[index].uTop[i] = oTop[which][i];
221			}
222
223			directorHead[index].orderTop = 1.5 * max;
224
225			ifail = 0;
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