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 (20 years, 1 month ago) by xsun
Content type:	text/plain
File size:	9909 byte(s)
Log Message:	this is the ripple codes
#	Content
1	#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	}