1 / 20

Can we use Monte Carlo to compute Z?

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 (↓). ⇒ 2 N microstates

halen
Download Presentation

Can we use Monte Carlo to compute Z?

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  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

  4. 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;

  5. 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;

  6. 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); } }

  7. 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);

  8. /* 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

  9. 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; }

  10. 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;

  11. 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;

  12. 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; }

  13. 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; }

  14. 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; }

  15. 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);

  16. 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; }

More Related