質問編集履歴

4

書式の改善を行いました

2018/05/11 11:51

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -134,11 +134,17 @@
134
134
 
135
135
 
136
136
 
137
+ # Intitialize the XY network
138
+
137
139
  def Init():
138
140
 
139
141
  return np.random.rand(L, L)*2*np.pi
140
142
 
141
-
143
+ #return np.ones([L, L])
144
+
145
+
146
+
147
+ # periodic neighbor
142
148
 
143
149
  def next(x):
144
150
 
@@ -152,6 +158,8 @@
152
158
 
153
159
 
154
160
 
161
+ # construct the bond lattice
162
+
155
163
  def FreezeBonds(Ising,T,S):
156
164
 
157
165
  iBondFrozen = np.zeros([L,L],dtype=int)
@@ -178,6 +186,8 @@
178
186
 
179
187
 
180
188
 
189
+ # H-K algorithm to identify clusters
190
+
181
191
  def properlabel(prp_label,i):
182
192
 
183
193
  while prp_label[i] != i:
@@ -188,6 +198,8 @@
188
198
 
189
199
 
190
200
 
201
+ # Swendsen-Wang cluster
202
+
191
203
  def clusterfind(iBondFrozen,jBondFrozen):
192
204
 
193
205
  cluster = np.zeros([L, L],dtype=int)
@@ -208,6 +220,8 @@
208
220
 
209
221
 
210
222
 
223
+ # check to (i-1,j)
224
+
211
225
  if (i > 0) and iBondFrozen[i-1][j]:
212
226
 
213
227
  ibonds[bonds] = i-1
@@ -216,7 +230,7 @@
216
230
 
217
231
  bonds += 1
218
232
 
219
-
233
+ # (i,j) at i edge, check to (i+1,j)
220
234
 
221
235
  if (i == L-1) and iBondFrozen[i][j]:
222
236
 
@@ -226,7 +240,7 @@
226
240
 
227
241
  bonds += 1
228
242
 
229
-
243
+ # check to (i,j-1)
230
244
 
231
245
  if (j > 0) and jBondFrozen[i][j-1]:
232
246
 
@@ -236,7 +250,7 @@
236
250
 
237
251
  bonds += 1
238
252
 
239
-
253
+ # (i,j) at j edge, check to (i,j+1)
240
254
 
241
255
  if (j == L-1) and jBondFrozen[i][j]:
242
256
 
@@ -248,7 +262,7 @@
248
262
 
249
263
 
250
264
 
251
-
265
+ # check and label clusters
252
266
 
253
267
  if bonds == 0:
254
268
 
@@ -274,7 +288,7 @@
274
288
 
275
289
  cluster[i][j] = minlabel
276
290
 
277
-
291
+ # link to the previous labels
278
292
 
279
293
  for b in np.arange(bonds):
280
294
 
@@ -282,7 +296,7 @@
282
296
 
283
297
  prp_label[plabel_n] = minlabel
284
298
 
285
-
299
+ # re-set the labels on connected sites
286
300
 
287
301
  cluster[ibonds[b]][jbonds[b]] = minlabel
288
302
 
@@ -290,7 +304,7 @@
290
304
 
291
305
 
292
306
 
293
-
307
+ # flip the cluster spins
294
308
 
295
309
  def flipCluster(Ising,cluster,prp_label):
296
310
 
@@ -298,7 +312,7 @@
298
312
 
299
313
  for j in np.arange(L):
300
314
 
301
-
315
+ # relabel all the cluster labels with the right ones
302
316
 
303
317
  cluster[i][j] = properlabel(prp_label,cluster[i][j])
304
318
 
@@ -316,7 +330,7 @@
316
330
 
317
331
  randn = np.random.rand()
318
332
 
319
-
333
+ # mark the flipped label, use this label to flip all the cluster elements with this label
320
334
 
321
335
  if (not sNewChosen[label]) and randn < 0.5:
322
336
 
@@ -344,7 +358,9 @@
344
358
 
345
359
 
346
360
 
347
-
361
+ # Swendsen-Wang Algorithm in Ising model (with coupling constant dependency on sites)
362
+
363
+ # One-step for Ising
348
364
 
349
365
  def oneMCstepIsing(Ising, S):
350
366
 
@@ -358,7 +374,7 @@
358
374
 
359
375
 
360
376
 
361
-
377
+ # Decompose XY network to two Ising networks with project direction proj
362
378
 
363
379
  def decompose(XY,proj):
364
380
 
@@ -382,7 +398,7 @@
382
398
 
383
399
 
384
400
 
385
-
401
+ # Compose two Ising networks to XY network
386
402
 
387
403
  def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
388
404
 
@@ -408,7 +424,11 @@
408
424
 
409
425
  Isingx_new = oneMCstepIsing(Isingx, S_x)
410
426
 
411
-
427
+ #Isingy_new = oneMCstepIsing(Isingy, S_y)
428
+
429
+ # Here we only use the flopping of Ising in projection direction, without the perpendicular direction
430
+
431
+ #XY_new = compose(Isingx_new, Isingy_new, proj, S_x, S_y)
412
432
 
413
433
  XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
414
434
 
@@ -416,7 +436,7 @@
416
436
 
417
437
 
418
438
 
419
-
439
+ # Calculate the energy for XY network
420
440
 
421
441
  def EnMag(XY):
422
442
 
@@ -426,7 +446,7 @@
426
446
 
427
447
  for j in np.arange(L):
428
448
 
429
-
449
+ # energy
430
450
 
431
451
  energy = energy - (np.cos(XY[i][j]-XY[(i-1)%L][j])+np.cos(XY[i][j]-XY[(i+1)%L][j])+np.cos(XY[i][j]-XY[i][(j-1)%L])+np.cos(XY[i][j]-XY[i][(j+1)%L]))
432
452
 
@@ -440,19 +460,19 @@
440
460
 
441
461
 
442
462
 
443
-
463
+ # Swendsen Wang method for XY model
444
464
 
445
465
  def SWang(T):
446
466
 
447
467
  XY = Init()
448
468
 
449
-
469
+ # thermal steps to get the equilibrium
450
470
 
451
471
  for step in np.arange(ESTEP):
452
472
 
453
473
  XY = oneMCstepXY(XY)
454
474
 
455
-
475
+ # finish with thermal equilibrium, and begin to calc observables
456
476
 
457
477
  E_sum = 0
458
478
 
@@ -518,6 +538,8 @@
518
538
 
519
539
 
520
540
 
541
+ # plot the figures
542
+
521
543
  T = Trange
522
544
 
523
545
 
@@ -582,4 +604,6 @@
582
604
 
583
605
 
584
606
 
607
+
608
+
585
609
  ```

3

書式を改善しました

2018/05/11 11:51

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -104,6 +104,10 @@
104
104
 
105
105
 
106
106
 
107
+ ```
108
+
109
+
110
+
107
111
  ```python
108
112
 
109
113
  import matplotlib
@@ -576,10 +580,6 @@
576
580
 
577
581
  np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
578
582
 
583
+
584
+
579
585
  ```
580
-
581
-
582
-
583
-
584
-
585
- ```

2

コードの全文を貼りました

2018/05/11 11:47

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -104,4 +104,482 @@
104
104
 
105
105
 
106
106
 
107
+ ```python
108
+
109
+ import matplotlib
110
+
111
+ matplotlib.use('Agg')
112
+
113
+ import numpy as np
114
+
115
+ from numpy import linalg as LA
116
+
117
+ import matplotlib.pyplot as plt
118
+
119
+
120
+
121
+ L = 8
122
+
123
+ ESTEP = 1000
124
+
125
+ STEP = 10000
126
+
127
+
128
+
129
+ J = 1 # J>0 to make it ferromagnetic
130
+
131
+
132
+
133
+ def Init():
134
+
135
+ return np.random.rand(L, L)*2*np.pi
136
+
137
+
138
+
139
+ def next(x):
140
+
141
+ if x == L-1:
142
+
143
+ return 0
144
+
145
+ else:
146
+
147
+ return x+1
148
+
149
+
150
+
151
+ def FreezeBonds(Ising,T,S):
152
+
153
+ iBondFrozen = np.zeros([L,L],dtype=int)
154
+
155
+ jBondFrozen = np.zeros([L,L],dtype=int)
156
+
157
+ for i in np.arange(L):
158
+
159
+ for j in np.arange(L):
160
+
161
+ freezProb_nexti = 1 - np.exp(-2 * J * S[i][j] * S[next(i)][j] / T)
162
+
163
+ freezProb_nextj = 1 - np.exp(-2 * J * S[i][j] * S[i][next(j)] / T)
164
+
165
+ if (Ising[i][j] == Ising[next(i)][j]) and (np.random.rand() < freezProb_nexti):
166
+
167
+ iBondFrozen[i][j] = 1
168
+
169
+ if (Ising[i][j] == Ising[i][next(j)]) and (np.random.rand() < freezProb_nextj):
170
+
171
+ jBondFrozen[i][j] = 1
172
+
173
+ return iBondFrozen, jBondFrozen
174
+
175
+
176
+
177
+ def properlabel(prp_label,i):
178
+
179
+ while prp_label[i] != i:
180
+
181
+ i = prp_label[i]
182
+
183
+ return i
184
+
185
+
186
+
187
+ def clusterfind(iBondFrozen,jBondFrozen):
188
+
189
+ cluster = np.zeros([L, L],dtype=int)
190
+
191
+ prp_label = np.zeros(L**2,dtype=int)
192
+
193
+ label = 0
194
+
195
+ for i in np.arange(L):
196
+
197
+ for j in np.arange(L):
198
+
199
+ bonds = 0
200
+
201
+ ibonds = np.zeros(4,dtype=int)
202
+
203
+ jbonds = np.zeros(4,dtype=int)
204
+
205
+
206
+
207
+ if (i > 0) and iBondFrozen[i-1][j]:
208
+
209
+ ibonds[bonds] = i-1
210
+
211
+ jbonds[bonds] = j
212
+
213
+ bonds += 1
214
+
215
+
216
+
217
+ if (i == L-1) and iBondFrozen[i][j]:
218
+
219
+ ibonds[bonds] = 0
220
+
221
+ jbonds[bonds] = j
222
+
223
+ bonds += 1
224
+
225
+
226
+
227
+ if (j > 0) and jBondFrozen[i][j-1]:
228
+
229
+ ibonds[bonds] = i
230
+
231
+ jbonds[bonds] = j-1
232
+
233
+ bonds += 1
234
+
235
+
236
+
237
+ if (j == L-1) and jBondFrozen[i][j]:
238
+
239
+ ibonds[bonds] = i
240
+
241
+ jbonds[bonds] = 0
242
+
243
+ bonds += 1
244
+
245
+
246
+
247
+
248
+
249
+ if bonds == 0:
250
+
251
+ cluster[i][j] = label
252
+
253
+ prp_label[label] = label
254
+
255
+ label += 1
256
+
257
+ else:
258
+
259
+ minlabel = label
260
+
261
+ for b in np.arange(bonds):
262
+
263
+ plabel = properlabel(prp_label,cluster[ibonds[b]][jbonds[b]])
264
+
265
+ if minlabel > plabel:
266
+
267
+ minlabel = plabel
268
+
269
+
270
+
271
+ cluster[i][j] = minlabel
272
+
273
+
274
+
275
+ for b in np.arange(bonds):
276
+
277
+ plabel_n = cluster[ibonds[b]][jbonds[b]]
278
+
279
+ prp_label[plabel_n] = minlabel
280
+
281
+
282
+
283
+ cluster[ibonds[b]][jbonds[b]] = minlabel
284
+
285
+ return cluster, prp_label
286
+
287
+
288
+
289
+
290
+
291
+ def flipCluster(Ising,cluster,prp_label):
292
+
293
+ for i in np.arange(L):
294
+
295
+ for j in np.arange(L):
296
+
297
+
298
+
299
+ cluster[i][j] = properlabel(prp_label,cluster[i][j])
300
+
301
+ sNewChosen = np.zeros(L**2)
302
+
303
+ sNew = np.zeros(L**2)
304
+
305
+ flips = 0 # get the number of flipped spins to calculate the Endiff and Magdiff
306
+
307
+ for i in np.arange(L):
308
+
309
+ for j in np.arange(L):
310
+
311
+ label = cluster[i][j]
312
+
313
+ randn = np.random.rand()
314
+
315
+
316
+
317
+ if (not sNewChosen[label]) and randn < 0.5:
318
+
319
+ sNew[label] = +1
320
+
321
+ sNewChosen[label] = True
322
+
323
+ elif (not sNewChosen[label]) and randn >= 0.5:
324
+
325
+ sNew[label] = -1
326
+
327
+ sNewChosen[label] = True
328
+
329
+
330
+
331
+ if Ising[i][j] != sNew[label]:
332
+
333
+ Ising[i][j] = sNew[label]
334
+
335
+ flips += 1
336
+
337
+
338
+
339
+ return Ising,flips
340
+
341
+
342
+
343
+
344
+
345
+ def oneMCstepIsing(Ising, S):
346
+
347
+ [iBondFrozen, jBondFrozen] = FreezeBonds(Ising, T, S)
348
+
349
+ [SWcluster, prp_label] = clusterfind(iBondFrozen, jBondFrozen)
350
+
351
+ [Ising, flips] = flipCluster(Ising, SWcluster, prp_label)
352
+
353
+ return Ising
354
+
355
+
356
+
357
+
358
+
359
+ def decompose(XY,proj):
360
+
361
+ x = np.cos(XY)
362
+
363
+ y = np.sin(XY)
364
+
365
+ x_rot = np.multiply(x,np.cos(proj))+np.multiply(y,np.sin(proj))
366
+
367
+ y_rot = -np.multiply(x,np.sin(proj))+np.multiply(y,np.cos(proj))
368
+
369
+ Isingx = np.sign(x_rot)
370
+
371
+ Isingy = np.sign(y_rot)
372
+
373
+ S_x = np.absolute(x_rot)
374
+
375
+ S_y = np.absolute(y_rot)
376
+
377
+ return Isingx, Isingy, S_x, S_y
378
+
379
+
380
+
381
+
382
+
383
+ def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
384
+
385
+ x_rot_new = np.multiply(Isingx_new,S_x)
386
+
387
+ y_rot_new = np.multiply(Isingy_new,S_y)
388
+
389
+ x_new = np.multiply(x_rot_new,np.cos(proj))-np.multiply(y_rot_new,np.sin(proj))
390
+
391
+ y_new = np.multiply(x_rot_new,np.sin(proj))+np.multiply(y_rot_new,np.cos(proj))
392
+
393
+ XY_new = np.arctan2(y_new,x_new)
394
+
395
+ return XY_new
396
+
397
+
398
+
107
- ここにより詳細な情報を記載してください。
399
+ def oneMCstepXY(XY):
400
+
401
+ proj = np.random.rand()
402
+
403
+ [Isingx, Isingy, S_x, S_y] = decompose(XY, proj)
404
+
405
+ Isingx_new = oneMCstepIsing(Isingx, S_x)
406
+
407
+
408
+
409
+ XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
410
+
411
+ return XY_new
412
+
413
+
414
+
415
+
416
+
417
+ def EnMag(XY):
418
+
419
+ energy = 0
420
+
421
+ for i in np.arange(L):
422
+
423
+ for j in np.arange(L):
424
+
425
+
426
+
427
+ energy = energy - (np.cos(XY[i][j]-XY[(i-1)%L][j])+np.cos(XY[i][j]-XY[(i+1)%L][j])+np.cos(XY[i][j]-XY[i][(j-1)%L])+np.cos(XY[i][j]-XY[i][(j+1)%L]))
428
+
429
+ magx = np.sum(np.cos(XY))
430
+
431
+ magy = np.sum(np.sin(XY))
432
+
433
+ mag = np.array([magx,magy])
434
+
435
+ return energy * 0.5, LA.norm(mag)/(L**2)
436
+
437
+
438
+
439
+
440
+
441
+ def SWang(T):
442
+
443
+ XY = Init()
444
+
445
+
446
+
447
+ for step in np.arange(ESTEP):
448
+
449
+ XY = oneMCstepXY(XY)
450
+
451
+
452
+
453
+ E_sum = 0
454
+
455
+ M_sum = 0
456
+
457
+ Esq_sum = 0
458
+
459
+ Msq_sum = 0
460
+
461
+ for step in np.arange(STEP):
462
+
463
+ XY = oneMCstepXY(XY)
464
+
465
+ [E,M] = EnMag(XY)
466
+
467
+
468
+
469
+ E_sum += E
470
+
471
+ M_sum += M
472
+
473
+ Esq_sum += E**2
474
+
475
+ Msq_sum += M**2
476
+
477
+
478
+
479
+ E_mean = E_sum/STEP/(L**2)
480
+
481
+ M_mean = M_sum/STEP
482
+
483
+ Esq_mean = Esq_sum/STEP/(L**4)
484
+
485
+ Msq_mean = Msq_sum/STEP
486
+
487
+
488
+
489
+ return XY, E_mean, M_mean, Esq_mean, Msq_mean
490
+
491
+
492
+
493
+ M = np.array([])fig = plt.gcf()
494
+
495
+ E = np.array([])
496
+
497
+ M_sus = np.array([])
498
+
499
+ SpcH = np.array([])
500
+
501
+ Trange = np.linspace(0.1, 2.5, 10)
502
+
503
+ for T in Trange:
504
+
505
+ [Ising, E_mean, M_mean, Esq_mean, Msq_mean] = SWang(T)
506
+
507
+ M = np.append(M, np.abs(M_mean))
508
+
509
+ E = np.append(E, E_mean)
510
+
511
+ M_sus = np.append(M_sus, 1/T*(Msq_mean-M_mean**2))
512
+
513
+ SpcH = np.append(SpcH, 1/T**2*(Esq_mean-E_mean**2))
514
+
515
+
516
+
517
+ T = Trange
518
+
519
+
520
+
521
+ plt.figure()
522
+
523
+ plt.plot(T, E, 'rx-')
524
+
525
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
526
+
527
+ plt.ylabel(r'$\langle E \rangle$ per site $(J)$')
528
+
529
+ plt.savefig("E8.pdf", format='pdf', bbox_inches='tight')
530
+
531
+
532
+
533
+ plt.figure()
534
+
535
+ plt.plot(T, SpcH, 'kx-')
536
+
537
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
538
+
539
+ plt.ylabel(r'$C_V$ per site $(\frac{J^2}{k_B^2})$')
540
+
541
+ plt.savefig("Cv8.pdf", format='pdf', bbox_inches='tight')
542
+
543
+
544
+
545
+ plt.figure()
546
+
547
+ plt.plot(T, M, 'bx-')
548
+
549
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
550
+
551
+ plt.ylabel(r'$\langle|M|\rangle$ per site $(\mu)$')
552
+
553
+ plt.savefig("M8.pdf", format='pdf', bbox_inches='tight')
554
+
555
+
556
+
557
+ plt.figure()
558
+
559
+ plt.plot(T, M_sus, 'gx-')
560
+
561
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
562
+
563
+ plt.ylabel(r'$\chi$ $(\frac{\mu}{k_B})$')
564
+
565
+ plt.savefig("chi8.pdf", format='pdf', bbox_inches='tight')
566
+
567
+
568
+
569
+ plt.tight_layout()
570
+
571
+ fig = plt.gcf()
572
+
573
+ plt.show()
574
+
575
+
576
+
577
+ np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
578
+
579
+ ```
580
+
581
+
582
+
583
+
584
+
585
+ ```

1

初心者マークを入れ忘れていました。python初心者です。

2018/05/11 11:45

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
File without changes