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

  free( directorHead );
  directorHead = NULL;
}


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:	1073
Committed:	Sat Feb 28 16:45:57 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	text/plain
File size:	6487 byte(s)
Log Message:	finished the rmsd coding. hasn't been debugged
#	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	mmeineke	1073
137			free( directorHead );
138			directorHead = NULL;
139	mmeineke	1060	}
140
141
142			void accumDHFrame( int index, struct atomCoord *atoms ){
143
144			// list of 'a priori' constants
145
146			const int nLipAtoms = NL_ATOMS;
147			const int nBonds = NBONDS;
148			const int nLipids = NLIPIDS;
149			const int nSSD = NSSD;
150			const int nAtoms = nLipAtoms * nLipids + nSSD;
151			const double oneThird = 1.0 / 3.0;
152
153	mmeineke	1067	int i,j,k,l,m;
154	mmeineke	1060	int nTop;
155			int nBot;
156			int lWork;
157			int nfilled;
158			int ndiag;
159			int ifail;
160
161			double oTop[3][3];
162			double oBottom[3][3];
163	mmeineke	1067	double matWrap[9];
164	mmeineke	1060	double evals[3];
165			double work[9];
166			double *u;
167			double max;
168			int which;
169
170			char job, uplo;
171
172			job = 'V';
173			uplo = 'U';
174			nfilled = 3;
175			ndiag = 3;
176			lWork = 9;
177
178			for(i=0;i<3;i++){
179			for(j=0;j<3;j++){
180			oTop[i][j] = 0.0;
181			oBottom[i][j] = 0.0;
182			}
183			}
184
185			nTop = 0;
186			nBot = 0;
187			for(i=0;i<nLipids;i++){
188
189			k = i*nLipAtoms;
190			u = atoms[k].u;
191			if(atoms[k].pos[2] > 0){
192
193			oTop[0][0] += u[0] * u[0] - oneThird;
194			oTop[0][1] += u[0] * u[1];
195			oTop[0][2] += u[0] * u[2];
196
197			oTop[1][0] += u[1] * u[0];
198			oTop[1][1] += u[1] * u[1] - oneThird;
199			oTop[1][2] += u[1] * u[2];
200
201			oTop[2][0] += u[2] * u[0];
202			oTop[2][1] += u[2] * u[1];
203			oTop[2][2] += u[2] * u[2] - oneThird;
204
205			nTop++;
206			}
207			else{
208
209			oBottom[0][0] += u[0] * u[0] - oneThird;
210			oBottom[0][1] += u[0] * u[1];
211			oBottom[0][2] += u[0] * u[2];
212
213			oBottom[1][0] += u[1] * u[0];
214			oBottom[1][1] += u[1] * u[1] - oneThird;
215			oBottom[1][2] += u[1] * u[2];
216
217			oBottom[2][0] += u[2] * u[0];
218			oBottom[2][1] += u[2] * u[1];
219			oBottom[2][2] += u[2] * u[2] - oneThird;
220
221			nBot++;
222			}
223			}
224
225			for(i=0;i<3;i++){
226			for(j=0;j<3;j++){
227			oTop[i][j] /= (double)nTop;
228			oBottom[i][j] /= (double)nBot;
229			}
230			}
231
232	mmeineke	1067	// matWrap to fortran convention
233
234			for(j=0;j<3;j++)
235			for(l=0;l<3;l++)
236			matWrap[l+3*j] = oTop[l][j];
237
238	mmeineke	1060	ifail = 0;
239
240	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
241	mmeineke	1060
242			if (ifail) {
243			fprintf(stderr, "dsyev screwed something up!\n");
244			exit(0);
245			}
246	mmeineke	1067
247			// matWrap from fortran convention
248	mmeineke	1060
249	mmeineke	1067	for(j=0;j<3;j++)
250			for(l=0;l<3;l++)
251			oTop[l][j] = matWrap[l+3*j];
252
253	mmeineke	1060	max = 0.0;
254			for (i=0; i<3;i++) {
255			if (fabs(evals[i]) > max) {
256			which = i;
257			max = fabs(evals[i]);
258			}
259			}
260
261			for (i = 0; i < 3; i++) {
262	mmeineke	1067	directorHead[index].uTop[i] = oTop[i][which];
263	mmeineke	1060	}
264
265			directorHead[index].orderTop = 1.5 * max;
266
267	mmeineke	1067	// matWrap to fortran convention
268
269			for(j=0;j<3;j++)
270			for(l=0;l<3;l++)
271			matWrap[l+3*j] = oBottom[l][j];
272
273	mmeineke	1060	ifail = 0;
274	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
275	mmeineke	1060
276			if (ifail) {
277			fprintf(stderr, "dsyev screwed something up!\n");
278			exit(0);
279			}
280
281	mmeineke	1067	// matWrap from fortran convention
282
283			for(j=0;j<3;j++)
284			for(l=0;l<3;l++)
285			oBottom[l][j] = matWrap[l+3*j];
286
287	mmeineke	1060	max = 0.0;
288			for (i=0; i<3;i++) {
289			if (fabs(evals[i]) > max) {
290			which = i;
291			max = fabs(evals[i]);
292			}
293			}
294
295			for (i = 0; i < 3; i++) {
296	mmeineke	1067	directorHead[index].uBottom[i] = oBottom[i][which];
297	mmeineke	1060	}
298
299			directorHead[index].orderBottom = 1.5 * max;
300
301			directorHead[index].time = frameTimes[index];
302
303	mmeineke	1067	// fprintf(stderr,
304			// "frame[%d] => orderTop = %6G; < %6G, %6G, %6G >\n"
305			// " orderBottom = %6G; < %6G, %6G, %6G >\n\n",
306			// index,
307			// directorHead[index].orderTop,
308			// directorHead[index].uTop[0],
309			// directorHead[index].uTop[1],
310			// directorHead[index].uTop[2],
311			// directorHead[index].orderBottom,
312			// directorHead[index].uBottom[0],
313			// directorHead[index].uBottom[1],
314			// directorHead[index].uBottom[2] );
315	mmeineke	1060	}
316