trunk/ripple/dp.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <time.h>
#include "mkl_vsl.h"

// inlined functions (for periodic box wrapping)

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

// Structures to store our data:

// 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 kr;               // force constant for z displacement
  double ktheta;           // force constant for theta displacement
  double t;                // the temperature of the system
  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
  double XYNNDIST;         // maximum distance of nearest neighbors in XY plane
};

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

// Function prototypes:
void usage(void);
double toterg(struct system* state);
void getmag(struct system* state, double mag[3]);
double eneri(struct system* state, struct coords iTemp,  int i, int jb);
void adjust(struct system* state);
void mcmove(struct system* state, VSLStreamStatePtr stream, double *en);
void store(struct system* state, FILE* out_file);
void invertMat2(double a[2][2], double b[2][2]);
void wrapVector( double thePos[2], double Hmat[2][2], double HmatI[2][2]);

// Defines for the MKL random Number generator:
#define SEED    1
#define BRNG    VSL_BRNG_MCG31
#define METHOD  0
#define N       1

int main(argc, argv) 
     int argc;
     char *argv[];
{

  FILE *in_file; /* the input file */
  FILE *out_file; /* the output file */
  char *out_prefix = NULL; /* the prefix of the output file 
                              (if different from the root_name) */
  char out_name[500]; /* the output name */
  char in_name[500]; // the input file name
  char *root_name = NULL; /* the root name */

  int have_outName = 0;
  int have_inName = 0;
  int have_rootName = 0;
  int restart_from_file = 0;

  int i, j, k, imove;

  char temp_name[500];
  int lineCount = 0;  // the line number
  char read_buffer[1000]; /* the line buffer for reading */
  char *eof_test; /* ptr to see when we reach the end of the file */
  char *foo; /* the pointer to the current string token */
  char atomName[10];

  double aLat;
  double bLat;
  int cells;

  int haveAlat = 0;
  int haveBlat = 0;
  int haveCells = 0;

  struct system* state;
  struct coords* r;

  int done;
  char current_flag;

  double dx, dy, uxi, uyi, uzi, myran, en, ent, magmag;
  double mag[3];
  int which;

  // Other useful defines:

  double twopi = 2.0 * M_PI;
  double kb = 0.0019872198;
  strcpy(atomName, "Ar");

  state = (struct system *)malloc(sizeof(struct system));

  // The parameters for the simulation:

  state->strength = 7.0;
  state->z0 = 0.0;
  state->kz = kb; 
  state->kr = kb;
  state->theta0 = M_PI / 2.0;
  state->ktheta = 0.0; 
  state->dtheta = 0.3;
  state->deltaz = 1.0;
  state->deltaphi = 0.5; 
  state->t = 300.0;
  state->beta = 1.0 / (kb * state->t);

  // these three should really be given as command line arguments, but we'll
  // set defaults here just in case:

  state->nCycles = 1000000;
  state->nMoves  = 100;

  // Store roughly 100 frames while we are testing
  state->nSample = state->nCycles / 100;

  // Stuff for initializing the random number generator:
  struct        timeval         now_time_val;
  struct        timezone        time_zone;
  struct        tm              *now_tm;
  time_t                        now;
  int                           mySeed;
  VSLStreamStatePtr             stream;

  // First, initialize the random number generator:
  gettimeofday(&now_time_val, &time_zone);  /* get the time now  */
  now = now_time_val.tv_sec;                /*  convert to epoch time */
  mySeed = (int) now;
  vslNewStream(&stream, BRNG, mySeed);

  // Now handle the arguments to the program:

  program_name = argv[0]; // save the program name in case we need it

  for( i = 1; i < argc; i++){

    if(argv[i][0] =='-'){
      
      // parse single character options
      
      done =0;
      j = 1;
      current_flag = argv[i][j];
      while( (current_flag != '\0') && (!done) ){
        
        switch(current_flag){
          
        case 'o':
          // -o <filename> => the output
          
          i++;
          strcpy( out_name, argv[i] );
          have_outName = 1;
          done = 1;
          break;
          
        case 'i':
          // -i <filename> => the input
          
          i++;
          strcpy( in_name, argv[i] );
          have_inName = 1;
          done = 1;
          break;

        case 'r':
          // -r <root> => root Name
          
          i++;
          strcpy( root_name, argv[i] );
          have_rootName = 1;
          done = 1;
          break;
          
        case 'n':
          // -n <#> => do Cycle MC cycles
          
          i++;
          state->nCycles = atoi(argv[i]);
          done = 1;
          break;

        case 's':
          // -s <#> => write out every nSample MC cycles
          printf("setting nSample\n"); 
          i++;
          state->nSample = atoi(argv[i]);
          printf("setting nSample to %d\n", state->nSample); 
          done = 1;
          break;

        case 'm':
          // -m <#> => each MC cycle consists of nMoves atomic moves
          
          i++;
          state->nMoves = atoi(argv[i]);
          done = 1;
          break;

        case 'k':
          // -k kz   sets the value of kz in units of kb
          i++;
          state->kz = kb * atof( argv[i] );
          state->kr = state->kz;
          done = 1;
          break;
          
        case 'q':
          // -q mu set the strength of the dipole
          i++;
          state->strength =  atof( argv[i] );
          done = 1;
          break;
 
        case 't':
          // -t set the temperature of the system
          i++;
          state->t =  atof( argv[i] );
          state->beta = 1.0 / (kb * state->t);
          done = 1;
          break;

        case 'h':
          // -h hexSpace  set the inter-atomic spacing for a regular
          //              hexagonal lattice
          haveAlat = 1;
          haveBlat = 1;
          i++;
          bLat = atof( argv[i] );
          aLat = sqrt(3.0) * bLat;
          done = 1;
          break;

        case 'a':
          // -a aLat  set the x lattice spacing for the distorted hexagonal 
          //          lattice  
          haveAlat = 1;
          i++;
          aLat = atof( argv[i] );
          done = 1;
          break;
          
        case 'b':
          // -b bLat  set the y lattice spacing for the distorted hexagonal 
          //          lattice when initializing the system
          haveBlat = 1;
          i++;
          bLat = atof( argv[i] );
          done = 1;
          break;

        case 'c':
          // -c cells  sets the number of unit cells along the x axis
          //           to use when initializing the system
          haveCells = 1;
          i++;
          cells = atoi( argv[i] );
          done = 1;
          break;

        case 'x':
          usage();
          done = 1;
          break;
          
        default:
          
          (void)fprintf(stderr, "Bad option \"-%s\"\n", current_flag);
          usage();
        }
        j++;
        current_flag = argv[i][j];
      }
    } else {
      
      if( root_name != NULL ){
        fprintf( stderr,
                 "The root name has already been set with an argument to -r\n"
                 "But another argument looks like the root name.  Whassup?\n");
        usage();
      }
      
      root_name = argv[i];
    }
  }
  
  if ( !have_outName && root_name == NULL) {
    fprintf( stderr, "No output or root filename was set, so look for your data in dp.dump!\n");   
    strcpy(out_name, "dp.dump");
    have_outName=1;
  }
  

  // Figure out if we are starting from an input file:

  if (have_inName) {

    restart_from_file = 1;
    
    in_file = fopen(in_name, "r");
    if(in_file == NULL){
      printf("Cannot open file \"%s\" for reading.\n", in_name);
      exit(8);
    }
    
    // start reading the first frame
    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;    

    while(eof_test != NULL){
      
      (void)sscanf(read_buffer, "%d", &state->nAtoms);
      state->r =
        (struct coords *)calloc(state->nAtoms, sizeof(struct coords));
      
      // 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" );
      (void)sscanf( read_buffer, "%d", &state->iCycle );
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading file at line: %d\n", lineCount);
        exit(8);
      }
      (void)sscanf(foo, "%d", &state->nx);

      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading file at line: %d\n", lineCount);
        exit(8);
      }
      (void)sscanf(foo, "%d", &state->ny);
     
      foo = strtok(NULL, " ,;\t");
      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");
      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");
      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");
      if(foo == NULL){
        printf("error in reading file at line: %d\n", lineCount);
        exit(8);
      }
      (void)sscanf(foo, "%lf",&state->Hmat[1][1]);
      
      // 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);
      
      if (0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)) > bLat) {
         state->XYNNDIST =  0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0));
      } else {
         state->XYNNDIST = bLat;
      }
      // slightly larger so we can use <  as a comparison 
      state->XYNNDIST = state->XYNNDIST * 1.01;

      // Find HmatI:
      
      invertMat2(state->Hmat, state->HmatI);
      
      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");
        (void)strcpy(state->r[i].name, foo); //copy the atom name
        
        // next we grab the positions
        
        foo = strtok(NULL, " ,;\t");
        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");
        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");
        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");
        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");
        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", &uxi );
        
        foo = strtok(NULL, " ,;\t");
        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", &uyi );
        
        foo = strtok(NULL, " ,;\t");
        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", &uzi );
        
        state->r[i].theta = acos(uzi);

        // The one parameter not stored in the dump file is the dipole strength
        state->r[i].mu = state->strength;

      }
      eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
      lineCount++;
    } 
    (void)fclose(in_file);
    
  } else {

    // not restarting from file, so use data we've got!

    if (!haveAlat) {
      printf("aLat has not been set!\n");
      exit(8);
    }

    if (!haveBlat) {
      printf("bLat has not been set!\n");
      exit(8);
    }

    if (!haveCells) {
      printf("The number of cells has not been set!\n");
      exit(8);
    }

    // Create lattice here:
    
    // Domains should be roughly square:
    
    state->nx = cells;
    state->ny = (int) ((double)cells * aLat / bLat );

    if (0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)) > bLat) {
       state->XYNNDIST =  0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0));
    } else {
       state->XYNNDIST = bLat;
    }
    // slightly larger so we can use <  as a comparison 
    state->XYNNDIST = state->XYNNDIST * 1.01;

    // each cell has 2 atoms (one at (0,0) and one at (a/2 , b/2))
    state->nAtoms = 2 * state->nx * state->ny;
    
    state->r =
      (struct coords *)calloc(state->nAtoms, sizeof(struct coords));

    which = 0;
    for(i=0; i < state->nx; i++) {            
      for(j=0; j < state->ny; j++) {     
        
        // First atom is at (0,0)
        
        (void)strcpy(state->r[which].name, atomName); 
        state->r[which].pos[0] = i * aLat;
        state->r[which].pos[1] = j * bLat;
        
        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        state->r[which].phi = myran * twopi;

        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        state->r[which].theta = acos(2.0*myran - 1.0);

        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        state->r[which].pos[2] = state->z0 + (2.0*myran-1.0) * state->deltaz;
        
        state->r[which].mu = state->strength;
        
        which++;
        // Second atom is at (a/2, b/2)

        state->r[which].pos[0] = aLat * (2.0 * i + 1.0) / 2.0;
        state->r[which].pos[1] = bLat * (2.0 * j + 1.0) / 2.0;

        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        state->r[which].phi = myran * twopi;

        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        state->r[which].theta = acos(2.0*myran - 1.0);

        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        state->r[which].pos[2] = state->z0 + (2.0*myran-1.0)*state->deltaz;

        state->r[which].mu = state->strength;
        which++;
      }
    }      
    state->Hmat[0][0] = state->nx * aLat;
    state->Hmat[0][1] = 0.0;
    state->Hmat[1][0] = 0.0;
    state->Hmat[1][1] = state->ny * bLat;  

    // Find HmatI:
    
    invertMat2(state->Hmat, state->HmatI);

  }  
  
  // Open the dump file for writing:
  
  if( !have_outName ) sprintf( out_name, "%s.dump", root_name ); 
  
  out_file = fopen( out_name, "w" );
  if( out_file == NULL ){
    printf("Cannot open file \"%s\" for writing.\n", out_name);
    exit(8);
  }
  
/*  if (state->iCycle >= state->nCycles) {
    printf("This configuration has already gone past the requested number\n"
           "of MC cycles!  Use the -n flag to request more MC cycles!\n");
    exit(8);
  }
*/
  // Do the MC simulation (finally!)
  
  en = toterg(state);

  printf("MC simulation starting for %d cycles\n", state->nCycles);
  printf("MC simulation starting with %d moves per cycle\n", state->nMoves);
  printf("MC simulation starting with sampling done every %d cycles\n", state->nSample);
  
  printf("The initial Energy is : \t%f\n\n", en);
  
  state->nAttempts = 0;
  state->nAccepts = 0;
  
  adjust(state);
  
  for(state->iCycle = 0; state->iCycle < state->nCycles; state->iCycle++) {
    for(imove=0; imove < state->nMoves; imove++) {
      mcmove(state, stream, &en);
    }
    
    if(((state->iCycle)%state->nSample) == 0) {
      store(state, out_file);
      /* Don't bother with the magnetization for now */
      /*
        getmag(state, mag);
        magmag = sqrt(pow(mag[0],2) + pow(mag[1],2) + pow(mag[2],2));
        printf("mag=%f\t%f\t%f\t%f\n", mag[0], mag[1], mag[2], magmag);
      */
    }
    
    if( ((state->iCycle+1)%(state->nCycles / 5)) ==0) {
      printf("\n=====> Completed\t%d\tout of\t%d cycles\n", 
             state->iCycle + 1, 
             state->nCycles);
      adjust(state);
    }
  }
  
  if(state->nCycles != 0) {
    if(state->nAttempts !=0) {
      printf("Number of attempts to modify a particle :%d\n"
             "Number of successful modifications: %d (=%lf%)\n\n", 
             state->nAttempts, state->nAccepts, 
             100*(double)(state->nAccepts)/(double)(state->nAttempts));
    }
    
    ent = toterg(state);
    
    if ( fabs(ent-en) > 1e-6) {
      printf("\n###### ENERGY PROBLEMS ###############\n\n");
      printf("Total Energy end of simulation : %lf\n"
             "Running Energy : %lf\n"
             "Energy Difference : %lf\n\n", ent, en, ent-en);
    }
  }
  fclose(out_file);
  vslDeleteStream( &stream );
  return 0;
}

double toterg(struct system* state) {
  
  struct coords iTemp; 
  
  int i, jb;
  double ener;
  double eni;

  ener = 0.0;

  for(i = 0; i < state->nAtoms; i++) {
    
    // Copy atom i's values into iTemp:
    
    iTemp.pos[0] = state->r[i].pos[0];
    iTemp.pos[1] = state->r[i].pos[1];
    iTemp.pos[2] = state->r[i].pos[2];
    iTemp.phi    = state->r[i].phi;
    iTemp.theta  = state->r[i].theta;
    iTemp.mu     = state->r[i].mu;
    
    // Pointer for eneri loop start is set to i
    // We'll skip the i->i pairing in the eneri routine
    
    jb = i;
    
    eni = eneri(state, iTemp, i, jb);
    ener += eni;
  }
  state->energy = ener;
  return ener;  
}

void getmag(struct system* state, double mag[3]) {
  
  double  thetai, phii, mui, magmag;
  int i;
  
  mag[0] = 0.0;
  mag[1] = 0.0;
  mag[2] = 0.0;
  
  for(i = 0; i < state->nAtoms; i++) {
    phii = state->r[i].phi;
    thetai = state->r[i].theta;
    mui = state->r[i].mu;
    mag[0] += mui * cos(phii) * sin(thetai);
    mag[1] += mui * sin(phii) * sin(thetai);
    mag[2] += mui * cos(thetai);
  }
  
  mag[0] /=  (double)(state->nAtoms);
  mag[1] /=  (double)(state->nAtoms);
  mag[2] /=  (double)(state->nAtoms);
  
}

double eneri(struct system *state, struct coords iTemp, int i, int jb) {
  
  double uxi, uyi, uzi, rxy, rij, r, r2, r3, r5, uxj, uyj, uzj, rcut;
  double udotu, rdotui, rdotuj, vij, pre, vint;
  double dx, dy, aLat, bLat;
  double eni, dz;
  double d[2];
  int j;

  vint = 0.0;
  
  pre = 14.38362;
  
  rcut = 30.0;

  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);
  
  eni = 0.0;
  uxi = iTemp.mu * cos(iTemp.phi) * sin(iTemp.theta);
  uyi = iTemp.mu * sin(iTemp.phi) * sin(iTemp.theta);
  uzi = iTemp.mu * cos(iTemp.theta);

  for(j = jb; j < state->nAtoms; j++) {
    if(j != i) {

      // 2-d wrapping on x and y:

      d[0] = state->r[j].pos[0] - iTemp.pos[0];
      d[1] = state->r[j].pos[1] - iTemp.pos[1];

      wrapVector(d, state->Hmat, state->HmatI);

      // z is unwrapped!

      dz = state->r[j].pos[2] - iTemp.pos[2];


      rxy = sqrt(pow(d[0], 2) + pow(d[1], 2));
      r2 = pow(d[0], 2) + pow(d[1], 2) + pow(dz, 2);
      r = sqrt(r2);

      if (rxy < state->XYNNDIST) {
         vint += 0.5 * state->kr * r2;
      }

      if(r < rcut) {

        r3 = r2*r;
        r5 = r2*r3;
        
        uxj = state->r[j].mu * cos(state->r[j].phi) * sin(state->r[j].theta);
        uyj = state->r[j].mu * sin(state->r[j].phi) * sin(state->r[j].theta);
        uzj = state->r[j].mu * cos(state->r[j].theta);
        udotu = uxi*uxj + uyi*uyj + uzi*uzj;
        rdotui = d[0]*uxi + d[1]*uyi + dz*uzi;
        rdotuj = d[0]*uxj + d[1]*uyj + dz*uzj;
        
        vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
        eni += vij;
      } 
    } 
  }
  
  vint += 0.5 * state->ktheta * pow((iTemp.theta - state->theta0), 2);
  eni += vint;
  return eni;
}

void adjust(struct system* state) {

  static int attempp, naccp;
  double dzo, dphio, dthetao, frac;
  
  if((state->nAttempts == 0) || (attempp >= state->nAttempts)) {
    naccp = state->nAccepts;
    attempp = state->nAttempts;
  } else {
    frac = (double)(state->nAccepts - naccp) /
      (double)(state->nAttempts - attempp);

    dthetao = state->dtheta;
    dzo = state->deltaz;
    dphio = state->deltaphi;
    
    state->dtheta   *= fabs(frac/0.5);
    state->deltaz   *= fabs(frac/0.5);
    state->deltaphi *= fabs(frac/0.5);
    
    if((state->dtheta/dthetao)>1.5) state->dtheta = dthetao*1.5;
    if((state->dtheta/dthetao)<0.5) state->dtheta = dthetao*0.5;

    if((state->deltaz/dzo)>1.5) state->deltaz = dzo * 1.5;
    if((state->deltaz/dzo)<0.5) state->deltaz = dzo * 0.5;
    if((state->deltaphi/dphio)>1.5) state->deltaphi = dphio * 1.5;
    if((state->deltaphi/dphio)<0.5) state->deltaphi = dphio * 0.5;
    printf("Max. displ. set to :\t%lf\t%lf\t%lf\n"
           "           (old was:\t%lf\t%lf\t%lf)\n"
           "Fractional acceptance:\t%lf\n"
           "nAttempts:\t%d\n"
           "nSuccesses:\t%d\n\n", 
           state->deltaz, state->deltaphi, state->dtheta, 
           dzo, dphio, dthetao, 
           frac, 
           state->nAttempts - attempp, 
           state->nAccepts - naccp);
    naccp = state->nAccepts;
    attempp=state->nAttempts;
  }  
}

void mcmove(struct system* state,  VSLStreamStatePtr stream, double *en) {
  
  int o, jb;
  struct coords oTemp;
  double eno, enn, myran;
  
  state->nAttempts++;
  
  jb = 0;
  
  vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
  
  o = (int)(state->nAtoms * myran);
  oTemp.pos[0] = state->r[o].pos[0];
  oTemp.pos[1] = state->r[o].pos[1];
  oTemp.pos[2] = state->r[o].pos[2];
  oTemp.phi    = state->r[o].phi;
  oTemp.theta  = state->r[o].theta;
  oTemp.mu     = state->r[o].mu;
  
  eno = eneri(state, oTemp, o, jb);
  
  vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
  
  oTemp.pos[2] += (2.0*myran-1.0) * state->deltaz;

  vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);

  oTemp.phi += (2.0*myran-1.0) * state->deltaphi;

  vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);

  oTemp.theta += (2.0*myran-1.0) * state->dtheta;

  enn = eneri(state, oTemp, o, jb);
  

  vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);

  if(myran <= (exp(- state->beta * (enn-eno)))) {
    state->nAccepts++; 
    *en += (enn-eno);
    state->r[o].pos[2] = oTemp.pos[2];
    state->r[o].phi    = oTemp.phi;
    state->r[o].theta  = oTemp.theta;
  }
}


void store(struct system* state, FILE* out_file) {

  double uxi, uyi, uzi;
  int i;
  
  fprintf(out_file,"%d\n",state->nAtoms);
 
  fprintf(out_file, "%d;\t%d\t%d;\t%f\t%f;\t%f\t%f;\n",
          state->iCycle,
          state->nx,
          state->ny,
          state->Hmat[0][0],
          state->Hmat[1][0],
          state->Hmat[0][1],
          state->Hmat[1][1]);
  
  for(i=0; i < state->nAtoms; i++) {
    uxi = cos(state->r[i].phi)*sin(state->r[i].theta);
    uyi = sin(state->r[i].phi)*sin(state->r[i].theta);
    uzi = cos(state->r[i].theta);
    fprintf(out_file, "%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",
            "Ar",
            state->r[i].pos[0],
            state->r[i].pos[1],
            state->r[i].pos[2],
            state->r[i].phi,
            uxi,
            uyi,
            uzi);
  }  

  fflush(out_file);
}

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

/***************************************************************************
 * prints out the usage for the command line arguments, then exits.
 ***************************************************************************/

void usage(){
  (void)fprintf(stderr,
                "The proper usage is: %s [options]\n"
                "\n"
                "Options:\n"
                "\n"
                "   -x              Display this message\n"
                "   -o <out_name>   The output file (Defaults to <root>.dump)\n"
                "   -i <in_name>    The input file (no default)\n"
                "   -r <root>       The root (Defaults to dp)\n"
                "   -n nCycles      The number of MC cycles to do\n"
                "   -s nSample      The number of MC cycles between samples\n"
                "   -m nMoves       The number of particle moves in each MC cycle\n"
                "   -a aLat         Set the lattice spacing along x axis\n"
                "   -b bLat         Set the lattice spacing along y axis\n"
                "   -c cells        Set the number of cells along x direction\n"
                "                   (Domains are kept nearly square)\n"
                "   -h hexSpace     Set up a hexagonal lattice\n"
                "                   aLat = sqrt(3) hexSpace\n"
                "                   bLat = hexSpace\n"
                "   -k kz           Sets the value of kz in units of kb\n"
                "   -t t           Sets the value of temperature in units of kelvin\n"
                "   -q strength    Sets the strength of the dipole\n"
                "\n",
                
                program_name);
  exit(8);
}
Revision:	1193
Committed:	Mon May 24 21:14:07 2004 UTC (20 years, 1 month ago) by xsun
Content type:	text/plain
File size:	28861 byte(s)
Log Message:	* empty log message *