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