1 |
xsun |
1185 |
#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 |
|
|
} |