編集履歴

質問編集履歴

書式の改善を行いました

2018/05/11 11:51

投稿

astromelt0416

スコア15

test CHANGED Viewed

File without changes

test CHANGED Viewed

@@ -134,11 +134,17 @@
+# Intitialize the XY network
 def Init():
     return np.random.rand(L, L)*2*np.pi
+    #return np.ones([L, L])
+# periodic neighbor
 def next(x):
@@ -152,6 +158,8 @@
+# construct the bond lattice
 def FreezeBonds(Ising,T,S):
     iBondFrozen = np.zeros([L,L],dtype=int)
@@ -178,6 +186,8 @@
+# H-K algorithm to identify clusters
 def properlabel(prp_label,i):
     while prp_label[i] != i:
@@ -188,6 +198,8 @@
+# Swendsen-Wang cluster
 def clusterfind(iBondFrozen,jBondFrozen):
     cluster = np.zeros([L, L],dtype=int)
@@ -208,6 +220,8 @@
+            # check to (i-1,j)
             if (i > 0) and iBondFrozen[i-1][j]:
                 ibonds[bonds] = i-1
@@ -216,7 +230,7 @@
                 bonds += 1
+            # (i,j) at i edge, check to (i+1,j)
             if (i == L-1) and iBondFrozen[i][j]:
@@ -226,7 +240,7 @@
                 bonds += 1
+            # check to (i,j-1)
             if (j > 0) and jBondFrozen[i][j-1]:
@@ -236,7 +250,7 @@
                 bonds += 1
+            # (i,j) at j edge, check to (i,j+1)
             if (j == L-1) and jBondFrozen[i][j]:
@@ -248,7 +262,7 @@
+            # check and label clusters
             if bonds == 0:
@@ -274,7 +288,7 @@
                 cluster[i][j] = minlabel
+                # link to the previous labels
                 for b in np.arange(bonds):
@@ -282,7 +296,7 @@
                     prp_label[plabel_n] = minlabel
+                    # re-set the labels on connected sites
                     cluster[ibonds[b]][jbonds[b]] = minlabel
@@ -290,7 +304,7 @@
+# flip the cluster spins
 def flipCluster(Ising,cluster,prp_label):
@@ -298,7 +312,7 @@
         for j in np.arange(L):
+            # relabel all the cluster labels with the right ones
             cluster[i][j] = properlabel(prp_label,cluster[i][j])
@@ -316,7 +330,7 @@
             randn = np.random.rand()
+            # mark the flipped label, use this label to flip all the cluster elements with this label
             if (not sNewChosen[label]) and randn < 0.5:
@@ -344,7 +358,9 @@
+# Swendsen-Wang Algorithm in Ising model (with coupling constant dependency on sites)
+# One-step for Ising
 def oneMCstepIsing(Ising, S):
@@ -358,7 +374,7 @@
+# Decompose XY network to two Ising networks with project direction proj
 def decompose(XY,proj):
@@ -382,7 +398,7 @@
+# Compose two Ising networks to XY network
 def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
@@ -408,7 +424,11 @@
     Isingx_new = oneMCstepIsing(Isingx, S_x)
+    #Isingy_new = oneMCstepIsing(Isingy, S_y)
+    # Here we only use the flopping of Ising in projection direction, without the perpendicular direction
+    #XY_new = compose(Isingx_new, Isingy_new, proj, S_x, S_y)
     XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
@@ -416,7 +436,7 @@
+# Calculate the energy for XY network
 def EnMag(XY):
@@ -426,7 +446,7 @@
         for j in np.arange(L):
+            # energy
             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]))
@@ -440,19 +460,19 @@
+# Swendsen Wang method for XY model
 def SWang(T):
     XY = Init()
+    # thermal steps to get the equilibrium
     for step in np.arange(ESTEP):
         XY = oneMCstepXY(XY)
+    # finish with thermal equilibrium, and begin to calc observables
     E_sum = 0
@@ -518,6 +538,8 @@
+# plot the figures
 T = Trange
@@ -582,4 +604,6 @@
 ```

書式を改善しました

2018/05/11 11:51

投稿

astromelt0416

スコア15

test CHANGED Viewed

File without changes

test CHANGED Viewed

@@ -104,6 +104,10 @@
+```
 ```python
 import matplotlib
@@ -576,10 +580,6 @@
 np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
 ```
-```

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

2018/05/11 11:47

投稿

astromelt0416

スコア15

test CHANGED Viewed

File without changes

test CHANGED Viewed

@@ -104,4 +104,482 @@
+```python
+import matplotlib
+matplotlib.use('Agg')
+import numpy as np
+from numpy import linalg as LA
+import matplotlib.pyplot as plt
+L = 8
+ESTEP = 1000
+STEP = 10000
+J = 1 # J>0 to make it ferromagnetic
+def Init():
+    return np.random.rand(L, L)*2*np.pi
+def next(x):
+    if x == L-1:
+        return 0
+    else:
+        return x+1
+def FreezeBonds(Ising,T,S):
+    iBondFrozen = np.zeros([L,L],dtype=int)
+    jBondFrozen = np.zeros([L,L],dtype=int)
+    for i in np.arange(L):
+        for j in np.arange(L):
+            freezProb_nexti = 1 - np.exp(-2 * J * S[i][j] * S[next(i)][j] / T)
+            freezProb_nextj = 1 - np.exp(-2 * J * S[i][j] * S[i][next(j)] / T)
+            if (Ising[i][j] == Ising[next(i)][j]) and (np.random.rand() < freezProb_nexti):
+                iBondFrozen[i][j] = 1
+            if (Ising[i][j] == Ising[i][next(j)]) and (np.random.rand() < freezProb_nextj):
+                jBondFrozen[i][j] = 1
+    return iBondFrozen, jBondFrozen
+def properlabel(prp_label,i):
+    while prp_label[i] != i:
+        i = prp_label[i]
+    return i
+def clusterfind(iBondFrozen,jBondFrozen):
+    cluster = np.zeros([L, L],dtype=int)
+    prp_label = np.zeros(L**2,dtype=int)
+    label = 0
+    for i in np.arange(L):
+        for j in np.arange(L):
+            bonds = 0
+            ibonds = np.zeros(4,dtype=int)
+            jbonds = np.zeros(4,dtype=int)
+            if (i > 0) and iBondFrozen[i-1][j]:
+                ibonds[bonds] = i-1
+                jbonds[bonds] = j
+                bonds += 1
+            if (i == L-1) and iBondFrozen[i][j]:
+                ibonds[bonds] = 0
+                jbonds[bonds] = j
+                bonds += 1
+            if (j > 0) and jBondFrozen[i][j-1]:
+                ibonds[bonds] = i
+                jbonds[bonds] = j-1
+                bonds += 1
+            if (j == L-1) and jBondFrozen[i][j]:
+                ibonds[bonds] = i
+                jbonds[bonds] = 0
+                bonds += 1
+            if bonds == 0:
+                cluster[i][j] = label
+                prp_label[label] = label
+                label += 1
+            else:
+                minlabel = label
+                for b in np.arange(bonds):
+                    plabel = properlabel(prp_label,cluster[ibonds[b]][jbonds[b]])
+                    if minlabel > plabel:
+                        minlabel = plabel
+                cluster[i][j] = minlabel
+                for b in np.arange(bonds):
+                    plabel_n = cluster[ibonds[b]][jbonds[b]]
+                    prp_label[plabel_n] = minlabel
+                    cluster[ibonds[b]][jbonds[b]] = minlabel
+    return cluster, prp_label
+def flipCluster(Ising,cluster,prp_label):
+    for i in np.arange(L):
+        for j in np.arange(L):
+            cluster[i][j] = properlabel(prp_label,cluster[i][j])
+    sNewChosen = np.zeros(L**2)
+    sNew = np.zeros(L**2)
+    flips = 0 # get the number of flipped spins to calculate the Endiff and Magdiff
+    for i in np.arange(L):
+        for j in np.arange(L):
+            label = cluster[i][j]
+            randn = np.random.rand()
+            if (not sNewChosen[label]) and randn < 0.5:
+                sNew[label] = +1
+                sNewChosen[label] = True
+            elif (not sNewChosen[label]) and randn >= 0.5:
+                sNew[label] = -1
+                sNewChosen[label] = True
+            if Ising[i][j] != sNew[label]:
+                Ising[i][j] = sNew[label]
+                flips += 1
+    return Ising,flips
+def oneMCstepIsing(Ising, S):
+    [iBondFrozen, jBondFrozen] = FreezeBonds(Ising, T, S)
+    [SWcluster, prp_label] = clusterfind(iBondFrozen, jBondFrozen)
+    [Ising, flips] = flipCluster(Ising, SWcluster, prp_label)
+    return Ising
+def decompose(XY,proj):
+    x = np.cos(XY)
+    y = np.sin(XY)
+    x_rot = np.multiply(x,np.cos(proj))+np.multiply(y,np.sin(proj))
+    y_rot = -np.multiply(x,np.sin(proj))+np.multiply(y,np.cos(proj))
+    Isingx = np.sign(x_rot)
+    Isingy = np.sign(y_rot)
+    S_x = np.absolute(x_rot)
+    S_y = np.absolute(y_rot)
+    return Isingx, Isingy, S_x, S_y
+def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
+    x_rot_new = np.multiply(Isingx_new,S_x)
+    y_rot_new = np.multiply(Isingy_new,S_y)
+    x_new = np.multiply(x_rot_new,np.cos(proj))-np.multiply(y_rot_new,np.sin(proj))
+    y_new = np.multiply(x_rot_new,np.sin(proj))+np.multiply(y_rot_new,np.cos(proj))
+    XY_new = np.arctan2(y_new,x_new)
+    return XY_new
-ここにより詳細な情報を記載してください。
+def oneMCstepXY(XY):
+    proj = np.random.rand()
+    [Isingx, Isingy, S_x, S_y] = decompose(XY, proj)
+    Isingx_new = oneMCstepIsing(Isingx, S_x)
+    XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
+    return XY_new
+def EnMag(XY):
+    energy = 0
+    for i in np.arange(L):
+        for j in np.arange(L):
+            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]))
+    magx = np.sum(np.cos(XY))
+    magy = np.sum(np.sin(XY))
+    mag = np.array([magx,magy])
+    return energy * 0.5, LA.norm(mag)/(L**2)
+def SWang(T):
+    XY = Init()
+    for step in np.arange(ESTEP):
+        XY = oneMCstepXY(XY)
+    E_sum = 0
+    M_sum = 0
+    Esq_sum = 0
+    Msq_sum = 0
+    for step in np.arange(STEP):
+        XY = oneMCstepXY(XY)
+        [E,M] = EnMag(XY)
+        E_sum += E
+        M_sum += M
+        Esq_sum += E**2
+        Msq_sum += M**2
+    E_mean = E_sum/STEP/(L**2)
+    M_mean = M_sum/STEP
+    Esq_mean = Esq_sum/STEP/(L**4)
+    Msq_mean = Msq_sum/STEP
+    return XY, E_mean, M_mean, Esq_mean, Msq_mean
+M = np.array([])fig = plt.gcf()
+E = np.array([])
+M_sus = np.array([])
+SpcH = np.array([])
+Trange = np.linspace(0.1, 2.5, 10)
+for T in Trange:
+    [Ising, E_mean, M_mean, Esq_mean, Msq_mean] = SWang(T)
+    M = np.append(M, np.abs(M_mean))
+    E = np.append(E, E_mean)
+    M_sus = np.append(M_sus, 1/T*(Msq_mean-M_mean**2))
+    SpcH = np.append(SpcH, 1/T**2*(Esq_mean-E_mean**2))
+T = Trange
+plt.figure()
+plt.plot(T, E, 'rx-')
+plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
+plt.ylabel(r'$\langle E \rangle$ per site $(J)$')
+plt.savefig("E8.pdf", format='pdf', bbox_inches='tight')
+plt.figure()
+plt.plot(T, SpcH, 'kx-')
+plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
+plt.ylabel(r'$C_V$ per site $(\frac{J^2}{k_B^2})$')
+plt.savefig("Cv8.pdf", format='pdf', bbox_inches='tight')
+plt.figure()
+plt.plot(T, M, 'bx-')
+plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
+plt.ylabel(r'$\langle|M|\rangle$ per site $(\mu)$')
+plt.savefig("M8.pdf", format='pdf', bbox_inches='tight')
+plt.figure()
+plt.plot(T, M_sus, 'gx-')
+plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
+plt.ylabel(r'$\chi$ $(\frac{\mu}{k_B})$')
+plt.savefig("chi8.pdf", format='pdf', bbox_inches='tight')
+plt.tight_layout()
+fig = plt.gcf()
+plt.show()
+np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
+```
+```

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

2018/05/11 11:45

投稿

astromelt0416

スコア15

test CHANGED Viewed

File without changes

test CHANGED Viewed

File without changes