ising/ising1d.cpp at main · UVaCompPhys/ising · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
/*
******************************************************************************
* ising1d.c                                                                  *
* =========                                                                  *
*    Uses Metropolis algorithm to generate thermal ensemble for the          *
*    one-dimensional ising chain with free boundary conditions               *
*    in a magnetic field h.                                                  *
*								             *
*    The state of the lattice is stored as spins +-1 in elements 1     *
*    through N of an array of length N+2.  The free boundary conditions      *
*    are handled by fixing elements 0 and N+1 of the array to be zero.       *
*  								             *
******************************************************************************
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>


/* update_spin (s, env)
** -----------
** Do a metropolis update on a spin *s whose "environment" is env.
** The environment should be set beforehand so that the dependence
** of the total beta*energy on the value of the selected spin is
** given by (selected spin)*env
*/


void update_spin(int *s, double env) {
  int spin = *s;
  int newspin = ( drand48() < 0.5 ? 1 : -1 );
  double DeltaBetaE = -(newspin-spin)*env;       /* beta*(E(new) - E(old)) */
  if ( DeltaBetaE <= 0 || drand48() < exp(-DeltaBetaE) ) *s = newspin;
}

/* sweep (spin, N, beta, h)
** -----
** Sweep once through all the sites of the lattice spin of length N,
** trying an update at each site with inverse temperature beta and external
** magnetic field parameter h.
*/
void sweep(int *spin, int N, double beta, double h) {
  for(int ns=1; ns<=N; ns++) {
    update_spin(&(spin[ns]), beta*(spin[ns-1] + spin[ns+1]) + h);
  }
}

/* InitializeAllSpinUp (spin, N)
** =============================
**   Initialize all the N spins of the lattice "spin" to be +1.  Also
** initialize the two fake "boundary" spins to zero to implement free
** boundary conditions.
*/

void InitializeAllSpinUp(int *spin, int N) {
  for(int ns=1; ns<=N; ns++) spin[ns] = 1;
  spin[0] = spin[N+1] = 0;             // no interaction w/ endpoints
}

void HotStart(int *spin, int N) {
  for(int ns=1; ns<=N; ns++)
    drand48()>0.5 ? spin[ns] = 1 : spin[ns] = -1;
  spin[0] = spin[N+1] = 0;             // no interaction w/ endpoints
}


int main(){
  int ntherm=10000;           // initial sweeps to "thermalize" the system
  const int VisualDisplay=1;  // 1 or 0 to turn on/off display of chains

  int nsweep, N, nmag=0, ntotal=0;
  int *spin;
  double beta, h;
  double expected;


  printf("Program generates a thermal ensemble for chain of N spins\n");
  printf("with free boundary conditions.\n\n");

  srand48((long)time(0));               // seed the number generator

  printf("Enter number of spins in chain:\n");
  scanf("%d", &N);

  printf("Enter total number of configurations (sweeps) generated:\n");
  scanf("%d", &nsweep);

  printf("Enter value of magnetic field parameter h:\n");
  scanf("%lf", &h);

  printf("Enter temperature parameter beta (= 1/kT):\n");
  scanf("%lf", &beta);

  // make room for N spins and two fake boundary spins
  spin = (int*)malloc((N+2)*sizeof(int));
  //InitializeAllSpinUp(spin, N);  // cold start
  HotStart(spin,N);   // random spins to start

  // sweep ntherm times to thermalize system
  for(int i=0; i<ntherm; i++) sweep(spin,N,beta,h);

  // now sweep through lattice nsweep times
  // we will average the magnitization of the system
  // over nsweeps to get the average "thermalized" value
  nmag = ntotal = 0;
  for(int n=0; n<nsweep; n++) {
    sweep(spin,N,beta,h);
    for(int ns=1; ns<=N; ns++) {
      nmag += spin[ns];
      ntotal++;
      if (VisualDisplay) {
        if(spin[ns]==1) printf("+");
        if(spin[ns]==-1) printf("-");
      }
    }
    if (VisualDisplay) printf("\n");
   }

   printf("Magnetization of thermalized system: <s> = %lf\n", (double)nmag/ntotal);
   expected = exp(beta)*sinh(h) /
     sqrt( exp(2*beta)*sinh(h)*sinh(h) + exp(-2*beta) );
   printf("Expected result for infinite length chain over long time: %lf\n",
	  expected);

   return(0);
}