ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/ripple/readh1.c
(Generate patch)

Comparing trunk/ripple/readh1.c (file contents):
Revision 1185 by xsun, Fri May 21 20:07:10 2004 UTC vs.
Revision 1333 by xsun, Fri Jul 16 18:05:12 2004 UTC

# Line 19 | Line 19 | inline double roundMe( double x ){
19   inline double roundMe( double x ){
20    return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
21   }
22 + inline double min( double a, double b ){
23 +  return (a < b ) ? a : b;
24 + }
25 + inline double max( double a, double b ){
26 +  return (a > b ) ? a : b;
27 + }
28  
29 +
30   // coords holds the data for a single tethered dipole:
31   struct coords{
32    double pos[3]; // cartesian coords
# Line 95 | Line 102 | int main(argc, argv)
102    double *x, *y, *z;
103    double dh2, dh, sumh2, sumh, averh2, averh, t, delta, gamma, hi, proj;
104    double *corrhist, *h2hist;
105 <  double vrhist[100][500];
106 <  double sum_vrhist[1000], ave_vrhist[1000];
107 <  int vrsamp[1000][1000];
105 >  double vrhist[1000];
106 >  double sum_vrhist[1000];
107 >  int vrsamp[1000];
108    int *ophist;
109    double d[2], hcorr;
110 <  double hsum, hsum_temp, h2sum, have, h2ave, fluc, bigL, smallA, areaPerMolecule, area, h, h2;
110 >  double hsum, hsum_frame, h2sum, have, h2ave, h_ave_frame;
111 >  double fluc, bigL, smallA, areaPerMolecule, area, h, h2;
112    int nbins, nbins2, opbin, whichbinx, whichbiny, whichbin2, n1, n2, n3, n4, m, selfx, selfy;
113  
114    int which;
115 +  int highestAtom;
116 +  double highestZ;
117    double omat[3][3];
118 +  double myPerp[3];
119 +  double myDir[3];
120 +  double myVec[2];
121 +  double lperp;
122 +  double ldir;
123 +  double dot;
124 +  double maxProj, maxProjOut;
125 +  double avgHeightAtProj;
126    double wrapMat[9];    
127    double onethird, ordvals[5000];
128 <  double max;
128 >  double maxEval;
129    double director[3][1000], vr[3][1000];
130    double sum_director[3], ave_director[3], sum_vr[3], ave_vr[3];
131    double orderpar[1000];
# Line 119 | Line 137 | int main(argc, argv)
137    int lwork;
138    double* work;
139    int ifail;
140 <  int done;
140 >  int done, lastData, firstData;
141    char current_flag;
142    
143    lineCount = 0;
# Line 162 | Line 180 | int main(argc, argv)
180    ifail = 0;
181    
182    nbins = 30;
183 <  nbins2 = 30;
183 >  nbins2 = 100;
184    binmin = 0.0;
185    binmax = 1.0;
186    delr = (binmax - binmin) / (double) nbins2;
# Line 185 | Line 203 | int main(argc, argv)
203    
204    for(i = 0; i < 1000; i++){
205      sum_vrhist[i] = 0.0;
188 –    ave_vrhist[i] = 0.0;
206    }
207  
208    for(i = 0; i < 3; i++){
# Line 202 | Line 219 | int main(argc, argv)
219      }
220    }
221  
222 <  for (i = 0; i < 100; i++) {
223 <    for (j = 0; j < 500; j++) {
207 <      vrhist[i][j] = 0.0;
208 <    }
222 >  for (i = 0; i < 1000; i++) {
223 >    vrhist[i] = 0.0;
224    }
225  
226    for(i = 0; i < nbins2; i++){
# Line 247 | Line 262 | int main(argc, argv)
262  
263    while(eof_test != NULL){
264      
265 +    highestAtom = -1;
266 +    highestZ = 0.0;
267 +
268      nframes++;
269      (void)sscanf(read_buffer, "%d", &state->nAtoms);
270      N = 2 * state->nAtoms;    
# Line 372 | Line 390 | int main(argc, argv)
390          exit(8);
391        }
392        (void)sscanf( foo, "%lf", &state->r[i].pos[2] );
393 +      if (state->r[i].pos[2] > highestZ) {
394 +        highestAtom = i;
395 +        highestZ = state->r[i].pos[2];
396 +      }
397        
398        foo = strtok(NULL, " ,;\t\0");
399        if(foo == NULL){
# Line 415 | Line 437 | int main(argc, argv)
437        state->r[i].mu = state->strength;
438      }
439      
440 <    hsum_temp = 0.0;
440 >    hsum_frame = 0.0;
441      for (i = 0; i < state->nAtoms; i++) {
442        
443        h = state->r[i].pos[2];
444        h2 = pow(h,2);
445 <
446 <      hsum_temp += h;
445 >      
446 >      hsum_frame += h;
447        hsum += h;
448        h2sum += h2;
449 <
449 >      
450        n1++;
451      }
452  
453 +    h_ave_frame = hsum_frame / (double) state->nAtoms;
454 +
455      for(i = 0; i < state->nAtoms; i++){
456      
457        omat[0][0] += ux[i] * ux[i] - onethird;
# Line 451 | Line 475 | int main(argc, argv)
475      //    temp_array = dsyev_ctof(omat, nfilled, nfilled);
476      
477      for(j=0;j<3;j++)
478 <            for(i=0;i<3;i++)
479 <                    wrapMat[i+j*3] = omat[i][j];
478 >      for(i=0;i<3;i++)
479 >        wrapMat[i+j*3] = omat[i][j];
480      
481      ifail = 0;
482      dsyev(&job, &uplo, &nfilled, wrapMat, &ndiag, evals, work, &lwork, &ifail);
483 <
483 >    
484      for(j=0;j<3;j++)
485 <            for(i=0;i<3;i++)
486 <                    omat[i][j] = wrapMat[i+j*3];
487 <                            
485 >      for(i=0;i<3;i++)
486 >        omat[i][j] = wrapMat[i+j*3];
487 >    
488      //dsyev_ftoc2(temp_array, omat, nfilled, nfilled);
489 <  
489 >    
490      //free(temp_array);
491    
492 <    max = 0.0;
492 >    maxEval = 0.0;
493      for(j = 0; j < 3; j++){
494 <      if(fabs(evals[j]) > max){
494 >      if(fabs(evals[j]) > maxEval){
495          which = j;
496 <        max = fabs(evals[j]);
496 >        maxEval = fabs(evals[j]);
497        }
498      }
475 –    
476 –    director[0][nframes-1] = omat[0][which];
477 –    director[1][nframes-1] = omat[1][which];
478 –    director[2][nframes-1] = omat[2][which];
479 –    
480 –    vr[0][nframes-1] = fabs(director[1][nframes-1]);
481 –    vr[1][nframes-1] = fabs(-director[0][nframes-1]);
482 –    vr[2][nframes-1] = 0;
483 –    
484 –    orderpar[nframes-1] = 1.5 * max;
499  
500 +    orderpar[nframes-1] = 1.5 * maxEval;    
501      opbin = (int) (orderpar[nframes-1] / delr);
502      if(opbin < nbins2) ophist[opbin] += 1;
503  
504 <    for(i = 0; i < state->nAtoms; i++){
505 <      proj = vr[0][nframes-1] * state->r[i].pos[0] + vr[1][nframes-1] * state->r[i].pos[1];
506 <      hi = state->r[i].pos[2]  - hsum_temp / (double) state->nAtoms;
504 >    for(i = 0; i < 3; i++){
505 >      myDir[i] = 0.0;
506 >      myPerp[i] = 0.0;
507 >    }
508 >
509 >    myDir[0] = omat[0][which];
510 >    myDir[1] = omat[1][which];
511 >    myDir[2] = omat[2][which];
512        
513 <      whichbin2 = (int) ((proj / (vr[0][nframes-1] * dx + vr[1][nframes-1] * dy)) * nbins2);
514 <      vrhist[whichbin2][nframes-1] += hi;
515 <      vrsamp[whichbin2][nframes-1]++;
513 >    if (myDir[0] < 0.0) {
514 >      myDir[0] = -myDir[0];
515 >      myDir[1] = -myDir[1];
516      }
517 +
518 +    ldir = sqrt(myDir[0]*myDir[0] + myDir[1]*myDir[1] + myDir[2]*myDir[2]);
519 +
520 +    myDir[0] /= ldir;
521 +    myDir[1] /= ldir;
522 +    myDir[2] /= ldir;
523 +
524 +    //printf("%f\t%f\t%f\n", myDir[0], myDir[1], myDir[2]);
525 +
526 +    director[0][nframes-1] = myDir[0];
527 +    director[1][nframes-1] = myDir[1];
528 +    director[2][nframes-1] = myDir[2];
529 +
530 +    //printf("%f\t%f\t%f\n", director[0][nframes-1], director[1][nframes-1], director[2][nframes-1]);
531 +
532 +    myPerp[0] = myDir[1];
533 +    myPerp[1] = -myDir[0];
534 +    myPerp[2] = 0.0;
535 +
536 +    if (myPerp[0] < 0.0) {
537 +      myPerp[0] = -myPerp[0];
538 +      myPerp[1] = -myPerp[1];
539 +    }
540 +
541 +    lperp = sqrt(myPerp[0]*myPerp[0] + myPerp[1]*myPerp[1] + myPerp[2]*myPerp[2]);
542 +
543 +    myPerp[0] /= lperp;
544 +    myPerp[1] /= lperp;
545 +    myPerp[2] /= lperp;
546 +
547 +    vr[0][nframes-1] = myPerp[0];
548 +    vr[1][nframes-1] = myPerp[1];
549 +    vr[2][nframes-1] = myPerp[2];
550 +
551 +    maxProj = 0.5 * sqrt(dx*dx + dy*dy);
552 +    //maxProj = myPerp[0]*dx + myPerp[1]*dy;
553 +        // for now, assume highest atom is atom 0
554 +    highestAtom = 0;
555 +    
556 +    for(i = 0; i < state->nAtoms; i++){
557 +
558 +      // difference vector from highest point:
559 +      myVec[0] = state->r[i].pos[0] - state->r[highestAtom].pos[0];
560 +      myVec[1] = state->r[i].pos[1] - state->r[highestAtom].pos[1];
561 +      // wrapped in periodic boundary conditions:
562 +      wrapVector(myVec, state->Hmat, state->HmatI);
563 +
564 +      // then projected onto myPerp:
565 +      proj = myPerp[0]*myVec[0] + myPerp[1]*myVec[1];
566 +      // and binned:
567 +
568 +      whichbin2 = (int)  (((1.0 + (proj/maxProj)) / 2.0) * (double) nbins2);
569 +      //  printf( "proj/maxProj = %f\twhichbin2 = %d\n", proj/maxProj, whichbin2);
570 +
571 +      hi = state->r[i].pos[2] - h_ave_frame;
572 +      //      for(i = 0; i < nbins2; i++) printf("%f\t%d\n", vrhist[i], vrsamp[i]);
573 +      vrhist[whichbin2] += hi;
574 +      vrsamp[whichbin2] ++;
575 +    }
576      
577      for(i = 0; i < state->nAtoms; i++){
578        for(j = 0; j < state->nAtoms; j++){
# Line 571 | Line 650 | int main(argc, argv)
650    ave_orderpar = sum_orderpar / (double) nframes;
651    ave2_orderpar = sum2_orderpar / (double) nframes;
652    err_orderpar = ave2_orderpar - pow(ave_orderpar, 2);
653 +
654 +  for(i = 0; i < 3; i++){
655 +    sum_director[i] = 0.0;
656 +    ave_director[i] = 0.0;
657 +    sum_vr[i] = 0.0;
658 +    ave_vr[i] = 0.0;
659 +  }
660    
661    for(i = 0; i < nframes; i++){
662      sum_director[0] += director[0][i];
# Line 594 | Line 680 | int main(argc, argv)
680    printf("# error = %lf\n", err_orderpar);
681    printf("# director axis is ( %lf\t%lf\t%lf )\n",
682                    ave_director[0], ave_director[1], ave_director[2]);  
683 <  printf("# vr axis is ( %lf\t%lf\t%lf )\n", ave_vr[0], ave_vr[1], ave_vr[2]);  
683 >  printf("# vr axis is ( %lf\t%lf\t%lf )\n", ave_vr[0], ave_vr[1], ave_vr[2]);
684  
599 –  for(i = 0; i < nbins2; i++){
600 –    for(j = 0; j < nframes; j++){
601 –      sum_vrhist[i] += vrhist[i][j];
602 –    }
603 –  }
604 –
685    dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
686    dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));
687    
# Line 617 | Line 697 | int main(argc, argv)
697  
698    printf("# first gamma estimate = %lf\n", gamma);
699  
700 <  printf("# \n");
700 >  firstData = -1;
701 >  lastData = -1;
702 >
703 >  for (i=0; i< nbins2; i++) {
704 >    if (vrsamp[i] > 0) {
705 >      if (firstData == -1) firstData = i;
706 >      lastData = i;
707 >    }
708 >  }    
709 >
710 >  maxProj = 0.5*sqrt(dx*dx + dy*dy);
711 >  maxProjOut = 2.0*maxProj*(double)(lastData - firstData)/(double)nbins2;
712    
713 <  for(i = 0; i < nbins2; i++){
714 <    ave_vrhist[i] = sum_vrhist[i] / (double) nframes;
715 <    printf("%lf\t%lf\n", (double) i * (ave_vr[0] * dx + ave_vr[1] * dy) / (double) nbins2,
716 <                   ave_vrhist[i]);
713 >  printf("# maximum projection = %lf\n", maxProjOut);
714 >  printf("# \n");
715 >
716 >  for (i = firstData ; i < lastData; i++) {
717 >    //proj = maxProj * ((double)i / (double)nbins2);      
718 >    
719 >    proj = maxProj*( (2.0 * (double)i / (double)nbins2) -1.0);
720 >    
721 >    if (vrsamp[i] > 0) {
722 >      avgHeightAtProj = vrhist[i]/(double)vrsamp[i];
723 >    } else {
724 >      avgHeightAtProj = 0.0;
725 >    }
726 >    printf("%lf\t%lf\n", proj, avgHeightAtProj);
727    }
728    
729    selfx = (int) ((double)nbins / 2.0);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines