trunk/ripple/readh.c

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

// 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 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:

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;
  double uxi, uyi, uzi, p1;
  double aLat, bLat;
  int cells;
  int i, j, 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 whichbinx, whichbiny, minwhichbinx, minwhichbiny, maxwhichbinx, maxwhichbiny, n;
  double hist[500][500], averdh2[500][500];
  int nbins, nx, ny;
  int *samples;
  double *bin;
  double dh2;

  int done;
  char current_flag;
  
  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;

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

  nframes = 0;
  n = 0;
   
  // start reading the first frame

  eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
  lineCount++;

  state = (struct system *) malloc(sizeof(struct system));
  state->next = NULL;
  state->strength = 10.0;

  dh2 = 0.0;
  for(i = 0; i < 500; i++){
    for(j = 0; j < 500; j++){
      averdh2[i][j] = 0.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));
        
    // 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]);

    
    // 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", &uxi );
      
      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", &uyi );
      
      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", &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;
    }

    minwhichbinx = 0;
    minwhichbiny = 0;
    maxwhichbinx = 0;
    maxwhichbiny = 0;

    for(i = 0; i< (state->nAtoms - 1); i++){
      for(j = i + 1; j < state->nAtoms; j++){
        
        dx = state->r[j].pos[0] - state->r[i].pos[0];
        dy = state->r[j].pos[1] - state->r[i].pos[1];

        whichbinx = (int) (dx * 2 / aLat) + 2 * state->nx;
        whichbiny = (int) (dy * 2 / bLat) + 2 * state->ny;

        maxwhichbinx = ((maxwhichbinx > whichbinx) ? maxwhichbinx : whichbinx);
        maxwhichbiny = ((maxwhichbiny > whichbiny) ? maxwhichbiny : whichbiny);

        minwhichbinx = ((minwhichbinx < whichbinx) ? minwhichbinx : whichbinx);
        minwhichbiny = ((minwhichbiny < whichbiny) ? minwhichbiny : whichbiny);

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

        hist[whichbinx][whichbiny] += dh2;
      }
      n++;
    }

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

  for(whichbinx = minwhichbinx; whichbinx < (maxwhichbinx + 1); whichbinx++){
    for(whichbiny = minwhichbiny; whichbiny < (maxwhichbiny + 1); whichbiny++){
      averdh2[whichbinx][whichbiny] = hist[whichbinx][whichbiny] / (double)n;
    }
  }

  for(whichbinx = minwhichbinx; whichbinx < (maxwhichbinx + 1); whichbinx++){
    for(whichbiny = minwhichbiny; whichbiny < (maxwhichbiny + 1); whichbiny++){
      printf("%lf\t", averdh2[whichbinx][whichbiny]);
    }
    printf("\n");
  }

  return 1;

}
Revision:	1185
Committed:	Fri May 21 20:07:10 2004 UTC (21 years, 5 months ago) by xsun
Content type:	text/plain
File size:	9909 byte(s)
Log Message:	this is the ripple codes
#	User	Rev	Content
1	xsun	1185	#include<stdio.h>
2			#include<string.h>
3			#include<stdlib.h>
4			#include<math.h>
5
6			// Structures to store our data:
7
8			// coords holds the data for a single tethered dipole:
9			struct coords{
10			double pos[3]; // cartesian coords
11			double theta; // orientational angle relative to z axis
12			double phi; // orientational angle in x-y plane
13			double mu; // dipole strength
14			char name[30]; // an identifier for the type of atom
15			};
16
17			// state holds the current "configuration" of the entire system
18			struct system {
19			int nAtoms; // Number of Atoms in this configuration
20			struct coords *r; // The set of coordinates for all atoms
21			double beta; // beta = 1 /(kb*T)
22			double strength; // strength of the dipoles (Debye)
23			double z0; // default z axis position
24			double theta0; // default theta angle
25			double kz; // force constant for z displacement
26			double ktheta; // force constant for theta displacement
27			int nCycles; // How many cycles to do in total
28			int iCycle; // How many cycles have we done?
29			int nMoves; // How many MC moves in each cycle
30			int nSample; // How many cycles between samples
31			double Hmat[2][2]; // The information about the size of the per. box
32			double HmatI[2][2]; // The inverse box
33			double energy; // The current Energy
34			double dtheta; // maximum size of a theta move
35			double deltaz; // maximum size of a z move
36			double deltaphi; // maximum size of a phi move
37			int nAttempts; // number of MC moves that have been attempted
38			int nAccepts; // number of MC moves that have been accepted
39			int nx; // number of unit cells in x direction
40			int ny; // number of unit cells in y direction
41			struct system *next; // Next frame in the linked list
42			};
43
44			char program_name; / the name of the program */
45
46			// Function prototypes:
47
48			int main(argc, argv)
49			int argc;
50			char *argv[];
51			{
52			FILE *in_file;
53			char in_name[500];
54			char eof_test, foo;
55			char read_buffer[1000];
56			int lineCount = 0;
57			double uxi, uyi, uzi, p1;
58			double aLat, bLat;
59			int cells;
60			int i, j, nloops;
61			int newx, newy, newindex, index;
62			int new_i, new_j, new_index;
63			int N, nframes;
64			double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
65			double maxfreqx, maxfreqy, maxfreq, dfreq;
66			double dx, dy, dx1, dy1, xTemp, yTemp, pt1x, pt1y, pt2x, pt2y;
67			int whichbinx, whichbiny, minwhichbinx, minwhichbiny, maxwhichbinx, maxwhichbiny, n;
68			double hist[500][500], averdh2[500][500];
69			int nbins, nx, ny;
70			int *samples;
71			double *bin;
72			double dh2;
73
74			int done;
75			char current_flag;
76
77			program_name = argv[0];
78			if (argc >= 2)
79			strcpy(in_name, argv[1]);
80			/*
81			for(i = 1; i < argc; i++){
82			if(argv[i][0] == '-'){
83			done = 0;
84			j = 1;
85			current_flag = argv[i][j];
86			while( (current_flag != '\0') && ( !done ) ){
87			switch(current_flag){
88			case 'i':
89			i++;
90			strcpy( in_name, argv[i] );
91			done = 1;
92			break;
93			}
94			}
95			}
96			}
97			*/
98			struct system* state;
99			struct system* temp_state;
100			struct coords* r;
101
102			in_file = fopen(in_name, "r");
103			if(in_file == NULL){
104			printf("Cannot open file \"%s\" for reading.\n", in_name);
105			exit(8);
106			}
107
108			nframes = 0;
109			n = 0;
110
111			// start reading the first frame
112
113			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
114			lineCount++;
115
116			state = (struct system *) malloc(sizeof(struct system));
117			state->next = NULL;
118			state->strength = 10.0;
119
120			dh2 = 0.0;
121			for(i = 0; i < 500; i++){
122			for(j = 0; j < 500; j++){
123			averdh2[i][j] = 0.0;
124			}
125			}
126
127			while(eof_test != NULL){
128
129			nframes++;
130			(void)sscanf(read_buffer, "%d", &state->nAtoms);
131			N = 2 * state->nAtoms;
132
133			state->r = (struct coords *)calloc(N, sizeof(struct coords));
134
135			// read and the comment line and grab the time and box dimensions
136
137			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
138			lineCount++;
139			if(eof_test == NULL){
140			printf("error in reading file at line: %d\n", lineCount);
141			exit(8);
142			}
143
144			foo = strtok( read_buffer, " ,;\t\n" );
145			(void)sscanf( read_buffer, "%d", &state->iCycle );
146
147			foo = strtok(NULL, " ,;\t\0");
148			if(foo == NULL){
149			printf("error in reading file at line: %d\n", lineCount);
150			exit(8);
151			}
152			(void)sscanf(foo, "%d", &state->nx);
153
154			nx = state->nx;
155
156			foo = strtok(NULL, " ,;\t\0");
157			if(foo == NULL){
158			printf("error in reading file at line: %d\n", lineCount);
159			exit(8);
160			}
161			(void)sscanf(foo, "%d", &state->ny);
162
163			ny = state->ny;
164
165			foo = strtok(NULL, " ,;\t\0");
166			if(foo == NULL){
167			printf("error in reading file at line: %d\n", lineCount);
168			exit(8);
169			}
170			(void)sscanf(foo, "%lf",&state->Hmat[0][0]);
171
172			foo = strtok(NULL, " ,;\t\0");
173			if(foo == NULL){
174			printf("error in reading file at line: %d\n", lineCount);
175			exit(8);
176			}
177			(void)sscanf(foo, "%lf",&state->Hmat[1][0]);
178
179			foo = strtok(NULL, " ,;\t\0");
180			if(foo == NULL){
181			printf("error in reading file at line: %d\n", lineCount);
182			exit(8);
183			}
184			(void)sscanf(foo, "%lf",&state->Hmat[0][1]);
185
186			foo = strtok(NULL, " ,;\t\0");
187			if(foo == NULL){
188			printf("error in reading file at line: %d\n", lineCount);
189			exit(8);
190			}
191			(void)sscanf(foo, "%lf",&state->Hmat[1][1]);
192
193
194			// Length of the two box vectors:
195
196			dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
197			dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));
198
199			aLat = dx / (double)(state->nx);
200			bLat = dy / (double)(state->ny);
201
202
203			// FFT stuff depends on nx and ny, so delay allocation until we have
204			// that information
205
206			for (i=0;i<state->nAtoms;i++){
207
208			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
209			lineCount++;
210			if(eof_test == NULL){
211			printf("error in reading file at line: %d\n", lineCount);
212			exit(8);
213			}
214
215			foo = strtok(read_buffer, " ,;\t\0");
216			(void)strcpy(state->r[i].name, foo); //copy the atom name
217
218			// next we grab the positions
219
220			foo = strtok(NULL, " ,;\t\0");
221			if(foo == NULL){
222			printf("error in reading postition x from %s\n"
223			"natoms = %d, line = %d\n",
224			in_name, state->nAtoms, lineCount );
225			exit(8);
226			}
227			(void)sscanf( foo, "%lf", &state->r[i].pos[0] );
228
229			foo = strtok(NULL, " ,;\t\0");
230			if(foo == NULL){
231			printf("error in reading postition y from %s\n"
232			"natoms = %d, line = %d\n",
233			in_name, state->nAtoms, lineCount );
234			exit(8);
235			}
236			(void)sscanf( foo, "%lf", &state->r[i].pos[1] );
237
238			foo = strtok(NULL, " ,;\t\0");
239			if(foo == NULL){
240			printf("error in reading postition z from %s\n"
241			"natoms = %d, line = %d\n",
242			in_name, state->nAtoms, lineCount );
243			exit(8);
244			}
245			(void)sscanf( foo, "%lf", &state->r[i].pos[2] );
246
247			foo = strtok(NULL, " ,;\t\0");
248			if(foo == NULL){
249			printf("error in reading angle phi from %s\n"
250			"natoms = %d, line = %d\n",
251			in_name, state->nAtoms, lineCount );
252			exit(8);
253			}
254			(void)sscanf( foo, "%lf", &state->r[i].phi );
255
256			foo = strtok(NULL, " ,;\t\0");
257			if(foo == NULL){
258			printf("error in reading unit vector x from %s\n"
259			"natoms = %d, line = %d\n",
260			in_name, state->nAtoms, lineCount );
261			exit(8);
262			}
263			(void)sscanf( foo, "%lf", &uxi );
264
265			foo = strtok(NULL, " ,;\t\0");
266			if(foo == NULL){
267			printf("error in reading unit vector y from %s\n"
268			"natoms = %d, line = %d\n",
269			in_name, state->nAtoms, lineCount );
270			exit(8);
271			}
272			(void)sscanf( foo, "%lf", &uyi );
273
274			foo = strtok(NULL, " ,;\t\0");
275			if(foo == NULL){
276			printf("error in reading unit vector z from %s\n"
277			"natoms = %d, line = %d\n",
278			in_name, state->nAtoms, lineCount );
279			exit(8);
280			}
281			(void)sscanf( foo, "%lf", &uzi );
282
283			state->r[i].theta = acos(uzi);
284
285			// The one parameter not stored in the dump file is the dipole strength
286			state->r[i].mu = state->strength;
287			}
288
289			minwhichbinx = 0;
290			minwhichbiny = 0;
291			maxwhichbinx = 0;
292			maxwhichbiny = 0;
293
294			for(i = 0; i< (state->nAtoms - 1); i++){
295			for(j = i + 1; j < state->nAtoms; j++){
296
297			dx = state->r[j].pos[0] - state->r[i].pos[0];
298			dy = state->r[j].pos[1] - state->r[i].pos[1];
299
300			whichbinx = (int) (dx * 2 / aLat) + 2 * state->nx;
301			whichbiny = (int) (dy * 2 / bLat) + 2 * state->ny;
302
303			maxwhichbinx = ((maxwhichbinx > whichbinx) ? maxwhichbinx : whichbinx);
304			maxwhichbiny = ((maxwhichbiny > whichbiny) ? maxwhichbiny : whichbiny);
305
306			minwhichbinx = ((minwhichbinx < whichbinx) ? minwhichbinx : whichbinx);
307			minwhichbiny = ((minwhichbiny < whichbiny) ? minwhichbiny : whichbiny);
308
309			dh2 = pow((state->r[j].pos[2] - state->r[i].pos[2]), 2);
310
311			hist[whichbinx][whichbiny] += dh2;
312			}
313			n++;
314			}
315
316
317			temp_state = state->next;
318			state->next = NULL;
319
320			if (temp_state != NULL) {
321			free(temp_state->r);
322			free(temp_state);
323			}
324
325			// Make a new frame
326
327			temp_state = (struct system *) malloc(sizeof(struct system));
328			temp_state->next = state;
329			state = temp_state;
330			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
331			lineCount++;
332			}
333
334			for(whichbinx = minwhichbinx; whichbinx < (maxwhichbinx + 1); whichbinx++){
335			for(whichbiny = minwhichbiny; whichbiny < (maxwhichbiny + 1); whichbiny++){
336			averdh2[whichbinx][whichbiny] = hist[whichbinx][whichbiny] / (double)n;
337			}
338			}
339
340			for(whichbinx = minwhichbinx; whichbinx < (maxwhichbinx + 1); whichbinx++){
341			for(whichbiny = minwhichbiny; whichbiny < (maxwhichbiny + 1); whichbiny++){
342			printf("%lf\t", averdh2[whichbinx][whichbiny]);
343			}
344			printf("\n");
345			}
346
347			return 1;
348
349			}