8 9 7 6 4 2 3 5 1 j x j j j + 1 l j l j = x j+1 x j, n x n x 1 = n 1 l j j=1 H j j + 1 l j l j E
a n 1 H = ae l j, j=1 l j = x j+1 x j, x n x 1 = n 1 j=1 l j, l j = ±l l > 0) n 1 H = ϵ l j, j=1 ϵ e x x = 0 e x = j=0 x j j! = 1 + x + 1 2 x2 + 1 6 x3 + + 1 j! xj +, 0! = 1 x = 1 e = 1 + 1 + 1 2 + 1 6 + + 1 j! + = j=0 1 j!,
e = n j=0 1 j!. n 1 15 n T j H j j Ψ j k Ψ j = 1 Z H ) j, Z = j H ) j, j T a b P a b H ) a P ab = H ) b P ba,
N N H 3N r N = r 1, r 2, r N ) 3N p N = p 1, p 2, p N ) H = N j=1 p 2 j 2m + Ur 1, r 2,, r N ), U r j = x j, y j, z j ) p j = p xj, p yj, p zj ) x j y j z j j x y z p xj p yj p zj p N r N r N = r 1, r 2,, r N ) p N = p 1, p 2,, p N ) r + r N = r 1 + r 1, r 2 + r 2,, r N + r N ) p + p N = p 1 + p 1, p 2 + p 2,, p N + p N ) H ) r N p N H ) Q r N Q < Q > Q Q H ) r N p N k < Q >= T H ) = r N p N U ) r N ), r N U N r N k t k t t m N 1 a x a y b z b
x a x a + δ x 1 2ξ 1 ), y a y a + δ y 1 2ξ 2 ), z a z a + δ z 1 2ξ 3 ). δ x, δ y, δ z x y z ξ j [0, 1] 8δ x δ y δ z t m t m t U t m U m U t U m k + 1 m U t < U m m [0, 1] η η < U ) m U t m k + 1 η > U ) m U t t k + 1 N l j = ±1 0
t m N l j 1 k t t m 1 n 1 a l a +1 1 l a +1 1 1 +1 l a 1 t m m t H t m H m H t H m k + 1 m H t < H m [0, 1] η m k + 1 η < H ) m H t η > H ) m H t t k + 1
1 #include <stdio.h> 2 #include <stdlib.h> 3 4 #define N 20) 5 6 double uran){ 7 returndouble)rand)/rand_max); 8 } 9 10 int main) 11 { 12 float rn[n]; 13 int l[n]; 14 int j; 15 16 forj=0;j<n;j++){ 17 rn[j] = uran); 18 ifrn[j] < 0.5){ 19 l[j] = -1; 20 }else{ 21 l[j] = 1; 22 } 23 24 printf"%f,%2d\n",rn[j],l[j]); 25 } 26 return0); 27 } 28
H ϵ = 6.0 10 22 k = 1.380662 23 1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 5 #define N 100) 6 #define EP 6.0e-22) /* [J] */ 7 #define KB 1.380662e-23) /* The Botzmann's constant [J/K] */ 8 #define MCTIME 1000) 9 10 double uran){ 11 returndouble)rand)/rand_max); 12 } 13 14 int main) 15 { 16 int l[n],p[n+1]; 17 int k,j,mct,lbefore,end; 18 float dene,t; 19 FILE *fp; 20 21 printf"temperature ="); 22 scanf"%f",&t); 23 printf"\n"); 24 25 forj=0;j<n;j++){ 26 l[j] = 1; 27 } 28 29 forj=0;j<n+1; j++){ 30 p[j] = 0;
31 } 32 33 formct=0;mct<mctime;mct++){ 34 fork=0;k<n;k++){ 35 j = N+1)*uran); 36 lbefore = l[j]; 37 l[j] *= -1); 38 dene = - EP*l[j] - -EP*lbefore); 39 ifdene > 0){ 40 ifexp-dene/kb*t)) < uran)){ 41 l[j] = lbefore; 42 } 43 } 44 } 45 46 forj=0,end=0;j<n;j++){ 47 end += l[j]; 48 } 49 p[end+n)/2]++; 50 51 printf"%d\n",mct); 52 } 53 54 fp = fopen"result.txt","w"); 55 56 forj=0;j<n+1);j++){ 57 fprintffp,"%d\t%e\n",2*j-n,float)p[j])/mctime); 58 } 59 60 fclosefp); 61 62 return0); 63 } 64
65 result.txt $ mv result.txt 10.txt result.txt 10.txt 100.txt 1000.txt
[0, 1] uran) double j int l[n] [ ]... [ ] p[n+1] [ ]... [ ] lbefor end float dene t FILE fp t for for { }
j 0 l[j] = 1 j j + 1 j j < N l[j] j 0 N p[j] 0 p[j] mct { } mct 1 mct MCTIME N N 1 j 0 N l[j] l[j] 1 1 lbefore l[j] l[j] l[j] 1 1 1 +1 l[j] 1 l[j] l[j] *= -1); l[j] = l[j]*-1); l j l j l j H H H H = ϵl 0 + + l j + + l N 1 ) ϵ)l 0 + + l j + + l N 1 ) = ϵl j ϵl j), dene
dene H ) H, H H [0, 1] l[j] x l j 1 x = N x = N 2 x 2 x = N, N + 2, N + 4,..., N 4, N 2, N j = x+n)/2 x = N, N +2, N +4,..., N 4, N 2, N j = 0, 1, 2,..., N 2, N 1, N x = 2j N p[j] 1 p[j] x p[x 1 x = 2j N result.txt 2j N p[j]/mctime n 1 H = ϵ l xj, j=1
l xj j x rothex.c 1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <GL/glut.h> 4 #include <math.h> 5 #include <time.h> 6 7 #define WIDTH 600 8 #define HIGHT 600 9 #define RADIUS 100 10 #define LENGTH 3 11 #define X0 300 12 #define Y0 300 13 #define DEGSTEP 10 14 15 #define N 100) 16 #define EP 6.0e-22) /* [J] */ 17 #define KB 1.380662e-23) /* The Botzmann's constant [J/K] */ 18 #define MCTIME 1000) 19 #define DL 0.5) 20 21 void displayvoid); 22 void reshapeint w, int h); 23 void keyboardunsigned char key, int x, int y); 24 void ginitint* pargc, char** argv); 25 void idlevoid); 26 27 int Degshift=0; 28 int Zoom=1,Stop=1; 29 float Ep,T, Lx[N], Ly[N]; 30 31
32 double uran){ 33 returndouble)rand)/rand_max); 34 } 35 36 37 int mcrun) 38 { 39 int j; 40 float dene; 41 float lxbefore, lybefore; 42 float dx, dy, r; 43 44 j = N-1)*uran)+1; 45 lxbefore = Lx[j]; 46 lybefore = Ly[j]; 47 48 do{ 49 dx = 2.0*DL*uran)-0.5); 50 dy = 2.0*DL*uran)-0.5); 51 }while dx*dx + dy*dy > DL*DL); 52 53 /* 54 Lx[j] = Lx[j] + dx; 55 Ly[j] = Ly[j] + dy; 56 */ 57 58 Lx[j] = dx; 59 Ly[j] = dy; 60 61 r = sqrtlx[j]*lx[j] + Ly[j]*Ly[j]); 62 Lx[j] = Lx[j]/r; 63 Ly[j] = Ly[j]/r; 64 65 dene = - Ep*Lx[j] - -Ep*lxbefore);
66 dene = dene*ep; 67 68 ifdene > 0){ 69 ifexp-dene/kb*t)) < uran)){ 70 Lx[j] = lxbefore; 71 Ly[j] = lybefore; 72 } 73 } 74 75 return0); 76 } 77 78 void idlevoid) 79 { 80 ifstop == 1)return; 81 mcrun); 82 glutpostredisplay); 83 } 84 85 void displayvoid) 86 { 87 int x, y, j; 88 89 glcleargl_color_buffer_bit); 90 glpointsize1); 91 glbegingl_line_strip); 92 glcolor3f1.0, 1.0, 1.0); /* white */ 93 x = 0; 94 y = 0; 95 glvertex2fx0+x, Y0+y); 96 forj=1;j<n;j++){ 97 ifj%2 == 0){ 98 glcolor3f1.0, 1.0, 1.0); 99 }else{
100 glcolor3f1.0, 1.0, 0.0); 101 } 102 x = x+zoom*length*lx[j]; 103 y = y+zoom*length*ly[j]; 104 glvertex2fx0+x, Y0+y); 105 } 106 107 glend); 108 glutswapbuffers); 109 } 110 111 void reshapeint w, int h) 112 { 113 glviewport0, 0, GLsizei)w, GLsizei)h ); 114 glmatrixmodegl_projection); 115 glloadidentity); 116 gluortho2d0.0, GLdouble)w, 0.0, GLdouble) h); 117 } 118 119 void keyboardunsigned char key, int x, int y) 120 { 121 switchkey) { 122 case 's': 123 Stop *= -1; 124 break; 125 case 'i': 126 ifzoom == 0){ 127 Zoom=1; 128 }else{ 129 Zoom *= 2; 130 } 131 break; 132 case 'o': 133 Zoom /= 2;
134 break; 135 case 'e': 136 printf"\nelectric field ="); 137 scanf"%f",&ep); 138 Stop = 1; 139 break; 140 case 't': 141 printf"\ntemperature ="); 142 scanf"%f",&t); 143 Stop = 1; 144 break; 145 case 'q': 146 case 'Q': 147 exit0); 148 default: 149 break; 150 } 151 } 152 153 void ginitint* pargc, char** argv) 154 { 155 glutinitpargc,argv); 156 glutinitdisplaymodeglut_double GLUT_RGB); 157 glutinitwindowsizewidth,hight); 158 glutinitwindowposition100,50); 159 glutcreatewindow"numerical Recipes for Polymer Science"); 160 glclearcolor0.0,0.0,0.0,0.0); /* background color */ 161 glshademodelgl_flat); 162 glutdisplayfuncdisplay); 163 glutreshapefuncreshape); 164 glutkeyboardfunckeyboard); 165 glutidlefuncidle); 166 } 167
168 int mainint argc, char* argv[]) 169 { 170 int j; 171 ginit&argc,argv); 172 173 printf"temperature ="); 174 scanf"%f",&t); 175 printf"\n"); 176 printf"electric field ="); 177 scanf"%f",&ep); 178 179 forj=1;j<n;j++){ 180 Lx[j] = 0.0; 181 Ly[j] = 1.0; 182 } 183 184 glutmainloop); 185 return 0; 186 } 187