| 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 |
+ |
double t; // the temperature of the system |
| 41 |
|
int nCycles; // How many cycles to do in total |
| 42 |
|
int iCycle; // How many cycles have we done? |
| 43 |
|
int nMoves; // How many MC moves in each cycle |
| 52 |
|
int nAccepts; // number of MC moves that have been accepted |
| 53 |
|
int nx; // number of unit cells in x direction |
| 54 |
|
int ny; // number of unit cells in y direction |
| 55 |
+ |
double XYNNDIST; // maximum distance of nearest neighbors in XY plane |
| 56 |
|
}; |
| 57 |
|
|
| 58 |
|
char *program_name; /* the name of the program */ |
| 132 |
|
state->strength = 7.0; |
| 133 |
|
state->z0 = 0.0; |
| 134 |
|
state->kz = kb; |
| 135 |
+ |
state->kr = kb; |
| 136 |
|
state->theta0 = M_PI / 2.0; |
| 137 |
|
state->ktheta = 0.0; |
| 138 |
|
state->dtheta = 0.3; |
| 139 |
|
state->deltaz = 1.0; |
| 140 |
|
state->deltaphi = 0.5; |
| 141 |
< |
state->beta = 1.0 / (kb * 300.0); |
| 141 |
> |
state->t = 300.0; |
| 142 |
> |
state->beta = 1.0 / (kb * state->t); |
| 143 |
|
|
| 144 |
|
// these three should really be given as command line arguments, but we'll |
| 145 |
|
// set defaults here just in case: |
| 237 |
|
// -k kz sets the value of kz in units of kb |
| 238 |
|
i++; |
| 239 |
|
state->kz = kb * atof( argv[i] ); |
| 240 |
+ |
state->kr = state->kz; |
| 241 |
|
done = 1; |
| 242 |
|
break; |
| 243 |
|
|
| 244 |
+ |
case 'q': |
| 245 |
+ |
// -q mu set the strength of the dipole |
| 246 |
+ |
i++; |
| 247 |
+ |
state->strength = atof( argv[i] ); |
| 248 |
+ |
done = 1; |
| 249 |
+ |
break; |
| 250 |
+ |
|
| 251 |
+ |
case 't': |
| 252 |
+ |
// -t set the temperature of the system |
| 253 |
+ |
i++; |
| 254 |
+ |
state->t = atof( argv[i] ); |
| 255 |
+ |
state->beta = 1.0 / (kb * state->t); |
| 256 |
+ |
done = 1; |
| 257 |
+ |
break; |
| 258 |
+ |
|
| 259 |
|
case 'h': |
| 260 |
|
// -h hexSpace set the inter-atomic spacing for a regular |
| 261 |
|
// hexagonal lattice |
| 412 |
|
aLat = dx / (double)(state->nx); |
| 413 |
|
bLat = dy / (double)(state->ny); |
| 414 |
|
|
| 415 |
+ |
if (0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)) > bLat) { |
| 416 |
+ |
state->XYNNDIST = 0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)); |
| 417 |
+ |
} else { |
| 418 |
+ |
state->XYNNDIST = bLat; |
| 419 |
+ |
} |
| 420 |
+ |
// slightly larger so we can use < as a comparison |
| 421 |
+ |
state->XYNNDIST = state->XYNNDIST * 1.01; |
| 422 |
+ |
|
| 423 |
|
// Find HmatI: |
| 424 |
|
|
| 425 |
|
invertMat2(state->Hmat, state->HmatI); |
| 538 |
|
state->nx = cells; |
| 539 |
|
state->ny = (int) ((double)cells * aLat / bLat ); |
| 540 |
|
|
| 541 |
+ |
if (0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)) > bLat) { |
| 542 |
+ |
state->XYNNDIST = 0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)); |
| 543 |
+ |
} else { |
| 544 |
+ |
state->XYNNDIST = bLat; |
| 545 |
+ |
} |
| 546 |
+ |
// slightly larger so we can use < as a comparison |
| 547 |
+ |
state->XYNNDIST = state->XYNNDIST * 1.01; |
| 548 |
+ |
|
| 549 |
|
// each cell has 2 atoms (one at (0,0) and one at (a/2 , b/2)) |
| 550 |
|
state->nAtoms = 2 * state->nx * state->ny; |
| 551 |
|
|
| 613 |
|
exit(8); |
| 614 |
|
} |
| 615 |
|
|
| 616 |
< |
if (state->iCycle >= state->nCycles) { |
| 616 |
> |
/* if (state->iCycle >= state->nCycles) { |
| 617 |
|
printf("This configuration has already gone past the requested number\n" |
| 618 |
|
"of MC cycles! Use the -n flag to request more MC cycles!\n"); |
| 619 |
|
exit(8); |
| 620 |
|
} |
| 621 |
< |
|
| 621 |
> |
*/ |
| 622 |
|
// Do the MC simulation (finally!) |
| 623 |
|
|
| 624 |
|
en = toterg(state); |
| 738 |
|
|
| 739 |
|
double eneri(struct system *state, struct coords iTemp, int i, int jb) { |
| 740 |
|
|
| 741 |
< |
double uxi, uyi, uzi, r, r2, r3, r5, uxj, uyj, uzj, rcut; |
| 741 |
> |
double uxi, uyi, uzi, rxy, rij, r, r2, r3, r5, uxj, uyj, uzj, rcut; |
| 742 |
|
double udotu, rdotui, rdotuj, vij, pre, vint; |
| 743 |
+ |
double dx, dy, aLat, bLat; |
| 744 |
|
double eni, dz; |
| 745 |
|
double d[2]; |
| 746 |
|
int j; |
| 747 |
+ |
|
| 748 |
+ |
vint = 0.0; |
| 749 |
|
|
| 750 |
|
pre = 14.38362; |
| 751 |
|
|
| 752 |
|
rcut = 30.0; |
| 753 |
+ |
|
| 754 |
+ |
dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2)); |
| 755 |
+ |
dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2)); |
| 756 |
+ |
|
| 757 |
+ |
aLat = dx / (double)(state->nx); |
| 758 |
+ |
bLat = dy / (double)(state->ny); |
| 759 |
|
|
| 760 |
|
eni = 0.0; |
| 761 |
|
uxi = iTemp.mu * cos(iTemp.phi) * sin(iTemp.theta); |
| 775 |
|
// z is unwrapped! |
| 776 |
|
|
| 777 |
|
dz = state->r[j].pos[2] - iTemp.pos[2]; |
| 778 |
< |
|
| 778 |
> |
|
| 779 |
> |
|
| 780 |
> |
rxy = sqrt(pow(d[0], 2) + pow(d[1], 2)); |
| 781 |
|
r2 = pow(d[0], 2) + pow(d[1], 2) + pow(dz, 2); |
| 782 |
|
r = sqrt(r2); |
| 783 |
|
|
| 784 |
+ |
if (rxy < state->XYNNDIST) { |
| 785 |
+ |
vint += 0.5 * state->kr * r2; |
| 786 |
+ |
} |
| 787 |
+ |
|
| 788 |
|
if(r < rcut) { |
| 789 |
|
|
| 790 |
|
r3 = r2*r; |
| 803 |
|
} |
| 804 |
|
} |
| 805 |
|
|
| 754 |
– |
vint = 0.5 * state->kz * pow((iTemp.pos[2] - state->z0), 2); |
| 806 |
|
vint += 0.5 * state->ktheta * pow((iTemp.theta - state->theta0), 2); |
| 807 |
|
eni += vint; |
| 808 |
|
return eni; |
| 1011 |
|
" aLat = sqrt(3) hexSpace\n" |
| 1012 |
|
" bLat = hexSpace\n" |
| 1013 |
|
" -k kz Sets the value of kz in units of kb\n" |
| 1014 |
+ |
" -t t Sets the value of temperature in units of kelvin\n" |
| 1015 |
+ |
" -q strength Sets the strength of the dipole\n" |
| 1016 |
|
"\n", |
| 1017 |
|
|
| 1018 |
|
program_name); |
