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 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 * 950.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 |
|
|
} |