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.dirHead", outPrefix );
//   outFile = fopen( outName, "w" );


//   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 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] );
      
      
    }

    ifail = 0;
    dsyev(&job, &uplo, &nfilled, iMat, &ndiag, evals, work, &lWork, &ifail);
    
    if (ifail) {
      fprintf(stderr, "dsyev screwed something up!\n");
      exit(0);
    }

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

    fprintf(stderr, "iMat = %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];
    
    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;

  ifail = 0;

  dsyev(&job, &uplo, &nfilled, oMat, &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++) {
    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:	1062
Committed:	Fri Feb 20 21:20:37 2004 UTC (20 years, 6 months ago) by mmeineke
Content type:	text/plain
File size:	5858 byte(s)
Log Message:	debugging the director code
#	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			// fprintf( stdout,
78			// "\rDirector bilayer corr %3d%% complete.",
79			// percentComplete );
80			// fflush( stdout );
81
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			// sprintf( outName, "%s.dirHead", outPrefix );
92			// outFile = fopen( outName, "w" );
93
94
95
96
97			// fflush(outFile);
98			// fclose(outFile);
99
100			percentComplete =
101			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
102
103			fprintf( stdout,
104			"\rDirector bilayer corr %3d%% complete.\n"
105			"done.\n",
106			percentComplete );
107			fflush( stdout );
108			}
109
110
111			void accumDWFrame( int index, struct atomCoord *atoms ){
112
113			// list of 'a priori' constants
114
115			const int nLipAtoms = NL_ATOMS;
116			const int nBonds = NBONDS;
117			const int nLipids = NLIPIDS;
118			const int nSSD = NSSD;
119			const int nAtoms = nLipAtoms * nLipids + nSSD;
120			const double oneThird = 1.0 / 3.0;
121
122			int i,j,k,l,m,n;
123			int lWork;
124			int nfilled;
125			int ndiag;
126			int ifail;
127
128			double oMat[3][3];
129			double iMat[3][3];
130			double com[3];
131			double totMass;
132			double lenPrinc, distl;
133
134			struct uVect{
135			double u[3];
136			};
137			struct uVect principal[nLipids];
138
139
140			double evals[3];
141			double work[9];
142			double *u;
143			double smallest, max;
144			int which;
145
146			char job, uplo;
147
148			job = 'V';
149			uplo = 'U';
150			nfilled = 3;
151			ndiag = 3;
152			lWork = 9;
153
154			for(i=0;i<3;i++)
155			for(j=0;j<3;j++)
156			oMat[i][j] = 0.0;
157
158			for(i=0;i<nLipids;i++){
159
160			k = i*nLipAtoms;
161
162			// find the moment of inertia
163
164			for(j=0;j<3;j++)
165			for(l=0;l<3;l++)
166			iMat[j][l] = 0.0;
167
168			// com calc
169
170			totMass = 0.0;
171			for(m=0;m<3;m++) com[m] = 0.0;
172			for(j=0;j<nLipAtoms;j++){
173			l=k+j;
174
175			for(m=0;m<3;m++){
176			com[m] += atoms[l].mass * atoms[l].pos[m];
177			totMass += atoms[l].mass;
178			}
179			}
180			for(m=0;m<3;m++) com[m] /= totMass;
181
182			for(j=0;j<nLipAtoms;j++){
183			l=k+j;
184
185
186			distl = 0.0;
187			for(m=0; m<3;m++)
188			distl += (atoms[l].pos[m]-com[m])*(atoms[l].pos[m]-com[m]);
189
190			for(m=0; m<3; m++)
191			iMat[m][m] += atoms[l].mass * distl;
192
193			for(m=0;m<3;m++)
194			for(n=0;n<3;n++)
195			iMat[m][n] -= atoms[l].mass * ( atoms[l].pos[m] - com[m] ) *
196			( atoms[l].pos[n] - com[n] );
197
198
199			}
200
201			ifail = 0;
202			dsyev(&job, &uplo, &nfilled, iMat, &ndiag, evals, work, &lWork, &ifail);
203
204			if (ifail) {
205			fprintf(stderr, "dsyev screwed something up!\n");
206			exit(0);
207			}
208
209			fprintf(stderr, "evals = %6G, %6G, %6G\n\n",
210			evals[0], evals[1], evals[2]);
211
212			fprintf(stderr, "iMat = %6G, %6G, %6G\n"
213			" %6G, %6G, %6G\n"
214			" %6G, %6G, %6G\n\n",
215			iMat[0][0], iMat[0][1], iMat[0][2],
216			iMat[1][0], iMat[1][1], iMat[1][2],
217			iMat[2][0], iMat[2][1], iMat[2][2]);
218
219
220
221			smallest = fabs(evals[0]);
222			which = 0;
223			for(j=0;j<3;j++){
224			if( fabs(evals[j]) < smallest ){
225			which = j;
226			smallest = fabs(evals[j]);
227			}
228			}
229
230			lenPrinc = 0.0;
231			for(j=0;j<3;j++)
232			lenPrinc += iMat[j][which] * iMat[j][which];
233
234			for(j=0;j<3;j++)
235			principal[i].u[j] = iMat[j][which]/lenPrinc;
236
237
238
239
240
241			// find the director of the bilayer
242
243			u = principal[i].u;
244
245			oMat[0][0] += u[0] * u[0] - oneThird;
246			oMat[0][1] += u[0] * u[1];
247			oMat[0][2] += u[0] * u[2];
248
249			oMat[1][0] += u[1] * u[0];
250			oMat[1][1] += u[1] * u[1] - oneThird;
251			oMat[1][2] += u[1] * u[2];
252
253			oMat[2][0] += u[2] * u[0];
254			oMat[2][1] += u[2] * u[1];
255			oMat[2][2] += u[2] * u[2] - oneThird;
256			}
257
258			for(i=0;i<3;i++)
259			for(j=0;j<3;j++)
260			oMat[i][j] /= (double)nLipids;
261
262			ifail = 0;
263
264			dsyev(&job, &uplo, &nfilled, oMat, &ndiag, evals, work, &lWork, &ifail);
265
266			if (ifail) {
267			fprintf(stderr, "dsyev screwed something up!\n");
268			exit(0);
269			}
270
271			max = 0.0;
272			for (i=0; i<3;i++) {
273			if (fabs(evals[i]) > max) {
274			which = i;
275			max = fabs(evals[i]);
276			}
277			}
278
279			for (i = 0; i < 3; i++) {
280			directorWhole[index].u[i] = oMat[i][which];
281			}
282
283			directorWhole[index].order = 1.5 * max;
284
285			fprintf(stderr,
286			"frame[%d] => order = %6G; < %6G, %6G, %6G >\n",
287			index,
288			directorWhole[index].order,
289			directorWhole[index].u[0],
290			directorWhole[index].u[1],
291			directorWhole[index].u[2] );
292
293			// fprintf(stdout,
294			// "%6g\t%6G\n",
295			// directorWhole[index].time,
296			// directorWhole[index].order);
297			}