1 |
#include <stdio.h> |
2 |
#include <stdlib.h> |
3 |
#include <string.h> |
4 |
#include <math.h> |
5 |
#include <unistd.h> |
6 |
#include <sys/types.h> |
7 |
#include <sys/stat.h> |
8 |
#include <sys/time.h> |
9 |
#include <time.h> |
10 |
#include "mkl_vsl.h" |
11 |
|
12 |
// inlined functions (for periodic box wrapping) |
13 |
|
14 |
inline double roundMe( double x ){ |
15 |
return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); |
16 |
} |
17 |
|
18 |
// Structures to store our data: |
19 |
|
20 |
// coords holds the data for a single tethered dipole: |
21 |
struct coords{ |
22 |
double pos[3]; // cartesian coords |
23 |
double theta; // orientational angle relative to z axis |
24 |
double phi; // orientational angle in x-y plane |
25 |
double mu; // dipole strength |
26 |
char name[30]; // an identifier for the type of atom |
27 |
}; |
28 |
|
29 |
// state holds the current "configuration" of the entire system |
30 |
struct system { |
31 |
int nAtoms; // Number of Atoms in this configuration |
32 |
struct coords *r; // The set of coordinates for all atoms |
33 |
double beta; // beta = 1 /(kb*T) |
34 |
double strength; // strength of the dipoles (Debye) |
35 |
double z0; // default z axis position |
36 |
double theta0; // default theta angle |
37 |
double kz; // force constant for z displacement |
38 |
double ktheta; // force constant for theta displacement |
39 |
int nCycles; // How many cycles to do in total |
40 |
int iCycle; // How many cycles have we done? |
41 |
int nMoves; // How many MC moves in each cycle |
42 |
int nSample; // How many cycles between samples |
43 |
double Hmat[2][2]; // The information about the size of the per. box |
44 |
double HmatI[2][2]; // The inverse box |
45 |
double energy; // The current Energy |
46 |
double dtheta; // maximum size of a theta move |
47 |
double deltaz; // maximum size of a z move |
48 |
double deltaphi; // maximum size of a phi move |
49 |
int nAttempts; // number of MC moves that have been attempted |
50 |
int nAccepts; // number of MC moves that have been accepted |
51 |
int nx; // number of unit cells in x direction |
52 |
int ny; // number of unit cells in y direction |
53 |
}; |
54 |
|
55 |
char *program_name; /* the name of the program */ |
56 |
|
57 |
// Function prototypes: |
58 |
void usage(void); |
59 |
double toterg(struct system* state); |
60 |
void getmag(struct system* state, double mag[3]); |
61 |
double eneri(struct system* state, struct coords iTemp, int i, int jb); |
62 |
void adjust(struct system* state); |
63 |
void mcmove(struct system* state, VSLStreamStatePtr stream, double *en); |
64 |
void store(struct system* state, FILE* out_file); |
65 |
void invertMat2(double a[2][2], double b[2][2]); |
66 |
void wrapVector( double thePos[2], double Hmat[2][2], double HmatI[2][2]); |
67 |
|
68 |
// Defines for the MKL random Number generator: |
69 |
#define SEED 1 |
70 |
#define BRNG VSL_BRNG_MCG31 |
71 |
#define METHOD 0 |
72 |
#define N 1 |
73 |
|
74 |
int main(argc, argv) |
75 |
int argc; |
76 |
char *argv[]; |
77 |
{ |
78 |
|
79 |
FILE *in_file; /* the input file */ |
80 |
FILE *out_file; /* the output file */ |
81 |
char *out_prefix = NULL; /* the prefix of the output file |
82 |
(if different from the root_name) */ |
83 |
char out_name[500]; /* the output name */ |
84 |
char in_name[500]; // the input file name |
85 |
char *root_name = NULL; /* the root name */ |
86 |
|
87 |
int have_outName = 0; |
88 |
int have_inName = 0; |
89 |
int have_rootName = 0; |
90 |
int restart_from_file = 0; |
91 |
|
92 |
int i, j, k, imove; |
93 |
|
94 |
char temp_name[500]; |
95 |
int lineCount = 0; // the line number |
96 |
char read_buffer[1000]; /* the line buffer for reading */ |
97 |
char *eof_test; /* ptr to see when we reach the end of the file */ |
98 |
char *foo; /* the pointer to the current string token */ |
99 |
char atomName[10]; |
100 |
|
101 |
double aLat; |
102 |
double bLat; |
103 |
int cells; |
104 |
|
105 |
int haveAlat = 0; |
106 |
int haveBlat = 0; |
107 |
int haveCells = 0; |
108 |
|
109 |
struct system* state; |
110 |
struct coords* r; |
111 |
|
112 |
int done; |
113 |
char current_flag; |
114 |
|
115 |
double dx, dy, uxi, uyi, uzi, myran, en, ent, magmag; |
116 |
double mag[3]; |
117 |
int which; |
118 |
|
119 |
// Other useful defines: |
120 |
|
121 |
double twopi = 2.0 * M_PI; |
122 |
double kb = 0.0019872198; |
123 |
strcpy(atomName, "Ar"); |
124 |
|
125 |
state = (struct system *)malloc(sizeof(struct system)); |
126 |
|
127 |
// The parameters for the simulation: |
128 |
|
129 |
state->strength = 7.0; |
130 |
state->z0 = 0.0; |
131 |
state->kz = kb; |
132 |
state->theta0 = M_PI / 2.0; |
133 |
state->ktheta = 0.0; |
134 |
state->dtheta = 0.3; |
135 |
state->deltaz = 1.0; |
136 |
state->deltaphi = 0.5; |
137 |
state->beta = 1.0 / (kb * 300.0); |
138 |
|
139 |
// these three should really be given as command line arguments, but we'll |
140 |
// set defaults here just in case: |
141 |
|
142 |
state->nCycles = 1000000; |
143 |
state->nMoves = 100; |
144 |
|
145 |
// Store roughly 100 frames while we are testing |
146 |
state->nSample = state->nCycles / 100; |
147 |
|
148 |
// Stuff for initializing the random number generator: |
149 |
struct timeval now_time_val; |
150 |
struct timezone time_zone; |
151 |
struct tm *now_tm; |
152 |
time_t now; |
153 |
int mySeed; |
154 |
VSLStreamStatePtr stream; |
155 |
|
156 |
// First, initialize the random number generator: |
157 |
gettimeofday(&now_time_val, &time_zone); /* get the time now */ |
158 |
now = now_time_val.tv_sec; /* convert to epoch time */ |
159 |
mySeed = (int) now; |
160 |
vslNewStream(&stream, BRNG, mySeed); |
161 |
|
162 |
// Now handle the arguments to the program: |
163 |
|
164 |
program_name = argv[0]; // save the program name in case we need it |
165 |
|
166 |
for( i = 1; i < argc; i++){ |
167 |
|
168 |
if(argv[i][0] =='-'){ |
169 |
|
170 |
// parse single character options |
171 |
|
172 |
done =0; |
173 |
j = 1; |
174 |
current_flag = argv[i][j]; |
175 |
while( (current_flag != '\0') && (!done) ){ |
176 |
|
177 |
switch(current_flag){ |
178 |
|
179 |
case 'o': |
180 |
// -o <filename> => the output |
181 |
|
182 |
i++; |
183 |
strcpy( out_name, argv[i] ); |
184 |
have_outName = 1; |
185 |
done = 1; |
186 |
break; |
187 |
|
188 |
case 'i': |
189 |
// -i <filename> => the input |
190 |
|
191 |
i++; |
192 |
strcpy( in_name, argv[i] ); |
193 |
have_inName = 1; |
194 |
done = 1; |
195 |
break; |
196 |
|
197 |
case 'r': |
198 |
// -r <root> => root Name |
199 |
|
200 |
i++; |
201 |
strcpy( root_name, argv[i] ); |
202 |
have_rootName = 1; |
203 |
done = 1; |
204 |
break; |
205 |
|
206 |
case 'n': |
207 |
// -n <#> => do Cycle MC cycles |
208 |
|
209 |
i++; |
210 |
state->nCycles = atoi(argv[i]); |
211 |
done = 1; |
212 |
break; |
213 |
|
214 |
case 's': |
215 |
// -s <#> => write out every nSample MC cycles |
216 |
printf("setting nSample\n"); |
217 |
i++; |
218 |
state->nSample = atoi(argv[i]); |
219 |
printf("setting nSample to %d\n", state->nSample); |
220 |
done = 1; |
221 |
break; |
222 |
|
223 |
case 'm': |
224 |
// -m <#> => each MC cycle consists of nMoves atomic moves |
225 |
|
226 |
i++; |
227 |
state->nMoves = atoi(argv[i]); |
228 |
done = 1; |
229 |
break; |
230 |
|
231 |
case 'k': |
232 |
// -k kz sets the value of kz in units of kb |
233 |
i++; |
234 |
state->kz = kb * atof( argv[i] ); |
235 |
done = 1; |
236 |
break; |
237 |
|
238 |
case 'h': |
239 |
// -h hexSpace set the inter-atomic spacing for a regular |
240 |
// hexagonal lattice |
241 |
haveAlat = 1; |
242 |
haveBlat = 1; |
243 |
i++; |
244 |
bLat = atof( argv[i] ); |
245 |
aLat = sqrt(3.0) * bLat; |
246 |
done = 1; |
247 |
break; |
248 |
|
249 |
case 'a': |
250 |
// -a aLat set the x lattice spacing for the distorted hexagonal |
251 |
// lattice |
252 |
haveAlat = 1; |
253 |
i++; |
254 |
aLat = atof( argv[i] ); |
255 |
done = 1; |
256 |
break; |
257 |
|
258 |
case 'b': |
259 |
// -b bLat set the y lattice spacing for the distorted hexagonal |
260 |
// lattice when initializing the system |
261 |
haveBlat = 1; |
262 |
i++; |
263 |
bLat = atof( argv[i] ); |
264 |
done = 1; |
265 |
break; |
266 |
|
267 |
case 'c': |
268 |
// -c cells sets the number of unit cells along the x axis |
269 |
// to use when initializing the system |
270 |
haveCells = 1; |
271 |
i++; |
272 |
cells = atoi( argv[i] ); |
273 |
done = 1; |
274 |
break; |
275 |
|
276 |
case 'x': |
277 |
usage(); |
278 |
done = 1; |
279 |
break; |
280 |
|
281 |
default: |
282 |
|
283 |
(void)fprintf(stderr, "Bad option \"-%s\"\n", current_flag); |
284 |
usage(); |
285 |
} |
286 |
j++; |
287 |
current_flag = argv[i][j]; |
288 |
} |
289 |
} else { |
290 |
|
291 |
if( root_name != NULL ){ |
292 |
fprintf( stderr, |
293 |
"The root name has already been set with an argument to -r\n" |
294 |
"But another argument looks like the root name. Whassup?\n"); |
295 |
usage(); |
296 |
} |
297 |
|
298 |
root_name = argv[i]; |
299 |
} |
300 |
} |
301 |
|
302 |
if ( !have_outName && root_name == NULL) { |
303 |
fprintf( stderr, "No output or root filename was set, so look for your data in dp.dump!\n"); |
304 |
strcpy(out_name, "dp.dump"); |
305 |
have_outName=1; |
306 |
} |
307 |
|
308 |
|
309 |
// Figure out if we are starting from an input file: |
310 |
|
311 |
if (have_inName) { |
312 |
|
313 |
restart_from_file = 1; |
314 |
|
315 |
in_file = fopen(in_name, "r"); |
316 |
if(in_file == NULL){ |
317 |
printf("Cannot open file \"%s\" for reading.\n", in_name); |
318 |
exit(8); |
319 |
} |
320 |
|
321 |
// start reading the first frame |
322 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
323 |
lineCount++; |
324 |
|
325 |
while(eof_test != NULL){ |
326 |
|
327 |
(void)sscanf(read_buffer, "%d", &state->nAtoms); |
328 |
state->r = |
329 |
(struct coords *)calloc(state->nAtoms, sizeof(struct coords)); |
330 |
|
331 |
// read and the comment line and grab the time and box dimensions |
332 |
|
333 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
334 |
lineCount++; |
335 |
if(eof_test == NULL){ |
336 |
printf("error in reading file at line: %d\n", lineCount); |
337 |
exit(8); |
338 |
} |
339 |
|
340 |
foo = strtok( read_buffer, " ,;\t" ); |
341 |
(void)sscanf( read_buffer, "%d", &state->iCycle ); |
342 |
|
343 |
foo = strtok(NULL, " ,;\t"); |
344 |
if(foo == NULL){ |
345 |
printf("error in reading file at line: %d\n", lineCount); |
346 |
exit(8); |
347 |
} |
348 |
(void)sscanf(foo, "%d", &state->nx); |
349 |
|
350 |
foo = strtok(NULL, " ,;\t"); |
351 |
if(foo == NULL){ |
352 |
printf("error in reading file at line: %d\n", lineCount); |
353 |
exit(8); |
354 |
} |
355 |
(void)sscanf(foo, "%d", &state->ny); |
356 |
|
357 |
foo = strtok(NULL, " ,;\t"); |
358 |
if(foo == NULL){ |
359 |
printf("error in reading file at line: %d\n", lineCount); |
360 |
exit(8); |
361 |
} |
362 |
(void)sscanf(foo, "%lf",&state->Hmat[0][0]); |
363 |
|
364 |
foo = strtok(NULL, " ,;\t"); |
365 |
if(foo == NULL){ |
366 |
printf("error in reading file at line: %d\n", lineCount); |
367 |
exit(8); |
368 |
} |
369 |
(void)sscanf(foo, "%lf",&state->Hmat[1][0]); |
370 |
|
371 |
|
372 |
foo = strtok(NULL, " ,;\t"); |
373 |
if(foo == NULL){ |
374 |
printf("error in reading file at line: %d\n", lineCount); |
375 |
exit(8); |
376 |
} |
377 |
(void)sscanf(foo, "%lf",&state->Hmat[0][1]); |
378 |
|
379 |
foo = strtok(NULL, " ,;\t"); |
380 |
if(foo == NULL){ |
381 |
printf("error in reading file at line: %d\n", lineCount); |
382 |
exit(8); |
383 |
} |
384 |
(void)sscanf(foo, "%lf",&state->Hmat[1][1]); |
385 |
|
386 |
// Length of the two box vectors: |
387 |
|
388 |
dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2)); |
389 |
dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2)); |
390 |
|
391 |
aLat = dx / (double)(state->nx); |
392 |
bLat = dy / (double)(state->ny); |
393 |
|
394 |
// Find HmatI: |
395 |
|
396 |
invertMat2(state->Hmat, state->HmatI); |
397 |
|
398 |
for( i=0; i < state->nAtoms; i++){ |
399 |
|
400 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
401 |
lineCount++; |
402 |
if(eof_test == NULL){ |
403 |
printf("error in reading file at line: %d\n", lineCount); |
404 |
exit(8); |
405 |
} |
406 |
|
407 |
foo = strtok(read_buffer, " ,;\t"); |
408 |
(void)strcpy(state->r[i].name, foo); //copy the atom name |
409 |
|
410 |
// next we grab the positions |
411 |
|
412 |
foo = strtok(NULL, " ,;\t"); |
413 |
if(foo == NULL){ |
414 |
printf("error in reading postition x from %s\n" |
415 |
"natoms = %d, line = %d\n", |
416 |
in_name, state->nAtoms, lineCount ); |
417 |
exit(8); |
418 |
} |
419 |
(void)sscanf( foo, "%lf", &state->r[i].pos[0] ); |
420 |
|
421 |
foo = strtok(NULL, " ,;\t"); |
422 |
if(foo == NULL){ |
423 |
printf("error in reading postition y from %s\n" |
424 |
"natoms = %d, line = %d\n", |
425 |
in_name, state->nAtoms, lineCount ); |
426 |
exit(8); |
427 |
} |
428 |
(void)sscanf( foo, "%lf", &state->r[i].pos[1] ); |
429 |
|
430 |
foo = strtok(NULL, " ,;\t"); |
431 |
if(foo == NULL){ |
432 |
printf("error in reading postition z from %s\n" |
433 |
"natoms = %d, line = %d\n", |
434 |
in_name, state->nAtoms, lineCount ); |
435 |
exit(8); |
436 |
} |
437 |
(void)sscanf( foo, "%lf", &state->r[i].pos[2] ); |
438 |
|
439 |
foo = strtok(NULL, " ,;\t"); |
440 |
if(foo == NULL){ |
441 |
printf("error in reading angle phi from %s\n" |
442 |
"natoms = %d, line = %d\n", |
443 |
in_name, state->nAtoms, lineCount ); |
444 |
exit(8); |
445 |
} |
446 |
(void)sscanf( foo, "%lf", &state->r[i].phi ); |
447 |
|
448 |
foo = strtok(NULL, " ,;\t"); |
449 |
if(foo == NULL){ |
450 |
printf("error in reading unit vector x from %s\n" |
451 |
"natoms = %d, line = %d\n", |
452 |
in_name, state->nAtoms, lineCount ); |
453 |
exit(8); |
454 |
} |
455 |
(void)sscanf( foo, "%lf", &uxi ); |
456 |
|
457 |
foo = strtok(NULL, " ,;\t"); |
458 |
if(foo == NULL){ |
459 |
printf("error in reading unit vector y from %s\n" |
460 |
"natoms = %d, line = %d\n", |
461 |
in_name, state->nAtoms, lineCount ); |
462 |
exit(8); |
463 |
} |
464 |
(void)sscanf( foo, "%lf", &uyi ); |
465 |
|
466 |
foo = strtok(NULL, " ,;\t"); |
467 |
if(foo == NULL){ |
468 |
printf("error in reading unit vector z from %s\n" |
469 |
"natoms = %d, line = %d\n", |
470 |
in_name, state->nAtoms, lineCount ); |
471 |
exit(8); |
472 |
} |
473 |
(void)sscanf( foo, "%lf", &uzi ); |
474 |
|
475 |
state->r[i].theta = acos(uzi); |
476 |
|
477 |
// The one parameter not stored in the dump file is the dipole strength |
478 |
state->r[i].mu = state->strength; |
479 |
|
480 |
} |
481 |
eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); |
482 |
lineCount++; |
483 |
} |
484 |
(void)fclose(in_file); |
485 |
|
486 |
} else { |
487 |
|
488 |
// not restarting from file, so use data we've got! |
489 |
|
490 |
if (!haveAlat) { |
491 |
printf("aLat has not been set!\n"); |
492 |
exit(8); |
493 |
} |
494 |
|
495 |
if (!haveBlat) { |
496 |
printf("bLat has not been set!\n"); |
497 |
exit(8); |
498 |
} |
499 |
|
500 |
if (!haveCells) { |
501 |
printf("The number of cells has not been set!\n"); |
502 |
exit(8); |
503 |
} |
504 |
|
505 |
// Create lattice here: |
506 |
|
507 |
// Domains should be roughly square: |
508 |
|
509 |
state->nx = cells; |
510 |
state->ny = (int) ((double)cells * aLat / bLat ); |
511 |
|
512 |
// each cell has 2 atoms (one at (0,0) and one at (a/2 , b/2)) |
513 |
state->nAtoms = 2 * state->nx * state->ny; |
514 |
|
515 |
state->r = |
516 |
(struct coords *)calloc(state->nAtoms, sizeof(struct coords)); |
517 |
|
518 |
which = 0; |
519 |
for(i=0; i < state->nx; i++) { |
520 |
for(j=0; j < state->ny; j++) { |
521 |
|
522 |
// First atom is at (0,0) |
523 |
|
524 |
(void)strcpy(state->r[which].name, atomName); |
525 |
state->r[which].pos[0] = i * aLat; |
526 |
state->r[which].pos[1] = j * bLat; |
527 |
|
528 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
529 |
state->r[which].phi = myran * twopi; |
530 |
|
531 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
532 |
state->r[which].theta = acos(2.0*myran - 1.0); |
533 |
|
534 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
535 |
state->r[which].pos[2] = state->z0 + (2.0*myran-1.0) * state->deltaz; |
536 |
|
537 |
state->r[which].mu = state->strength; |
538 |
|
539 |
which++; |
540 |
// Second atom is at (a/2, b/2) |
541 |
|
542 |
state->r[which].pos[0] = aLat * (2.0 * i + 1.0) / 2.0; |
543 |
state->r[which].pos[1] = bLat * (2.0 * j + 1.0) / 2.0; |
544 |
|
545 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
546 |
state->r[which].phi = myran * twopi; |
547 |
|
548 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
549 |
state->r[which].theta = acos(2.0*myran - 1.0); |
550 |
|
551 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
552 |
state->r[which].pos[2] = state->z0 + (2.0*myran-1.0)*state->deltaz; |
553 |
|
554 |
state->r[which].mu = state->strength; |
555 |
which++; |
556 |
} |
557 |
} |
558 |
state->Hmat[0][0] = state->nx * aLat; |
559 |
state->Hmat[0][1] = 0.0; |
560 |
state->Hmat[1][0] = 0.0; |
561 |
state->Hmat[1][1] = state->ny * bLat; |
562 |
|
563 |
// Find HmatI: |
564 |
|
565 |
invertMat2(state->Hmat, state->HmatI); |
566 |
|
567 |
} |
568 |
|
569 |
// Open the dump file for writing: |
570 |
|
571 |
if( !have_outName ) sprintf( out_name, "%s.dump", root_name ); |
572 |
|
573 |
out_file = fopen( out_name, "w" ); |
574 |
if( out_file == NULL ){ |
575 |
printf("Cannot open file \"%s\" for writing.\n", out_name); |
576 |
exit(8); |
577 |
} |
578 |
|
579 |
if (state->iCycle >= state->nCycles) { |
580 |
printf("This configuration has already gone past the requested number\n" |
581 |
"of MC cycles! Use the -n flag to request more MC cycles!\n"); |
582 |
exit(8); |
583 |
} |
584 |
|
585 |
// Do the MC simulation (finally!) |
586 |
|
587 |
en = toterg(state); |
588 |
|
589 |
printf("MC simulation starting for %d cycles\n", state->nCycles); |
590 |
printf("MC simulation starting with %d moves per cycle\n", state->nMoves); |
591 |
printf("MC simulation starting with sampling done every %d cycles\n", state->nSample); |
592 |
|
593 |
printf("The initial Energy is : \t%f\n\n", en); |
594 |
|
595 |
state->nAttempts = 0; |
596 |
state->nAccepts = 0; |
597 |
|
598 |
adjust(state); |
599 |
|
600 |
for(state->iCycle = 0; state->iCycle < state->nCycles; state->iCycle++) { |
601 |
for(imove=0; imove < state->nMoves; imove++) { |
602 |
mcmove(state, stream, &en); |
603 |
} |
604 |
|
605 |
if(((state->iCycle)%state->nSample) == 0) { |
606 |
store(state, out_file); |
607 |
/* Don't bother with the magnetization for now */ |
608 |
/* |
609 |
getmag(state, mag); |
610 |
magmag = sqrt(pow(mag[0],2) + pow(mag[1],2) + pow(mag[2],2)); |
611 |
printf("mag=%f\t%f\t%f\t%f\n", mag[0], mag[1], mag[2], magmag); |
612 |
*/ |
613 |
} |
614 |
|
615 |
if( ((state->iCycle+1)%(state->nCycles / 5)) ==0) { |
616 |
printf("\n=====> Completed\t%d\tout of\t%d cycles\n", |
617 |
state->iCycle + 1, |
618 |
state->nCycles); |
619 |
adjust(state); |
620 |
} |
621 |
} |
622 |
|
623 |
if(state->nCycles != 0) { |
624 |
if(state->nAttempts !=0) { |
625 |
printf("Number of attempts to modify a particle :%d\n" |
626 |
"Number of successful modifications: %d (=%lf%)\n\n", |
627 |
state->nAttempts, state->nAccepts, |
628 |
100*(double)(state->nAccepts)/(double)(state->nAttempts)); |
629 |
} |
630 |
|
631 |
ent = toterg(state); |
632 |
|
633 |
if ( fabs(ent-en) > 1e-6) { |
634 |
printf("\n###### ENERGY PROBLEMS ###############\n\n"); |
635 |
printf("Total Energy end of simulation : %lf\n" |
636 |
"Running Energy : %lf\n" |
637 |
"Energy Difference : %lf\n\n", ent, en, ent-en); |
638 |
} |
639 |
} |
640 |
fclose(out_file); |
641 |
vslDeleteStream( &stream ); |
642 |
return 0; |
643 |
} |
644 |
|
645 |
double toterg(struct system* state) { |
646 |
|
647 |
struct coords iTemp; |
648 |
|
649 |
int i, jb; |
650 |
double ener; |
651 |
double eni; |
652 |
|
653 |
ener = 0.0; |
654 |
|
655 |
for(i = 0; i < state->nAtoms; i++) { |
656 |
|
657 |
// Copy atom i's values into iTemp: |
658 |
|
659 |
iTemp.pos[0] = state->r[i].pos[0]; |
660 |
iTemp.pos[1] = state->r[i].pos[1]; |
661 |
iTemp.pos[2] = state->r[i].pos[2]; |
662 |
iTemp.phi = state->r[i].phi; |
663 |
iTemp.theta = state->r[i].theta; |
664 |
iTemp.mu = state->r[i].mu; |
665 |
|
666 |
// Pointer for eneri loop start is set to i |
667 |
// We'll skip the i->i pairing in the eneri routine |
668 |
|
669 |
jb = i; |
670 |
|
671 |
eni = eneri(state, iTemp, i, jb); |
672 |
ener += eni; |
673 |
} |
674 |
state->energy = ener; |
675 |
return ener; |
676 |
} |
677 |
|
678 |
void getmag(struct system* state, double mag[3]) { |
679 |
|
680 |
double thetai, phii, mui, magmag; |
681 |
int i; |
682 |
|
683 |
mag[0] = 0.0; |
684 |
mag[1] = 0.0; |
685 |
mag[2] = 0.0; |
686 |
|
687 |
for(i = 0; i < state->nAtoms; i++) { |
688 |
phii = state->r[i].phi; |
689 |
thetai = state->r[i].theta; |
690 |
mui = state->r[i].mu; |
691 |
mag[0] += mui * cos(phii) * sin(thetai); |
692 |
mag[1] += mui * sin(phii) * sin(thetai); |
693 |
mag[2] += mui * cos(thetai); |
694 |
} |
695 |
|
696 |
mag[0] /= (double)(state->nAtoms); |
697 |
mag[1] /= (double)(state->nAtoms); |
698 |
mag[2] /= (double)(state->nAtoms); |
699 |
|
700 |
} |
701 |
|
702 |
double eneri(struct system *state, struct coords iTemp, int i, int jb) { |
703 |
|
704 |
double uxi, uyi, uzi, r, r2, r3, r5, uxj, uyj, uzj, rcut; |
705 |
double udotu, rdotui, rdotuj, vij, pre, vint; |
706 |
double eni, dz; |
707 |
double d[2]; |
708 |
int j; |
709 |
|
710 |
pre = 14.38362; |
711 |
|
712 |
rcut = 30.0; |
713 |
|
714 |
eni = 0.0; |
715 |
uxi = iTemp.mu * cos(iTemp.phi) * sin(iTemp.theta); |
716 |
uyi = iTemp.mu * sin(iTemp.phi) * sin(iTemp.theta); |
717 |
uzi = iTemp.mu * cos(iTemp.theta); |
718 |
|
719 |
for(j = jb; j < state->nAtoms; j++) { |
720 |
if(j != i) { |
721 |
|
722 |
// 2-d wrapping on x and y: |
723 |
|
724 |
d[0] = state->r[j].pos[0] - iTemp.pos[0]; |
725 |
d[1] = state->r[j].pos[1] - iTemp.pos[1]; |
726 |
|
727 |
wrapVector(d, state->Hmat, state->HmatI); |
728 |
|
729 |
// z is unwrapped! |
730 |
|
731 |
dz = state->r[j].pos[2] - iTemp.pos[2]; |
732 |
|
733 |
r2 = pow(d[0], 2) + pow(d[1], 2) + pow(dz, 2); |
734 |
r = sqrt(r2); |
735 |
|
736 |
if(r < rcut) { |
737 |
|
738 |
r3 = r2*r; |
739 |
r5 = r2*r3; |
740 |
|
741 |
uxj = state->r[j].mu * cos(state->r[j].phi) * sin(state->r[j].theta); |
742 |
uyj = state->r[j].mu * sin(state->r[j].phi) * sin(state->r[j].theta); |
743 |
uzj = state->r[j].mu * cos(state->r[j].theta); |
744 |
udotu = uxi*uxj + uyi*uyj + uzi*uzj; |
745 |
rdotui = d[0]*uxi + d[1]*uyi + dz*uzi; |
746 |
rdotuj = d[0]*uxj + d[1]*uyj + dz*uzj; |
747 |
|
748 |
vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5); |
749 |
eni += vij; |
750 |
} |
751 |
} |
752 |
} |
753 |
|
754 |
vint = 0.5 * state->kz * pow((iTemp.pos[2] - state->z0), 2); |
755 |
vint += 0.5 * state->ktheta * pow((iTemp.theta - state->theta0), 2); |
756 |
eni += vint; |
757 |
return eni; |
758 |
} |
759 |
|
760 |
void adjust(struct system* state) { |
761 |
|
762 |
static int attempp, naccp; |
763 |
double dzo, dphio, dthetao, frac; |
764 |
|
765 |
if((state->nAttempts == 0) || (attempp >= state->nAttempts)) { |
766 |
naccp = state->nAccepts; |
767 |
attempp = state->nAttempts; |
768 |
} else { |
769 |
frac = (double)(state->nAccepts - naccp) / |
770 |
(double)(state->nAttempts - attempp); |
771 |
|
772 |
dthetao = state->dtheta; |
773 |
dzo = state->deltaz; |
774 |
dphio = state->deltaphi; |
775 |
|
776 |
state->dtheta *= fabs(frac/0.5); |
777 |
state->deltaz *= fabs(frac/0.5); |
778 |
state->deltaphi *= fabs(frac/0.5); |
779 |
|
780 |
if((state->dtheta/dthetao)>1.5) state->dtheta = dthetao*1.5; |
781 |
if((state->dtheta/dthetao)<0.5) state->dtheta = dthetao*0.5; |
782 |
|
783 |
if((state->deltaz/dzo)>1.5) state->deltaz = dzo * 1.5; |
784 |
if((state->deltaz/dzo)<0.5) state->deltaz = dzo * 0.5; |
785 |
if((state->deltaphi/dphio)>1.5) state->deltaphi = dphio * 1.5; |
786 |
if((state->deltaphi/dphio)<0.5) state->deltaphi = dphio * 0.5; |
787 |
printf("Max. displ. set to :\t%lf\t%lf\t%lf\n" |
788 |
" (old was:\t%lf\t%lf\t%lf)\n" |
789 |
"Fractional acceptance:\t%lf\n" |
790 |
"nAttempts:\t%d\n" |
791 |
"nSuccesses:\t%d\n\n", |
792 |
state->deltaz, state->deltaphi, state->dtheta, |
793 |
dzo, dphio, dthetao, |
794 |
frac, |
795 |
state->nAttempts - attempp, |
796 |
state->nAccepts - naccp); |
797 |
naccp = state->nAccepts; |
798 |
attempp=state->nAttempts; |
799 |
} |
800 |
} |
801 |
|
802 |
void mcmove(struct system* state, VSLStreamStatePtr stream, double *en) { |
803 |
|
804 |
int o, jb; |
805 |
struct coords oTemp; |
806 |
double eno, enn, myran; |
807 |
|
808 |
state->nAttempts++; |
809 |
|
810 |
jb = 0; |
811 |
|
812 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
813 |
|
814 |
o = (int)(state->nAtoms * myran); |
815 |
oTemp.pos[0] = state->r[o].pos[0]; |
816 |
oTemp.pos[1] = state->r[o].pos[1]; |
817 |
oTemp.pos[2] = state->r[o].pos[2]; |
818 |
oTemp.phi = state->r[o].phi; |
819 |
oTemp.theta = state->r[o].theta; |
820 |
oTemp.mu = state->r[o].mu; |
821 |
|
822 |
eno = eneri(state, oTemp, o, jb); |
823 |
|
824 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
825 |
|
826 |
oTemp.pos[2] += (2.0*myran-1.0) * state->deltaz; |
827 |
|
828 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
829 |
|
830 |
oTemp.phi += (2.0*myran-1.0) * state->deltaphi; |
831 |
|
832 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
833 |
|
834 |
oTemp.theta += (2.0*myran-1.0) * state->dtheta; |
835 |
|
836 |
enn = eneri(state, oTemp, o, jb); |
837 |
|
838 |
|
839 |
vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0); |
840 |
|
841 |
if(myran <= (exp(- state->beta * (enn-eno)))) { |
842 |
state->nAccepts++; |
843 |
*en += (enn-eno); |
844 |
state->r[o].pos[2] = oTemp.pos[2]; |
845 |
state->r[o].phi = oTemp.phi; |
846 |
state->r[o].theta = oTemp.theta; |
847 |
} |
848 |
} |
849 |
|
850 |
|
851 |
void store(struct system* state, FILE* out_file) { |
852 |
|
853 |
double uxi, uyi, uzi; |
854 |
int i; |
855 |
|
856 |
fprintf(out_file,"%d\n",state->nAtoms); |
857 |
|
858 |
fprintf(out_file, "%d;\t%d\t%d;\t%f\t%f;\t%f\t%f;\n", |
859 |
state->iCycle, |
860 |
state->nx, |
861 |
state->ny, |
862 |
state->Hmat[0][0], |
863 |
state->Hmat[1][0], |
864 |
state->Hmat[0][1], |
865 |
state->Hmat[1][1]); |
866 |
|
867 |
for(i=0; i < state->nAtoms; i++) { |
868 |
uxi = cos(state->r[i].phi)*sin(state->r[i].theta); |
869 |
uyi = sin(state->r[i].phi)*sin(state->r[i].theta); |
870 |
uzi = cos(state->r[i].theta); |
871 |
fprintf(out_file, "%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", |
872 |
"Ar", |
873 |
state->r[i].pos[0], |
874 |
state->r[i].pos[1], |
875 |
state->r[i].pos[2], |
876 |
state->r[i].phi, |
877 |
uxi, |
878 |
uyi, |
879 |
uzi); |
880 |
} |
881 |
|
882 |
fflush(out_file); |
883 |
} |
884 |
|
885 |
double matDet2(double a[2][2]) { |
886 |
|
887 |
double determinant; |
888 |
|
889 |
determinant = (a[0][0] * a[1][1]) - (a[0][1] * a[1][0]); |
890 |
|
891 |
return determinant; |
892 |
} |
893 |
|
894 |
|
895 |
void invertMat2(double a[2][2], double b[2][2]) { |
896 |
|
897 |
double determinant; |
898 |
|
899 |
determinant = matDet2( a ); |
900 |
|
901 |
if (determinant == 0.0) { |
902 |
printf("Can't invert a matrix with a zero determinant!\n"); |
903 |
} |
904 |
|
905 |
b[0][0] = a[1][1] / determinant; |
906 |
b[0][1] = -a[0][1] / determinant; |
907 |
b[1][0] = -a[1][0] / determinant; |
908 |
b[1][1] = a[0][0] / determinant; |
909 |
} |
910 |
|
911 |
void matVecMul2(double m[2][2], double inVec[2], double outVec[2]) { |
912 |
double a0, a1, a2; |
913 |
|
914 |
a0 = inVec[0]; a1 = inVec[1]; |
915 |
|
916 |
outVec[0] = m[0][0]*a0 + m[0][1]*a1; |
917 |
outVec[1] = m[1][0]*a0 + m[1][1]*a1; |
918 |
} |
919 |
|
920 |
void wrapVector( double thePos[2], double Hmat[2][2], double HmatInv[2][2]){ |
921 |
|
922 |
int i; |
923 |
double scaled[2]; |
924 |
|
925 |
// calc the scaled coordinates. |
926 |
|
927 |
matVecMul2(HmatInv, thePos, scaled); |
928 |
|
929 |
for(i=0; i<2; i++) |
930 |
scaled[i] -= roundMe(scaled[i]); |
931 |
|
932 |
// calc the wrapped real coordinates from the wrapped scaled coordinates |
933 |
|
934 |
matVecMul2(Hmat, scaled, thePos); |
935 |
|
936 |
} |
937 |
|
938 |
/*************************************************************************** |
939 |
* prints out the usage for the command line arguments, then exits. |
940 |
***************************************************************************/ |
941 |
|
942 |
void usage(){ |
943 |
(void)fprintf(stderr, |
944 |
"The proper usage is: %s [options]\n" |
945 |
"\n" |
946 |
"Options:\n" |
947 |
"\n" |
948 |
" -x Display this message\n" |
949 |
" -o <out_name> The output file (Defaults to <root>.dump)\n" |
950 |
" -i <in_name> The input file (no default)\n" |
951 |
" -r <root> The root (Defaults to dp)\n" |
952 |
" -n nCycles The number of MC cycles to do\n" |
953 |
" -s nSample The number of MC cycles between samples\n" |
954 |
" -m nMoves The number of particle moves in each MC cycle\n" |
955 |
" -a aLat Set the lattice spacing along x axis\n" |
956 |
" -b bLat Set the lattice spacing along y axis\n" |
957 |
" -c cells Set the number of cells along x direction\n" |
958 |
" (Domains are kept nearly square)\n" |
959 |
" -h hexSpace Set up a hexagonal lattice\n" |
960 |
" aLat = sqrt(3) hexSpace\n" |
961 |
" bLat = hexSpace\n" |
962 |
" -k kz Sets the value of kz in units of kb\n" |
963 |
"\n", |
964 |
|
965 |
program_name); |
966 |
exit(8); |
967 |
} |