trunk/tcProps/directorWhole.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 "directorWhole.h"

struct directStr{
  double u[3];
  double order;
  double time;
};

struct directStr* directorWhole;

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

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

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

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

    readFrame( i, atoms, Hmat );
    
    directorWhole[i-startFrame].time = frameTimes[i];

    accumDWFrame( i-startFrame, atoms );

    framesFinished++;
  }

  sprintf( outName, "%s.dirWhole", 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",
             directorWhole[i].time,
             directorWhole[i].order,
             directorWhole[i].u[0],
             directorWhole[i].u[1],
             directorWhole[i].u[2]);
  }

  
  fflush(outFile);
  fclose(outFile);

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


void accumDWFrame( 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,n;
  int lWork;
  int nfilled;
  int ndiag;
  int ifail;

  double oMat[3][3];   
  double iMat[3][3];
  double matWrap[9];
  double com[3];
  double totMass;
  double lenPrinc, distl;
  
  struct uVect{
    double u[3];
  };
  struct uVect principal[nLipids];


  double evals[3];
  double work[9];
  double *u;
  double smallest, 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++)
      oMat[i][j] = 0.0;
  
  for(i=0;i<nLipids;i++){
    
    k = i*nLipAtoms;
    
    // find the moment of inertia

    for(j=0;j<3;j++)
      for(l=0;l<3;l++)
        iMat[j][l] = 0.0;
    
    // com calc
    
    totMass = 0.0;
    for(m=0;m<3;m++) com[m] = 0.0;
    for(j=0;j<nLipAtoms;j++){
      l=k+j;
      
      for(m=0;m<3;m++)
        com[m] += atoms[l].mass * atoms[l].pos[m];
      
      totMass += atoms[l].mass;
    }
    for(m=0;m<3;m++) com[m] /= totMass;

    for(j=0;j<nLipAtoms;j++){
      l=k+j;
      

      distl = 0.0;
      for(m=0; m<3;m++) 
        distl += (atoms[l].pos[m]-com[m])*(atoms[l].pos[m]-com[m]);
      
      for(m=0; m<3; m++)
        iMat[m][m] += atoms[l].mass * distl;

      for(m=0;m<3;m++)
        for(n=0;n<3;n++)
          iMat[m][n] -= atoms[l].mass * ( atoms[l].pos[m] - com[m] ) *
            ( atoms[l].pos[n] - com[n] );
      
      
    }
    
//     fprintf(stderr, "iMat (before) = %6G, %6G, %6G\n"
//                  "                %6G, %6G, %6G\n"
//                  "                %6G, %6G, %6G\n\n",
//          iMat[0][0], iMat[0][1], iMat[0][2],
//          iMat[1][0], iMat[1][1], iMat[1][2],
//          iMat[2][0], iMat[2][1], iMat[2][2]);

    // matWrap to fortran convention

    for(j=0;j<3;j++)
      for(l=0;l<3;l++)
        matWrap[l+3*j] = iMat[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++)
        iMat[l][j] = matWrap[l+3*j];

//     fprintf(stderr, "evals = %6G, %6G, %6G\n\n",
//          evals[0], evals[1], evals[2]);

//     fprintf(stderr, "iMat (after) = %6G, %6G, %6G\n"
//                  "               %6G, %6G, %6G\n"
//                  "               %6G, %6G, %6G\n\n",
//          iMat[0][0], iMat[0][1], iMat[0][2],
//          iMat[1][0], iMat[1][1], iMat[1][2],
//          iMat[2][0], iMat[2][1], iMat[2][2]);


    smallest = fabs(evals[0]);
    which = 0;
    for(j=0;j<3;j++){
      if( fabs(evals[j]) < smallest ){
        which = j;
        smallest = fabs(evals[j]);
      }
    }

    lenPrinc = 0.0;
    for(j=0;j<3;j++)
      lenPrinc +=  iMat[j][which] * iMat[j][which];
    lenPrinc = sqrt( lenPrinc );
    
    for(j=0;j<3;j++)
      principal[i].u[j] = iMat[j][which]/lenPrinc;


    // find the director of the bilayer
    
    u = principal[i].u;
    
    oMat[0][0] += u[0] * u[0] - oneThird;
    oMat[0][1] += u[0] * u[1];
    oMat[0][2] += u[0] * u[2];
    
    oMat[1][0] += u[1] * u[0];
    oMat[1][1] += u[1] * u[1] - oneThird;
    oMat[1][2] += u[1] * u[2];
    
    oMat[2][0] += u[2] * u[0];
    oMat[2][1] += u[2] * u[1];
    oMat[2][2] += u[2] * u[2] - oneThird;
  }

  for(i=0;i<3;i++)
    for(j=0;j<3;j++)
      oMat[i][j] /= (double)nLipids;

  // matWrap to fortran convention
  
  for(j=0;j<3;j++)
    for(l=0;l<3;l++)
      matWrap[l+3*j] = oMat[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++)
      oMat[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++) {
    directorWhole[index].u[i] = oMat[i][which];
  }

  directorWhole[index].order = 1.5 * max;

//   fprintf(stderr,
//        "frame[%d] => order = %6G; < %6G, %6G, %6G >\n",
//        index,
//        directorWhole[index].order,
//        directorWhole[index].u[0],
//        directorWhole[index].u[1],
//        directorWhole[index].u[2] );

//   fprintf(stdout,
//        "%6g\t%6G\n",
//        directorWhole[index].time,
//        directorWhole[index].order);
}
Revision:	1067
Committed:	Tue Feb 24 17:13:06 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	text/plain
File size:	7021 byte(s)
Log Message:	worked out the diagonalization bugs in head and whole
#	User	Rev	Content
1	mmeineke	1062	#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 "directorWhole.h"
13
14			struct directStr{
15			double u[3];
16			double order;
17			double time;
18			};
19
20			struct directStr* directorWhole;
21
22			void accumDWFrame( int index, struct atomCoord *atoms );
23
24			void calcDirWholeCorr(double startTime, struct atomCoord* atoms,
25			char* outPrefix ){
26
27			// list of 'a priori' constants
28
29			const int nLipAtoms = NL_ATOMS;
30			const int nBonds = NBONDS;
31			const int nLipids = NLIPIDS;
32			const int nSSD = NSSD;
33			const int nAtoms = nLipAtoms * nLipids + nSSD;
34
35			// variables
36
37			char outName[500];
38			FILE* outFile;
39			int i, j, k;
40			int startFrame, corrFrames, framesFinished;
41			int startFound, percentComplete;
42
43			double Hmat[3][3];
44
45			framesFinished = 0;
46
47			startFound = 0;
48			startFrame = -1;
49			while( !startFound ){
50
51			startFrame++;
52
53			if(startFrame >= nFrames){
54
55			fprintf( stderr,
56			"Start Time, %G, was not found in the dump file.\n",
57			startTime );
58			exit(0);
59			}
60
61			if(startTime <= frameTimes[startFrame])
62			startFound = 1;
63
64
65			}
66
67			corrFrames = nFrames - startFrame;
68			directorWhole = (struct directStr*)calloc(corrFrames,
69			sizeof(struct directStr));
70
71
72			for(i=startFrame; i<nFrames; i++){
73
74			percentComplete =
75			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
76
77	mmeineke	1067	fprintf( stdout,
78			"\rDirector bilayer corr %3d%% complete.",
79			percentComplete );
80			fflush( stdout );
81	mmeineke	1062
82			readFrame( i, atoms, Hmat );
83
84			directorWhole[i-startFrame].time = frameTimes[i];
85
86			accumDWFrame( i-startFrame, atoms );
87
88			framesFinished++;
89			}
90
91	mmeineke	1067	sprintf( outName, "%s.dirWhole", outPrefix );
92			outFile = fopen( outName, "w" );
93	mmeineke	1062
94	mmeineke	1067	fprintf( outFile,
95			"#time\torderParam\tx\ty\tz\n");
96
97			for(i=0;i<corrFrames;i++){
98			fprintf( outFile,
99			"%6G\t%6G\t%6G\t%6G\t%6G\n",
100			directorWhole[i].time,
101			directorWhole[i].order,
102			directorWhole[i].u[0],
103			directorWhole[i].u[1],
104			directorWhole[i].u[2]);
105			}
106	mmeineke	1062
107
108
109	mmeineke	1067	fflush(outFile);
110			fclose(outFile);
111	mmeineke	1062
112			percentComplete =
113			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
114
115			fprintf( stdout,
116			"\rDirector bilayer corr %3d%% complete.\n"
117			"done.\n",
118			percentComplete );
119			fflush( stdout );
120			}
121
122
123			void accumDWFrame( int index, struct atomCoord *atoms ){
124
125			// list of 'a priori' constants
126
127			const int nLipAtoms = NL_ATOMS;
128			const int nBonds = NBONDS;
129			const int nLipids = NLIPIDS;
130			const int nSSD = NSSD;
131			const int nAtoms = nLipAtoms * nLipids + nSSD;
132			const double oneThird = 1.0 / 3.0;
133
134			int i,j,k,l,m,n;
135			int lWork;
136			int nfilled;
137			int ndiag;
138			int ifail;
139
140			double oMat[3][3];
141			double iMat[3][3];
142	mmeineke	1067	double matWrap[9];
143	mmeineke	1062	double com[3];
144			double totMass;
145			double lenPrinc, distl;
146
147			struct uVect{
148			double u[3];
149			};
150			struct uVect principal[nLipids];
151
152
153			double evals[3];
154			double work[9];
155			double *u;
156			double smallest, max;
157			int which;
158
159			char job, uplo;
160
161			job = 'V';
162			uplo = 'U';
163			nfilled = 3;
164			ndiag = 3;
165			lWork = 9;
166
167			for(i=0;i<3;i++)
168			for(j=0;j<3;j++)
169			oMat[i][j] = 0.0;
170
171			for(i=0;i<nLipids;i++){
172
173			k = i*nLipAtoms;
174
175			// find the moment of inertia
176
177			for(j=0;j<3;j++)
178			for(l=0;l<3;l++)
179			iMat[j][l] = 0.0;
180
181			// com calc
182
183			totMass = 0.0;
184			for(m=0;m<3;m++) com[m] = 0.0;
185			for(j=0;j<nLipAtoms;j++){
186			l=k+j;
187
188	mmeineke	1067	for(m=0;m<3;m++)
189	mmeineke	1062	com[m] += atoms[l].mass * atoms[l].pos[m];
190	mmeineke	1067
191			totMass += atoms[l].mass;
192	mmeineke	1062	}
193			for(m=0;m<3;m++) com[m] /= totMass;
194
195			for(j=0;j<nLipAtoms;j++){
196			l=k+j;
197
198
199			distl = 0.0;
200			for(m=0; m<3;m++)
201			distl += (atoms[l].pos[m]-com[m])*(atoms[l].pos[m]-com[m]);
202
203			for(m=0; m<3; m++)
204			iMat[m][m] += atoms[l].mass * distl;
205
206			for(m=0;m<3;m++)
207			for(n=0;n<3;n++)
208			iMat[m][n] -= atoms[l].mass * ( atoms[l].pos[m] - com[m] ) *
209			( atoms[l].pos[n] - com[n] );
210
211
212			}
213	mmeineke	1067
214			// fprintf(stderr, "iMat (before) = %6G, %6G, %6G\n"
215			// " %6G, %6G, %6G\n"
216			// " %6G, %6G, %6G\n\n",
217			// iMat[0][0], iMat[0][1], iMat[0][2],
218			// iMat[1][0], iMat[1][1], iMat[1][2],
219			// iMat[2][0], iMat[2][1], iMat[2][2]);
220	mmeineke	1062
221	mmeineke	1067	// matWrap to fortran convention
222
223			for(j=0;j<3;j++)
224			for(l=0;l<3;l++)
225			matWrap[l+3*j] = iMat[l][j];
226
227	mmeineke	1062	ifail = 0;
228	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
229	mmeineke	1062
230			if (ifail) {
231			fprintf(stderr, "dsyev screwed something up!\n");
232			exit(0);
233			}
234
235	mmeineke	1067	// matWrap from fortran convention
236	mmeineke	1062
237	mmeineke	1067	for(j=0;j<3;j++)
238			for(l=0;l<3;l++)
239			iMat[l][j] = matWrap[l+3*j];
240	mmeineke	1062
241	mmeineke	1067	// fprintf(stderr, "evals = %6G, %6G, %6G\n\n",
242			// evals[0], evals[1], evals[2]);
243	mmeineke	1062
244	mmeineke	1067	// fprintf(stderr, "iMat (after) = %6G, %6G, %6G\n"
245			// " %6G, %6G, %6G\n"
246			// " %6G, %6G, %6G\n\n",
247			// iMat[0][0], iMat[0][1], iMat[0][2],
248			// iMat[1][0], iMat[1][1], iMat[1][2],
249			// iMat[2][0], iMat[2][1], iMat[2][2]);
250
251
252	mmeineke	1062
253			smallest = fabs(evals[0]);
254			which = 0;
255			for(j=0;j<3;j++){
256			if( fabs(evals[j]) < smallest ){
257			which = j;
258			smallest = fabs(evals[j]);
259			}
260			}
261
262			lenPrinc = 0.0;
263			for(j=0;j<3;j++)
264			lenPrinc += iMat[j][which] * iMat[j][which];
265	mmeineke	1067	lenPrinc = sqrt( lenPrinc );
266	mmeineke	1062
267			for(j=0;j<3;j++)
268			principal[i].u[j] = iMat[j][which]/lenPrinc;
269
270
271
272
273
274			// find the director of the bilayer
275
276			u = principal[i].u;
277
278			oMat[0][0] += u[0] * u[0] - oneThird;
279			oMat[0][1] += u[0] * u[1];
280			oMat[0][2] += u[0] * u[2];
281
282			oMat[1][0] += u[1] * u[0];
283			oMat[1][1] += u[1] * u[1] - oneThird;
284			oMat[1][2] += u[1] * u[2];
285
286			oMat[2][0] += u[2] * u[0];
287			oMat[2][1] += u[2] * u[1];
288			oMat[2][2] += u[2] * u[2] - oneThird;
289			}
290
291			for(i=0;i<3;i++)
292			for(j=0;j<3;j++)
293			oMat[i][j] /= (double)nLipids;
294
295	mmeineke	1067	// matWrap to fortran convention
296
297			for(j=0;j<3;j++)
298			for(l=0;l<3;l++)
299			matWrap[l+3*j] = oMat[l][j];
300
301	mmeineke	1062	ifail = 0;
302	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
303	mmeineke	1062
304			if (ifail) {
305			fprintf(stderr, "dsyev screwed something up!\n");
306			exit(0);
307			}
308	mmeineke	1067
309			// matWrap from fortran convention
310	mmeineke	1062
311	mmeineke	1067	for(j=0;j<3;j++)
312			for(l=0;l<3;l++)
313			oMat[l][j] = matWrap[l+3*j];
314
315	mmeineke	1062	max = 0.0;
316			for (i=0; i<3;i++) {
317			if (fabs(evals[i]) > max) {
318			which = i;
319			max = fabs(evals[i]);
320			}
321			}
322
323			for (i = 0; i < 3; i++) {
324			directorWhole[index].u[i] = oMat[i][which];
325			}
326
327			directorWhole[index].order = 1.5 * max;
328
329	mmeineke	1067	// fprintf(stderr,
330			// "frame[%d] => order = %6G; < %6G, %6G, %6G >\n",
331			// index,
332			// directorWhole[index].order,
333			// directorWhole[index].u[0],
334			// directorWhole[index].u[1],
335			// directorWhole[index].u[2] );
336	mmeineke	1062
337			// fprintf(stdout,
338			// "%6g\t%6G\n",
339			// directorWhole[index].time,
340			// directorWhole[index].order);
341			}