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

  free(directorWhole);
  directorWhole = NULL;
}


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:	1073
Committed:	Sat Feb 28 16:45:57 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	text/plain
File size:	7069 byte(s)
Log Message:	finished the rmsd coding. hasn't been debugged
#	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	mmeineke	1073
121			free(directorWhole);
122			directorWhole = NULL;
123	mmeineke	1062	}
124
125
126			void accumDWFrame( int index, struct atomCoord *atoms ){
127
128			// list of 'a priori' constants
129
130			const int nLipAtoms = NL_ATOMS;
131			const int nBonds = NBONDS;
132			const int nLipids = NLIPIDS;
133			const int nSSD = NSSD;
134			const int nAtoms = nLipAtoms * nLipids + nSSD;
135			const double oneThird = 1.0 / 3.0;
136
137			int i,j,k,l,m,n;
138			int lWork;
139			int nfilled;
140			int ndiag;
141			int ifail;
142
143			double oMat[3][3];
144			double iMat[3][3];
145	mmeineke	1067	double matWrap[9];
146	mmeineke	1062	double com[3];
147			double totMass;
148			double lenPrinc, distl;
149
150			struct uVect{
151			double u[3];
152			};
153			struct uVect principal[nLipids];
154
155
156			double evals[3];
157			double work[9];
158			double *u;
159			double smallest, max;
160			int which;
161
162			char job, uplo;
163
164			job = 'V';
165			uplo = 'U';
166			nfilled = 3;
167			ndiag = 3;
168			lWork = 9;
169
170			for(i=0;i<3;i++)
171			for(j=0;j<3;j++)
172			oMat[i][j] = 0.0;
173
174			for(i=0;i<nLipids;i++){
175
176			k = i*nLipAtoms;
177
178			// find the moment of inertia
179
180			for(j=0;j<3;j++)
181			for(l=0;l<3;l++)
182			iMat[j][l] = 0.0;
183
184			// com calc
185
186			totMass = 0.0;
187			for(m=0;m<3;m++) com[m] = 0.0;
188			for(j=0;j<nLipAtoms;j++){
189			l=k+j;
190
191	mmeineke	1067	for(m=0;m<3;m++)
192	mmeineke	1062	com[m] += atoms[l].mass * atoms[l].pos[m];
193	mmeineke	1067
194			totMass += atoms[l].mass;
195	mmeineke	1062	}
196			for(m=0;m<3;m++) com[m] /= totMass;
197
198			for(j=0;j<nLipAtoms;j++){
199			l=k+j;
200
201
202			distl = 0.0;
203			for(m=0; m<3;m++)
204			distl += (atoms[l].pos[m]-com[m])*(atoms[l].pos[m]-com[m]);
205
206			for(m=0; m<3; m++)
207			iMat[m][m] += atoms[l].mass * distl;
208
209			for(m=0;m<3;m++)
210			for(n=0;n<3;n++)
211			iMat[m][n] -= atoms[l].mass * ( atoms[l].pos[m] - com[m] ) *
212			( atoms[l].pos[n] - com[n] );
213
214
215			}
216	mmeineke	1067
217			// fprintf(stderr, "iMat (before) = %6G, %6G, %6G\n"
218			// " %6G, %6G, %6G\n"
219			// " %6G, %6G, %6G\n\n",
220			// iMat[0][0], iMat[0][1], iMat[0][2],
221			// iMat[1][0], iMat[1][1], iMat[1][2],
222			// iMat[2][0], iMat[2][1], iMat[2][2]);
223	mmeineke	1062
224	mmeineke	1067	// matWrap to fortran convention
225
226			for(j=0;j<3;j++)
227			for(l=0;l<3;l++)
228			matWrap[l+3*j] = iMat[l][j];
229
230	mmeineke	1062	ifail = 0;
231	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
232	mmeineke	1062
233			if (ifail) {
234			fprintf(stderr, "dsyev screwed something up!\n");
235			exit(0);
236			}
237
238	mmeineke	1067	// matWrap from fortran convention
239	mmeineke	1062
240	mmeineke	1067	for(j=0;j<3;j++)
241			for(l=0;l<3;l++)
242			iMat[l][j] = matWrap[l+3*j];
243	mmeineke	1062
244	mmeineke	1067	// fprintf(stderr, "evals = %6G, %6G, %6G\n\n",
245			// evals[0], evals[1], evals[2]);
246	mmeineke	1062
247	mmeineke	1067	// fprintf(stderr, "iMat (after) = %6G, %6G, %6G\n"
248			// " %6G, %6G, %6G\n"
249			// " %6G, %6G, %6G\n\n",
250			// iMat[0][0], iMat[0][1], iMat[0][2],
251			// iMat[1][0], iMat[1][1], iMat[1][2],
252			// iMat[2][0], iMat[2][1], iMat[2][2]);
253
254
255	mmeineke	1062
256			smallest = fabs(evals[0]);
257			which = 0;
258			for(j=0;j<3;j++){
259			if( fabs(evals[j]) < smallest ){
260			which = j;
261			smallest = fabs(evals[j]);
262			}
263			}
264
265			lenPrinc = 0.0;
266			for(j=0;j<3;j++)
267			lenPrinc += iMat[j][which] * iMat[j][which];
268	mmeineke	1067	lenPrinc = sqrt( lenPrinc );
269	mmeineke	1062
270			for(j=0;j<3;j++)
271			principal[i].u[j] = iMat[j][which]/lenPrinc;
272
273
274
275
276
277			// find the director of the bilayer
278
279			u = principal[i].u;
280
281			oMat[0][0] += u[0] * u[0] - oneThird;
282			oMat[0][1] += u[0] * u[1];
283			oMat[0][2] += u[0] * u[2];
284
285			oMat[1][0] += u[1] * u[0];
286			oMat[1][1] += u[1] * u[1] - oneThird;
287			oMat[1][2] += u[1] * u[2];
288
289			oMat[2][0] += u[2] * u[0];
290			oMat[2][1] += u[2] * u[1];
291			oMat[2][2] += u[2] * u[2] - oneThird;
292			}
293
294			for(i=0;i<3;i++)
295			for(j=0;j<3;j++)
296			oMat[i][j] /= (double)nLipids;
297
298	mmeineke	1067	// matWrap to fortran convention
299
300			for(j=0;j<3;j++)
301			for(l=0;l<3;l++)
302			matWrap[l+3*j] = oMat[l][j];
303
304	mmeineke	1062	ifail = 0;
305	mmeineke	1067	dsyev(&job, &uplo, &nfilled, matWrap, &ndiag, evals, work, &lWork, &ifail);
306	mmeineke	1062
307			if (ifail) {
308			fprintf(stderr, "dsyev screwed something up!\n");
309			exit(0);
310			}
311	mmeineke	1067
312			// matWrap from fortran convention
313	mmeineke	1062
314	mmeineke	1067	for(j=0;j<3;j++)
315			for(l=0;l<3;l++)
316			oMat[l][j] = matWrap[l+3*j];
317
318	mmeineke	1062	max = 0.0;
319			for (i=0; i<3;i++) {
320			if (fabs(evals[i]) > max) {
321			which = i;
322			max = fabs(evals[i]);
323			}
324			}
325
326			for (i = 0; i < 3; i++) {
327			directorWhole[index].u[i] = oMat[i][which];
328			}
329
330			directorWhole[index].order = 1.5 * max;
331
332	mmeineke	1067	// fprintf(stderr,
333			// "frame[%d] => order = %6G; < %6G, %6G, %6G >\n",
334			// index,
335			// directorWhole[index].order,
336			// directorWhole[index].u[0],
337			// directorWhole[index].u[1],
338			// directorWhole[index].u[2] );
339	mmeineke	1062
340			// fprintf(stdout,
341			// "%6g\t%6G\n",
342			// directorWhole[index].time,
343			// directorWhole[index].order);
344			}