Can we use Monte Carlo to compute Z?

1. Can we use Monte Carlo to compute Z? (even for a very simple, minimum-size, discrete case) • Consider a model of spins on a 2D lattice (idealized magnetic model). • ⇒ a discrete system • Each of N spins can take 2 states: up (↑) and down (↓). • ⇒ 2N microstates • Each spin interacts only with its nearest neighbors (nn). • ⇒ 4 neighbors for a 2D square lattice • Suppose that it takes 10-6 (or 10-15) s to compute the interaction of a spin with its neighbors. • ⇒ Time to calculate the energy Ei of a microstate i = N x 10-6 (or N x 10-15) s • For N = 100 spins: • - 2100 ~ 1030 microstates • - 10-3 (or 10-12) s to calculate the energy of a microstate • ⇒ 1027 (or 1018) s to estimate Z! (~ age of the universe ~ 13.8 billion years ~ 4 x 1017 s) • The situation gets worse with a real (continuous, larger, with beyond-nn interaction) system! a microstate ⇒ We cannot compute Z (and the absolute free energy kT ln Z) for real systems! (However, we can compute a relative free energy between two states.)

2. MC simulation (beyond-1-dimension-1-body): 2-Dimensional square-lattice “Ising model”for ferromagnetic phase transition • One of the simplest and best studied of statistical mechanical models • A simple model of a magnet • Proposed to interpret ferromagnetism in materials by Wilhelm Lenz (1920) • and analytically solved by Ernst Ising (1925; 1D) and Lars Onsager (1944; 2D) • Magnetic dipoles or atomic spinsσ(↑ or ↓, +1 or -1) are placed on a lattice. • For N sites on the lattice (N-atom model), the system can e in 2N microstates. • The total energy of a microstate is given by the Ising Hamiltonian: • J = interaction energy between nearest-neighbor spin <i, j> • H or B = external magnetic field (= 0 at the moment) •  = magnetic moment •  Let’s estimate quantities like magnetization m and specific heat c at a given temperature. a microstate

3. Specificheat = variance of energy Magneticsusceptibility = variance of magnetization

8. ising_2d.c #include "test.h" #include "build_ferro.h" #include "build_para.h" #include "init_glx_2d.h" #include "visu_glx_2d.h" #include "compute_pot.h" #include "mc_cycle.h" GLfloat xAngle = 90.0, yAngle = 0.0, zAngle = .0, xTrans = 0.95, yTrans = 0.0, zTrans = -24, xAngle_init, yAngle_init, zAngle_init, xPlane = 100; /*! main() for driver_test. */ int main(int argc, char *argv[]) { int i, j, n_dim, n_atoms, init, n_x, n_y, i_step, n_steps, n_blocks,n_data_block; double sigma, sigma2, a; double *h, **r, *s, **sl; int **nl; double temp; int *nn; long seed; double pot_energy, m, cv, pot_energy_block, pot_energy2_block, m_block, cv_block; int n_trial, n_accept; FILE *fp, *fpb;

9. ising_2d.c Display *dpy; Window win; GLboolean doubleBuffer = GL_TRUE; GLUquadricObj *qobj; int graph_switch, resph; int antialias_switch = 1; double radius; KeySym ks; char buff[31]; XEvent event; int keep_going; /* Get command-line input */ if (argc != 9) { fprintf(stderr, "Usage: %s init n_x n_y temp n_steps n_blocks seed graph_switch\n", argv[0]); exit(1); } init = atoi(argv[1]); n_x = atoi(argv[2]); n_y = atoi(argv[3]); temp = atof(argv[4]); n_steps = atoi(argv[5]); n_blocks = atoi(argv[6]); seed = atoi(argv[7]); graph_switch = atoi(argv[8]); n_dim = 2; sigma = 1.0;

10. ising_2d.c /* Build initial condition. */ if (init == 0) { /* build paramagnetic */ build_para(nn, a, r, s, nl, sl, &seed); } else { /* build ferromagnetic */ build_ferro(nn, a, r, s, nl, sl); } /* Initialize graphics. */ antialias_switch = 1; radius = sigma / 2.0; resph = 8; if (graph_switch > 0) { initialize_2d(&dpy, &win, doubleBuffer, &qobj, graph_switch,antialias_switch, h); keep_going = 1; while (keep_going) { XNextEvent(dpy, &event); switch (event.type) { case KeyPress: (void) XLookupString(&event.xkey, buff, sizeof(buff), &ks, NULL); switch (ks) { case XK_Escape: keep_going = 0; break; default: break; } } redraw_2d(dpy, win, doubleBuffer, qobj, n_atoms, r, s, h, radius, resph); } }

11. ising_2d.c /* Compute total energy. */ pot_energy = compute_pot(n_atoms, nn, nl, sl); /* Print energy. */ fprintf(stderr, "\n\n Initial\n"); fprintf(stderr, "\n\nPotential energy = %12.8f\n", pot_energy); /* Initialize block variables. */ m_block = 0.0; pot_energy_block = 0.0; pot_energy2_block = 0.0; /* Open files. */ fp = fopen("inst.dat", "w"); fpb = fopen("block.dat", "w"); /* MC loop. */ for (i_step = 1; i_step <= n_steps; ++i_step) { /* Perform a MC cycle. */ mc_cycle(nn, n_atoms, nl, sl, s, temp, &seed, &pot_energy, &n_trial, &n_accept); /* Calculate and write thermodynamic quantities in file every cycle. */ m = 0; for (i = 0; i < n_atoms; ++i) { m += s[i]; } m /= (double) n_atoms; m = fabs(m); fprintf(fp, "%d %12.8f %12.8f\n", i_step, pot_energy / n_atoms, m);

12. /* Accumulate blocks. */ m_block += m; pot_energy_block += pot_energy; pot_energy2_block += pot_energy * pot_energy; if (i_step % n_blocks == 0) { m_block /= (double) n_blocks; pot_energy_block /= (double) n_blocks; pot_energy2_block /= (double) n_blocks; cv_block = (pot_energy2_block - pot_energy_block * pot_energy_block) / (temp * temp * n_atoms); fprintf(fpb, "%d %12.8f %12.8f %12.8f\n", i_step, pot_energy_block / n_atoms, m_block, cv_block); m_block = 0.0; pot_energy_block = 0.0; pot_energy2_block = 0.0; } /* Draw system. */ if (graph_switch > 0 && i_step % 1000 == 0) redraw_2d(dpy, win, doubleBuffer, qobj, n_atoms, r, s, h, radius, resph); } /* Print energy. */ fprintf(stderr, "\n\n Final\n"); fprintf(stderr, "\n\nPotential energy = %12.8f\n", pot_energy); fprintf(stderr, "\n\nAccepted moves = %12.8f\n", (double)n_accept / (double)n_trial); /* Close files. */ fclose(fp); fclose(fpb); /* Exit */ return(0); } ising_2d.c

13. Build_ferro.c #include "build_ferro.h" void build_ferro(int *nn, double a, double **r, double *s, int **nl, double **sl) { int i, j, ii, n_atoms; n_atoms = nn[0] * nn[1]; /* Loop over atoms. */ ii = 0; for (j = 0; j < nn[1]; ++j) { for (i = 0; i < nn[0]; ++i) { /* Generate lattice and spin */ r[ii][0] = -nn[0] * 0.5 * a + i * a; r[ii][1] = -nn[1] * 0.5 * a + j * a; nl[ii][0] = i; nl[ii][1] = j; s[ii] = 1.0; sl[i][j] = s[ii]; ++ii; } } if (ii != n_atoms) { printf("problem in ferro building\n"); exit(1); } printf("\nConfiguration building completed\n"); fflush(NULL); /* Return zero if the initial condition is successfully constructed. */ return; }

14. Periodicboundary condition: Implementation (2d case) y Ly/2 1. Real coordinates /* (xi, yi) particle i coordinates */ if (xi > Lx/2) xi = xi – Lx; else if (xi < -Lx/2) xi = xi + Lx; if (yi > Ly/2) yi = yi – Ly; else if (yi < -Ly/2) yi = yi + Ly; i yi xi -Lx/2 0 Lx/2 xi-Lx x -Ly/2 2. Scaled (between [-0.5,0.5]) coordinates (better to handleanycellshape): orthorombiccell case #define NINT(x) ((x) < 0.0 ? (int) ((x) - 0.5) : (int) ((x) + 0.5)) sxi = xi / Lx; /* (sxi, syi) particle i scaled coordinates */ syi = yi / Ly; sxi = NINT(sxi); /* Apply PBC */ syi = NINT(syi); xi = sxi * Lx; /* (xi, yi) particle i folded real coordinates */ yi = syi * Ly;

15. build_para.c #include "test.h" #include "build_para.h" #include "ran3.h" void build_para(int *nn, double a, double **r, double *s, int **nl, double **sl, long *seed) { int i, j, ii, s_up, s_down, spin_up, spin_down, n_atoms, assign; double rr; n_atoms = nn[0] * nn[1]; spin_up = n_atoms / 2; spin_down = spin_up; /* Loop over atoms. */ ii = 0; s_up = 0; s_down = 0; for (j = 0; j < nn[1]; ++j) { for (i = 0; i < nn[0]; ++i) { /* Generate lattice and spin */ r[ii][0] = -nn[0] * 0.5 * a + i * a; r[ii][1] = -nn[1] * 0.5 * a + j * a; nl[ii][0] = i; nl[ii][1] = j;

16. build_para.c /* Assign spins */ do { assign = 0; rr = ran3(seed); if (rr < 0.5 && s_up < spin_up) { s[ii] = 1.0; s_up += 1; assign = 1; } else if (rr >= 0.5 && s_down < spin_down) { s[ii] = -1.0; s_down += 1; assign = 1; } } while (!assign); sl[i][j] = s[ii]; ++ii; } } if (ii != (spin_up + spin_down)) { printf("problem in para building\n"); exit(1); } printf("\nConfiguration building completed\n"); fflush(NULL); /* Return zero if the initial condition is successfully constructed. */ return; }

17. comput_pot.c pbc.c #include "test.h" #include "compute_pot.h" #include "pbc.h" double compute_pot(int n_atoms, int *nn, int **nl, double **sl) { int i, j, ip1, jp1, i_atom, n_x, n_y; double u, spin; n_x = nn[0]; n_y = nn[1]; /* Zero potential energy. */ u = 0.0; /* Loop over first atom. */ for (i_atom = 0; i_atom < n_atoms; ++i_atom) { i = nl[i_atom][0]; j = nl[i_atom][1]; spin = sl[i][j]; /* Compute 1/2 interaction on the square lattice. */ ip1 = pbc(i + 1, n_x); jp1 = pbc(j + 1, n_y); u -= spin * (sl[ip1][j] + sl[i][jp1]); } /* Return potential energy. */ return u; } #include<stdio.h> #include<stdlib.h> #include "pbc.h" int pbc(int i, int n) { if (i >= n) i -= n; else if (i < 0) i += n; return i; }

18. mc_cycle.c /* This routine carries out a single Monte Carlo cycle, comprising only of Monte Carlo atom move. Output: particle positions and potential energy are modified on output */ #include "test.h" #include "mc_cycle.h" #include "spin_flip.h" void mc_cycle(int *nn, int n_atoms, int **nl, double **sl, double *s, double temp, long *seed, double *u_pot, int *n_trial, int *n_accept) { int i_move; /* Loop over trial Monte Carlo moves (we only flip spins here). */ for (i_move = 0; i_move < n_atoms; ++i_move) { /* Make trial move. */ ++(*n_trial); *n_accept += spin_flip(nn, n_atoms, nl, sl, s, temp, seed, u_pot); } return; }

19. Spin_flip.c /* This routine generates a trial displacement of an interaction atom and accepts it based on the Metropolis criterion. Output: the return value is 1 if the move was accepted and 0 otherwise. */ #include "test.h" #include "pbc.h" #include "ran3.h" #include "spin_flip.h" int spin_flip(int *nn, int n_atoms, int **nl, double **sl, double *s, double temp, long *seed, double *u_pot) { int i, j, i_atom, accept; double spin, delta_u; int n_x, n_y, ip1, im1, jp1, jm1; /* Select an atom at random. */ i_atom = (int) (ran3(seed) * n_atoms); if (i_atom == n_atoms) --i_atom; n_x = nn[0]; n_y = nn[1]; i = nl[i_atom][0]; j = nl[i_atom][1]; spin = sl[i][j]; ip1 = pbc(i+1, n_x); im1 = pbc(i-1, n_x); jp1 = pbc(j+1, n_y); jm1 = pbc(j-1, n_y);

20. Spin_flip.c /* Compute change in potential energy and accept move using Metropolis criterion. */ delta_u = 2.0 * spin * (sl[im1][j] + sl[ip1][j] + sl[i][jm1] + sl[i][jp1]); if (delta_u <= 0.0) accept = 1; else accept = ran3(seed) < exp(- delta_u / temp); /* Update system potential energy if trial move is accepted. */ if (accept) { *u_pot += delta_u; sl[i][j] = -spin; s[i_atom] = -spin; } /* Return 1 if trial move is accepted, 0 otherwise. */ return accept; }