trunk/ripple/readh1.c

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<math.h>
#include<fftw.h>
#include<mkl_lapack64.h>

//extern void dsyev(char *jobz, char *uplo, int *n, double *a, int *lda,
//              double *w, double *work, int *lwork,int *info);

//void direct(double rcut, double box, int n, int nstep, int maxbin);

//double* dsyev_ctof(double **in, int rows, int cols);

//void dsyev_ftoc2(double *in, double **out, int rows, int cols);

// Structures to store our data:

inline double roundMe( double x ){
  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
}

// coords holds the data for a single tethered dipole:
struct coords{
  double pos[3]; // cartesian coords
  double theta;  // orientational angle relative to z axis
  double phi;    // orientational angle in x-y plane
  double mu;     // dipole strength
  char name[30]; // an identifier for the type of atom
};

// state holds the current "configuration" of the entire system
struct system {
  int    nAtoms;           // Number of Atoms in this configuration
  struct coords *r;        // The set of coordinates for all atoms
  double beta;             // beta = 1 /(kb*T)
  double strength;         // strength of the dipoles (Debye)
  double z0;               // default z axis position
  double theta0;           // default theta angle
  double kz;               // force constant for z displacement
  double ktheta;           // force constant for theta displacement
  int    nCycles;          // How many cycles to do in total
  int    iCycle;           // How many cycles have we done?
  int    nMoves;           // How many MC moves in each cycle
  int    nSample;          // How many cycles between samples
  double Hmat[2][2];       // The information about the size of the per. box
  double HmatI[2][2];      // The inverse box
  double energy;           // The current Energy
  double dtheta;           // maximum size of a theta move
  double deltaz;           // maximum size of a z move
  double deltaphi;         // maximum size of a phi move
  int    nAttempts;        // number of MC moves that have been attempted
  int    nAccepts;         // number of MC moves that have been accepted
  int    nx;               // number of unit cells in x direction
  int    ny;               // number of unit cells in y direction
  struct system *next;      // Next frame in the linked list
};

char *program_name; /* the name of the program */

// Function prototypes:
void usage(void);
void invertMat2(double a[2][2], double b[2][2]);
void wrapVector( double thePos[2], double Hmat[2][2], double HmatI[2][2]);

int main(argc, argv)
     int argc;
     char *argv[];
{
  FILE *in_file;
  char in_name[500];
  char *eof_test, *foo;
  char read_buffer[1000];
  //int lineCount = 0;
  int lineCount;
  int nAtoms;
  double *mag, *newmag;
  int *present_in_old;
  double *ux, *uy, *uz, p1;
  double aLat, bLat;
  int cells;
  double sumZ, sumUx, sumUy, sumUz, sumP;
  double interpsum, value;
  int ninterp, px, py, newp;
  int i, j, k, l, nloops;
  int newx, newy, newindex, index;
  int new_i, new_j, new_index;
  int N, nframes;
  double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
  double maxfreqx, maxfreqy, maxfreq, dfreq;
  double dx, dy, dx1, dy1, xTemp, yTemp, pt1x, pt1y, pt2x, pt2y;
  int nx, ny;
  int *samples;
  double *bin, binmin, binmax, delr;
  double *x, *y, *z;
  double dh2, dh, sumh2, sumh, averh2, averh, t, delta, gamma, hi, proj;
  double *corrhist, *h2hist;
  double vrhist[100][500];
  double sum_vrhist[1000], ave_vrhist[1000];
  int vrsamp[1000][1000];
  int *ophist;
  double d[2], hcorr;
  double hsum, hsum_temp, h2sum, have, h2ave, fluc, bigL, smallA, areaPerMolecule, area, h, h2;
  int nbins, nbins2, opbin, whichbinx, whichbiny, whichbin2, n1, n2, n3, n4, m, selfx, selfy;

  int which;
  double omat[3][3];
  double wrapMat[9];    
  double onethird, ordvals[5000];
  double max;
  double director[3][1000], vr[3][1000];
  double sum_director[3], ave_director[3], sum_vr[3], ave_vr[3];
  double orderpar[1000];
  double sum_orderpar, sum2_orderpar, ave_orderpar, ave2_orderpar, err_orderpar;
  char job, uplo;
  int ndiag;
  int nfilled;
  double evals[100];
  int lwork;
  double* work;
  int ifail;
  int done;
  char current_flag;
  
  lineCount = 0;
  
  program_name = argv[0];
  if (argc >= 2)
    strcpy(in_name, argv[1]);
/*  
  for(i = 1; i < argc; i++){
    if(argv[i][0] == '-'){
      done = 0;
      j = 1;
      current_flag = argv[i][j];
      while( (current_flag != '\0') && ( !done ) ){
        switch(current_flag){
          case 'i':
            i++;
            strcpy( in_name, argv[i] );
            done = 1;
            break;
        }
      }
    }
  }
*/

  struct system* state;
  struct system* temp_state;
  struct coords* r;

  lwork = 9;

  work = (double *) malloc(lwork * sizeof(double));

  onethird = 1.0 / 3.0;
  ndiag = 3;
  nfilled = 3;
  job = 'V';
  uplo = 'U';
  ifail = 0;
  
  nbins = 30;
  nbins2 = 30;
  binmin = 0.0;
  binmax = 1.0;
  delr = (binmax - binmin) / (double) nbins2;
  corrhist = (double *) calloc(nbins*nbins, sizeof(double));
  h2hist = (double *) calloc(nbins*nbins, sizeof(double));
  ophist = (int *) calloc(nbins2, sizeof(int));
  hsum = 0.0;
  h2sum = 0.0;
  sum_orderpar = 0.0;
  sum2_orderpar = 0.0;
  ave_orderpar = 0.0;
  ave2_orderpar = 0.0;

  for(i = 0; i < 3; i++){
    for(j = 0; j < 1000; j++){
      director[i][j] = 0.0;
      vr[i][j] = 0.0;
    }
  }
  
  for(i = 0; i < 1000; i++){
    sum_vrhist[i] = 0.0;
    ave_vrhist[i] = 0.0;
  }

  for(i = 0; i < 3; i++){
    sum_director[i] = 0.0;
    ave_director[i] = 0.0;
    sum_vr[i] = 0.0;
    ave_vr[i] = 0.0;
  }

  for (i = 0; i < nbins; i++) {
    for (j = 0; j < nbins; j++) {
      corrhist[nbins * i + j] = 0.0;
      h2hist[nbins * i + j] = 0.0;
    }
  }

  for (i = 0; i < 100; i++) {
    for (j = 0; j < 500; j++) {
      vrhist[i][j] = 0.0;
    }
  }

  for(i = 0; i < nbins2; i++){
    ophist[i] = 0;
  }

  t = 300;

  in_file = fopen(in_name, "r");
  if(in_file == NULL){
    printf("Cannot open file \"%s\" for reading.\n", in_name);
    exit(8);
  }

  nframes = 0;
  n1 = 0;
  n2 = 0;
  n3 = 0;
  n4 = 0;
   
  // start reading the first frame

  eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
  nAtoms = atoi(read_buffer);
  ux = (double *) calloc(nAtoms, sizeof(double));
  uy = (double *) calloc(nAtoms, sizeof(double));
  uz = (double *) calloc(nAtoms, sizeof(double));
  lineCount++;

  for(i = 0; i < nAtoms; i++){
    ux[i] = 0.0; 
    uy[i] = 0.0;
    uz[i] = 0.0;
  }
  
  state = (struct system *) malloc(sizeof(struct system));
  state->next = NULL;
  state->strength = 7.0;

  while(eof_test != NULL){
    
    nframes++;
    (void)sscanf(read_buffer, "%d", &state->nAtoms);
    N = 2 * state->nAtoms;    

    state->r = (struct coords *)calloc(N, sizeof(struct coords));

    for(i = 0; i < 3; i++){
      for(j = 0; j < 3; j++){
        omat[i][j] = 0.0;
      }
    }

    // read and the comment line and grab the time and box dimensions
    
    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;
    if(eof_test == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    
    foo = strtok( read_buffer, " ,;\t\n" );
    (void)sscanf( read_buffer, "%d", &state->iCycle );

    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%d", &state->nx);

    nx = state->nx;

    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%d", &state->ny);
   
    ny = state->ny;

    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[0][0]);
    
    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[1][0]);
    
    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[0][1]);
    
    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[1][1]);

    //Find HmatI:

    invertMat2(state->Hmat, state->HmatI);
    
    // Length of the two box vectors:

    dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
    dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));

    aLat = dx / (double)(state->nx);
    bLat = dy / (double)(state->ny);
 
    // FFT stuff depends on nx and ny, so delay allocation until we have
    // that information
    
    for (i=0;i<state->nAtoms;i++){
      
      eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
      lineCount++;
      if(eof_test == NULL){
        printf("error in reading file at line: %d\n", lineCount);
        exit(8);
      }
      
      foo = strtok(read_buffer, " ,;\t\0");
      (void)strcpy(state->r[i].name, foo); //copy the atom name
      
      // next we grab the positions
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading postition x from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &state->r[i].pos[0] );
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading postition y from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &state->r[i].pos[1] );
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading postition z from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &state->r[i].pos[2] );
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading angle phi from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &state->r[i].phi );
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading unit vector x from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &ux[i] );
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading unit vector y from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &uy[i] );
      
      foo = strtok(NULL, " ,;\t\0");
      if(foo == NULL){
        printf("error in reading unit vector z from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, state->nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &uz[i] );
      
      state->r[i].theta = acos(uz[i]);
      
      // The one parameter not stored in the dump file is the dipole strength
      state->r[i].mu = state->strength;
    }
    
    hsum_temp = 0.0;
    for (i = 0; i < state->nAtoms; i++) {
      
      h = state->r[i].pos[2];
      h2 = pow(h,2);

      hsum_temp += h;
      hsum += h;
      h2sum += h2;

      n1++;
    }

    for(i = 0; i < state->nAtoms; i++){
    
      omat[0][0] += ux[i] * ux[i] - onethird;
      omat[0][1] += ux[i] * uy[i];
      omat[0][2] += ux[i] * uz[i];
      omat[1][0] += uy[i] * ux[i];
      omat[1][1] += uy[i] * uy[i] - onethird;
      omat[1][2] += uy[i] * uz[i];
      omat[2][0] += uz[i] * ux[i];
      omat[2][1] += uz[i] * uy[i];
      omat[2][2] += uz[i] * uz[i] - onethird;

    }

    for(i = 0; i < 3; i++){
      for(j = 0; j < 3; j++){
        omat[i][j] /= (double) state->nAtoms;
      }
    }
    
    //    temp_array = dsyev_ctof(omat, nfilled, nfilled);
    
    for(j=0;j<3;j++)
            for(i=0;i<3;i++)
                    wrapMat[i+j*3] = omat[i][j];
    
    ifail = 0;
    dsyev(&job, &uplo, &nfilled, wrapMat, &ndiag, evals, work, &lwork, &ifail);

    for(j=0;j<3;j++)
            for(i=0;i<3;i++)
                    omat[i][j] = wrapMat[i+j*3];
                            
    //dsyev_ftoc2(temp_array, omat, nfilled, nfilled);
  
    //free(temp_array);
  
    max = 0.0;
    for(j = 0; j < 3; j++){
      if(fabs(evals[j]) > max){
        which = j;
        max = fabs(evals[j]);
      }
    }
    
    director[0][nframes-1] = omat[0][which];
    director[1][nframes-1] = omat[1][which];
    director[2][nframes-1] = omat[2][which];
    
    vr[0][nframes-1] = fabs(director[1][nframes-1]);
    vr[1][nframes-1] = fabs(-director[0][nframes-1]);
    vr[2][nframes-1] = 0;
    
    orderpar[nframes-1] = 1.5 * max;

    opbin = (int) (orderpar[nframes-1] / delr);
    if(opbin < nbins2) ophist[opbin] += 1;

    for(i = 0; i < state->nAtoms; i++){
      proj = vr[0][nframes-1] * state->r[i].pos[0] + vr[1][nframes-1] * state->r[i].pos[1];
      hi = state->r[i].pos[2]  - hsum_temp / (double) state->nAtoms;
      
      whichbin2 = (int) ((proj / (vr[0][nframes-1] * dx + vr[1][nframes-1] * dy)) * nbins2);
      vrhist[whichbin2][nframes-1] += hi;
      vrsamp[whichbin2][nframes-1]++;
    }
    
    for(i = 0; i < state->nAtoms; i++){
      for(j = 0; j < state->nAtoms; j++){
        
        d[0] = state->r[j].pos[0] - state->r[i].pos[0];
        d[1] = state->r[j].pos[1] - state->r[i].pos[1];
        
        wrapVector(d, state->Hmat, state->HmatI);
        
        whichbinx = (int) ((nbins-1) * (dx/2.0 + d[0]) / dx);
        whichbiny = (int) ((nbins-1) * (dy/2.0 + d[1]) / dy);

        //if (i == j) {
          //printf("whichbinx = %i, whichbiny = %i\n", whichbinx, whichbiny);
        //}

        //printf("d0 = %lf, d1 = %lf\n", d[0], d[1]);
        //printf("wx = %d, wy = %d\n", whichbinx, whichbiny);

        if (whichbinx >= nbins || whichbiny >= nbins) {
          printf("off by one error\n");
          printf("whichbinx = %i, whichbiny = %i\n", whichbinx, whichbiny);
          exit(0);
        }
        if (whichbinx < 0 || whichbiny < 0) {
          printf("off by one error\n");
          printf("whichbinx = %i, whichbiny = %i\n", whichbinx, whichbiny);
          exit(0);
        }

        hcorr = state->r[j].pos[2] * state->r[i].pos[2];

        corrhist[nbins * whichbinx + whichbiny] += hcorr;

        dh = state->r[j].pos[2] - state->r[i].pos[2];
        dh2 = pow(dh, 2);

        h2hist[nbins * whichbinx + whichbiny] += dh2;

        n2++;
      }
    }

    temp_state = state->next;
    state->next = NULL;
  
    if (temp_state != NULL) {
      free(temp_state->r);
      free(temp_state);
    }

    // Make a new frame
  
    temp_state = (struct system *) malloc(sizeof(struct system));
    temp_state->next = state;
    state = temp_state;
    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;
  }


  have = hsum / (double) n1;
  h2ave = h2sum / (double) n1;
  fluc = h2ave - pow(have, 2);

  printf("# <h> = %lf\n", have);
  printf("# <h2> = %lf\n", h2ave);
  printf("# sigma(h) = %lf\n", sqrt(h2ave - have * have));
  printf("# fluctuation = %lf\n", fluc);

  for(i = 0; i < nframes; i++){
    sum_orderpar += orderpar[i];
    sum2_orderpar += pow(orderpar[i], 2);
  }
  ave_orderpar = sum_orderpar / (double) nframes;
  ave2_orderpar = sum2_orderpar / (double) nframes;
  err_orderpar = ave2_orderpar - pow(ave_orderpar, 2);
  
  for(i = 0; i < nframes; i++){
    sum_director[0] += director[0][i];
    sum_director[1] += director[1][i];
    sum_director[2] += director[2][i];
  }
  ave_director[0] = sum_director[0] / (double) nframes;
  ave_director[1] = sum_director[1] / (double) nframes;
  ave_director[2] = sum_director[2] / (double) nframes;

  for(i = 0; i < nframes; i++){
    sum_vr[0] += vr[0][i];
    sum_vr[1] += vr[1][i];
    sum_vr[2] += vr[2][i];
  }
  ave_vr[0] = sum_vr[0] / (double) nframes;
  ave_vr[1] = sum_vr[1] / (double) nframes;
  ave_vr[2] = sum_vr[2] / (double) nframes;

  printf("# orderparameter = %lf\n", ave_orderpar);
  printf("# error = %lf\n", err_orderpar);
  printf("# director axis is ( %lf\t%lf\t%lf )\n",
                  ave_director[0], ave_director[1], ave_director[2]);  
  printf("# vr axis is ( %lf\t%lf\t%lf )\n", ave_vr[0], ave_vr[1], ave_vr[2]);  

  for(i = 0; i < nbins2; i++){
    for(j = 0; j < nframes; j++){
      sum_vrhist[i] += vrhist[i][j];
    }
  }

  dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
  dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));
  
  area = dx * dy;

  bigL = sqrt(area);

  areaPerMolecule = area / (double) state->nAtoms;

  smallA = sqrt(areaPerMolecule);

  gamma = t * log(bigL / smallA)  / (2.0 * M_PI * fluc);

  printf("# first gamma estimate = %lf\n", gamma);

  printf("# \n");
  
  for(i = 0; i < nbins2; i++){
    ave_vrhist[i] = sum_vrhist[i] / (double) nframes;
    printf("%lf\t%lf\n", (double) i * (ave_vr[0] * dx + ave_vr[1] * dy) / (double) nbins2,
                   ave_vrhist[i]);
  }
  
  selfx = (int) ((double)nbins / 2.0);
  selfy = (int) ((double)nbins / 2.0);

/*   for (i = 0; i < nbins; i++) { */
/*     for (j = 0; j < nbins; j++) { */
/*       printf("%lf\t", h2hist[nbins * i + j] / (double) n2); */
/*     } */
/*     printf("\n"); */
/*   } */

  free(work);
  free(corrhist);
  free(h2hist);
  free(ux);
  free(uy);
  free(uz);
  return 1;

}

double matDet2(double a[2][2]) {

  double determinant;

  determinant = (a[0][0] * a[1][1]) - (a[0][1] * a[1][0]);

  return determinant;
}


void invertMat2(double a[2][2], double b[2][2]) {

  double determinant;

  determinant = matDet2( a );

  if (determinant == 0.0) {
    printf("Can't invert a matrix with a zero determinant!\n");
  }

  b[0][0] =  a[1][1] / determinant;
  b[0][1] = -a[0][1] / determinant;
  b[1][0] = -a[1][0] / determinant;
  b[1][1] =  a[0][0] / determinant;
}

void matVecMul2(double m[2][2], double inVec[2], double outVec[2]) {
  double a0, a1, a2;

  a0 = inVec[0];  a1 = inVec[1];

  outVec[0] = m[0][0]*a0 + m[0][1]*a1;
  outVec[1] = m[1][0]*a0 + m[1][1]*a1;
}

void wrapVector( double thePos[2], double Hmat[2][2], double HmatInv[2][2]){

  int i;
  double scaled[2];

  // calc the scaled coordinates.

  matVecMul2(HmatInv, thePos, scaled);
  
  for(i=0; i<2; i++)
    scaled[i] -= roundMe(scaled[i]);
  
  // calc the wrapped real coordinates from the wrapped scaled coordinates
  
  matVecMul2(Hmat, scaled, thePos);
  
}

/* double* dsyev_ctof(double **in, int rows, int cols) */
/* { */
/*   double *out; */
/*   int i, j; */

/*   out = (double *) calloc(rows * cols,  sizeof(double)); */
  
/*   if (!out){ */
/*     printf("Fail to allocate memory\n"); */
/*     exit(1);  */
/*   }  */

/*   for (i = 0; i < rows; i++) */
/*     for (j = 0; j < cols; j++) */
/*       out[i+j*cols] = in[i][j]; */
  
/*   return(out); */
/* } */

/* void dsyev_ftoc2(double *in, double **out, int rows, int cols) */
/* { */
/*   int i, j; */

/*   for (i = 0; i < rows; i++) */
/*     for (j = 0; j < cols; j++) */
/*       out[i][j] = in[i+j*cols]; */
/* } */
/* void direct(double rcut, double box, int n, int nstep, int maxbin){ */

/*   int bin, maxbin, nbins, i, j, k, n, bind, which; */
/*   int hist[maxbin], histd[maxbin], startstep; */
/*   int nstep, stopstep, ndiag, nfilled, ifail, lwork; */
/*   double delr, rijsp, rxij, ryij, rzij, rij; */
/*   double delrd, rijsqd, rijd, orderpar[1000], nideal; */
/*   double omat[3][3], evals[100], cons, rlower; */
/*   double rcut, box, dens, director[3][1000], rupper; */
/*   double rx[1000][1000], ry[1000][1000], rz[1000][1000]; */
/*   double ux[1000][1000], uy[1000][1000], uz[1000][1000]; */
/*   double vx[1000][1000], vy[1000][1000], vz[1000][1000]; */
/*   double jx[1000][1000], jy[1000][1000], jz[1000][1000]; */
/*   double grd[5000], gr[5000], max, binmin, binmax; */
/*   double onethird, ordvals[5000]; */
/*   char job, uplo; */
/*   int ndiag; */
/*   int nfilled; */
/*   double omat[3][3]; */
/*   double evals[100]; */
/*   double work[]; */
/*   int lwork; */
/*   int ifail; */

/*   lwork = 9; */
/*   work = (double *)calloc(lwork, sizeof(double)); */

/*   onethird = 1.0 / 3.0; */
/*   ndiag = 3; */
/*   nfilled = 3; */
/*   job = 'V'; */
/*   uplo = 'U'; */
/*   ifail = 0; */

/*   binmin = 0.0; */
/*   binmax = 1.0; */
/*   delr = (binmax - binmin) / (double) maxbin; */

/*   for(i = 0; i < maxbin; i++){ */

/*     ordvals(i) = 0.0; */
/*     hist(i) = 0.0; */
    
/*   } */

/*   for(i = 0; i < stopstep; i++){ */

/*     for(j = 0; j < 3; j++){ */
/*       for(k = 0; k < 3; k++){ */
/*         omat[j][k] = 0.0; */
/*       } */
/*     } */
    
/*     for(j = 0; j < n; j++){ */
    
/*       omat[1][1] = omat[1][1] + ux[j][i] * ux[j][i] - onethird; */
/*       omat[1][2] = omat[1][2] + ux[j][i] * uy[j][i]; */
/*       omat[1][3] = omat[1][3] + ux[j][i] * uz[j][i]; */
/*       omat[2][1] = omat[2][1] + uy[j][i] * ux[j][i]; */
/*       omat[2][2] = omat[2][2] + uy[j][i] * uy[j][i] - onethird; */
/*       omat[2][3] = omat[2][3] + uy[j][i] * uz[j][i]; */
/*       omat[3][1] = omat[3][1] + uz[j][i] * ux[j][i]; */
/*       omat[3][2] = omat[3][2] + uz[j][i] * uy[j][i]; */
/*       omat[3][3] = omat[3][3] + uz[j][i] * uz[j][i] - onethird; */
    
/*     } */

  
/*     for(j = 0; j < 3; j++){ */
/*       for(k = 0; k < 3; k++){ */
/*         omat[j][k] = omat[j][k] / (double) n; */
/*       } */
/*     } */

/*     dsyev_(job, uplo, nfilled, omat, ndiag, evals, work, lwork, ifail); */

/*     max = 0.0; */
/*     for(j = 0; j < 3; j++){ */
/*       if(fabs(evals(j)) > max){ */
/*         which = j; */
/*         max = fabs(evals(j)); */
/*       } */
/*     } */
    
/*     director[1][i] = omat[1][which]; */
/*     director[2][i] = omat[2][which]; */
/*     director[3][i] = omat[3][which]; */

/*     orderpar(i) = 1.5 * max; */

/*     bin = (int) (orderpar(i) / delr) + 1; */

/*     if(bin < maxbin) */
/*       hist(bin) = hist(bin) + 1; */
/*   } */
  
/*   cons = 1.0; */
/*   for(bin = 0; bin < maxbin; bin++){ */
    
/*     rlower = (double) (bin - 1) * delr; */
/*     rupper = rlower + delr; */
/*     ordvals(bin) = rlower + (delr / 2.0); */
/*     nideal = cons * (pow(rupper, 3) - pow(rlower, 3)); */
/*     gr(bin) = (double) (hist(bin)) / (double) nstep /nideal; */
    
/*   } */
/* } */
Revision:	1185
Committed:	Fri May 21 20:07:10 2004 UTC (20 years, 3 months ago) by xsun
Content type:	text/plain
File size:	22297 byte(s)
Log Message:	this is the ripple codes