Skip to content

Commit b746f0d

Browse files
Chengcheng-XiaoChengcheng-Xiao
committed
add Ising model, add input file interface
1 parent 11a0cd6 commit b746f0d

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

54 files changed

+1625
-492
lines changed

‎.gitignore

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,2 @@
11
**/.DS_Store
22
Docs
3-
ising

‎LICENSE.md

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
Copyright (c) 2019 Chengcheng XIAO, [email protected]
2+
3+
Permission is hereby granted, free of charge, to any person obtaining
4+
a copy of this software and associated documentation files (the
5+
"Software"), to deal in the Software without restriction, including
6+
without limitation the rights to use, copy, modify, merge, publish,
7+
distribute, sublicense, and/or sell copies of the Software, and to
8+
permit persons to whom the Software is furnished to do so, subject to
9+
the following conditions:
10+
11+
The above copyright notice and this permission notice shall be
12+
included in all copies or substantial portions of the Software.
13+
14+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15+
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16+
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17+
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18+
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19+
OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20+
WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

‎MC_MPI.py

Lines changed: 45 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,40 @@
1515
#----------------------------------------------------------------------
1616
## BLOCK OF FUNCTIONS USED IN THE MAIN CODE
1717
#----------------------------------------------------------------------
18+
def read_input():
19+
20+
''' a subroutine to get AO and MO number and number of kpoints and otehr information from input.woops'''
21+
22+
dataset={}
23+
file = open('input.MC', "r")
24+
#Default
25+
data = file.readlines()
26+
nt = 18 # number of temperature points
27+
N = 16 # size of the lattice, N x N
28+
eqSteps = 2000 # number of MC sweeps for equilibration
29+
mcSteps = 2000 # number of MC sweeps for calculation
30+
A_data, B_data, C_data, D_data = -15.261, 1.8389, 1.7137, 8.0904
31+
ps = 1.59
32+
T_low = 0.01
33+
T_high = 400
34+
35+
for line in data:
36+
key, value = line.split("=")
37+
dataset[key.strip()] = value.strip()
38+
# Read data
39+
nt = int(dataset["nt"])
40+
N = int(dataset["N"])
41+
eqSteps = int(dataset["eqSteps"])
42+
mcSteps = int(dataset["mcSteps"])
43+
A_data = float(dataset["A_data"])
44+
B_data = float(dataset["B_data"])
45+
C_data = float(dataset["C_data"])
46+
D_data = float(dataset["D_data"])
47+
ps = float(dataset["ps"])
48+
T_low = float(dataset["T_low"])
49+
T_high = float(dataset["T_high"])
50+
return nt, N, eqSteps, mcSteps, A_data, B_data, C_data, D_data, ps, T_low, T_high
51+
1852
def initialstate(N,ps):
1953
''' generates a random spin configuration for initial condition'''
2054
state = ps*np.random.randint(1,2, size=(N,N)) #1.51*(2*np.random.randint(2, size=(N,N))-1) #np.random.uniform(-2.98,2.98, size=(N,N))
@@ -31,7 +65,7 @@ def mcmove(config, beta, ps, A_data, B_data, C_data, D_data):
3165
nb = (config[(a+1)%N,b] + config[a,(b+1)%N] + config[(a-1)%N,b] + config[a,(b-1)%N])/4.
3266
# site polarization before change
3367
s = config[a, b]
34-
# site polarization after change
68+
# site polarization after change, note here I constrain the max polarization to be 2*ps.
3569
s1 = np.random.uniform(-2*ps,2*ps)
3670
# unitcell energy before and after
3771
cost_unitcell_orig = A_data/2.*s**2.+B_data/4.*s**4.+C_data/6.*s**6.
@@ -53,24 +87,22 @@ def calcMag(config):
5387
mag = np.abs(np.sum(config))
5488
return mag
5589

56-
## change these parameters for a smaller (faster) simulation
57-
nt = 18 # number of temperature points
58-
N = 30 # size of the lattice, N x N
59-
eqSteps = 8000 # number of MC sweeps for equilibration
60-
mcSteps = 4000 # number of MC sweeps for calculation
61-
A_data, B_data, C_data, D_data = -15.261, 1.8389, 1.7137, 8.0904
62-
ps = 1.59
63-
T = np.linspace(0.01, 400, nt);
64-
65-
E,M,C,X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
90+
#########################################################################
91+
# Here comes the Model parameters
92+
#########################################################################
93+
nt, N, eqSteps, mcSteps, A_data, B_data, C_data, D_data, ps, T_low, T_high = read_input()
94+
#########################################################################
95+
T = np.linspace(T_low, T_high, nt);
96+
E,M,C,X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
6697
E_1,M_1,C_1,X_1 = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
67-
n1, n2 = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N)
98+
n1, n2 = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N)
6899
# divide by number of samples, and by system size to get intensive values
69100

70101
#----------------------------------------------------------------------
71102
# MAIN PART OF THE CODE
72103
#----------------------------------------------------------------------
73104

105+
# automatically convert number of temperature steps to be devideable by "size"
74106
for tt in range(int(nt*rank/size),int(nt*(rank+1)/size)):
75107
E1 = M1 = E2 = M2 = 0
76108
config = initialstate(N, ps)
@@ -102,8 +134,7 @@ def calcMag(config):
102134
f = plt.figure(figsize=(18, 10)); # plot the calculated values
103135

104136
plt.scatter(T, abs(M), s=50, marker='o', color='RoyalBlue')
105-
plt.xlabel("Temperature (T)", fontsize=20);
137+
plt.xlabel("Temperature (K)", fontsize=20);
106138
plt.ylabel("Magnetization ", fontsize=20); plt.axis('tight');
107139

108140
plt.savefig("MC.png")
109-

‎MC_MPI_Ising.py

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
from __future__ import print_function
2+
from __future__ import division
3+
import numpy as np
4+
from numpy.random import rand
5+
import matplotlib
6+
matplotlib.use('Agg')
7+
import matplotlib.pyplot as plt
8+
9+
from mpi4py import MPI
10+
comm = MPI.COMM_WORLD
11+
rank = comm.Get_rank()
12+
size = comm.Get_size()
13+
14+
#----------------------------------------------------------------------
15+
## BLOCK OF FUNCTIONS USED IN THE MAIN CODE
16+
#----------------------------------------------------------------------
17+
def read_input():
18+
19+
''' a subroutine to get AO and MO number and number of kpoints and otehr information from input.woops'''
20+
21+
dataset={}
22+
file = open('input.MC', "r")
23+
#Default
24+
data = file.readlines()
25+
nt = 18 # number of temperature points
26+
N = 16 # size of the lattice, N x N
27+
eqSteps = 8000 # number of MC sweeps for equilibration
28+
mcSteps = 4000 # number of MC sweeps for calculation
29+
D_data = 1.
30+
T_low = 1.53
31+
T_high = 3.28
32+
33+
for line in data:
34+
key, value = line.split("=")
35+
dataset[key.strip()] = value.strip()
36+
# Read data
37+
nt = int(dataset["nt"])
38+
N = int(dataset["N"])
39+
eqSteps = int(dataset["eqSteps"])
40+
mcSteps = int(dataset["mcSteps"])
41+
D_data = float(dataset["D_data"])
42+
T_low = float(dataset["T_low"])
43+
T_high = float(dataset["T_high"])
44+
return nt, N, eqSteps, mcSteps, D_data, T_low, T_high
45+
46+
def initialstate(N):
47+
''' generates a random spin configuration for initial condition'''
48+
state = 2*np.random.randint(2, size=(N,N))-1
49+
return state
50+
51+
52+
def mcmove(config, beta,D_data):
53+
'''Monte Carlo move using Metropolis algorithm '''
54+
for i in range(N):
55+
for j in range(N):
56+
a = np.random.randint(0, N)
57+
b = np.random.randint(0, N)
58+
s = config[a, b]
59+
nb = config[(a+1)%N,b] + config[a,(b+1)%N] + config[(a-1)%N,b] + config[a,(b-1)%N]
60+
nnb = config[(a+1)%N,(b+1)%N] + config[(a+1)%N,(b-1)%N] + config[(a-1)%N,(b-1)%N] + config[(a-1)%N,(b+1)%N]
61+
cost = 2*s*D_data*nb#+2*0.5*s*nnb
62+
if cost < 0:
63+
s *= -1
64+
elif rand() < np.exp(-cost*beta):
65+
s *= -1
66+
config[a, b] = s
67+
return config
68+
69+
70+
def calcEnergy(config):
71+
'''Energy of a given configuration'''
72+
energy = 0
73+
for i in range(len(config)):
74+
for j in range(len(config)):
75+
S = config[i,j]
76+
nb = config[(i+1)%N, j] + config[i,(j+1)%N] + config[(i-1)%N, j] + config[i,(j-1)%N]
77+
energy += -nb*S
78+
return energy/4.
79+
80+
81+
def calcMag(config):
82+
'''Magnetization of a given configuration'''
83+
mag = np.abs(np.sum(config))
84+
return mag
85+
86+
#########################################################################
87+
# Here comes the Model parameters
88+
#########################################################################
89+
## change these parameters for a smaller (faster) simulation
90+
nt, N, eqSteps, mcSteps, D_data, T_low, T_high = read_input()
91+
#########################################################################
92+
T = np.linspace(T_low, T_high, nt);
93+
E,M,C,X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
94+
E_1,M_1,C_1,X_1 = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
95+
n1, n2 = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N)
96+
# divide by number of samples, and by system size to get intensive values
97+
98+
#----------------------------------------------------------------------
99+
# MAIN PART OF THE CODE
100+
#----------------------------------------------------------------------
101+
102+
for tt in range(int(nt*rank/size),int(nt*(rank+1)/size)):
103+
E1 = M1 = E2 = M2 = 0
104+
config = initialstate(N)
105+
iT=1.0/T[tt]; iT2=iT*iT;
106+
107+
for i in range(eqSteps): # equilibrate
108+
mcmove(config, iT, D_data) # Monte Carlo moves
109+
110+
for i in range(mcSteps):
111+
mcmove(config, iT, D_data)
112+
Ene = calcEnergy(config) # calculate the energy
113+
Mag = calcMag(config) # calculate the magnetisation
114+
115+
E1 = E1 + Ene
116+
M1 = M1 + Mag
117+
M2 = M2 + Mag*Mag
118+
E2 = E2 + Ene*Ene
119+
120+
E_1[tt] = n1*E1
121+
M_1[tt] = n1*M1
122+
C_1[tt] = (n1*E2 - n2*E1*E1)*iT2
123+
X_1[tt] = (n1*M2 - n2*M1*M1)*iT
124+
125+
comm.send(E_1,dest=0,tag=rank)
126+
comm.send(M_1,dest=0,tag=rank)
127+
comm.send(C_1,dest=0,tag=rank)
128+
comm.send(X_1,dest=0,tag=rank)
129+
#print(rank,E_1)
130+
131+
if rank == 0:
132+
for i in range(0,size):
133+
E += comm.recv(source=i,tag=i)
134+
M += comm.recv(source=i,tag=i)
135+
C += comm.recv(source=i,tag=i)
136+
X += comm.recv(source=i,tag=i)
137+
with open('Energy.txt', 'w') as f:
138+
for i in range(nt):
139+
print("{0:4d} {1:5f}".format(i,E[i]),file=f)
140+
with open('Polarization.txt', 'w') as f:
141+
for i in range(nt):
142+
print("{0:4d} {1:5f}".format(i,M[i]),file=f)
143+
with open('Specific_Heat.txt', 'w') as f:
144+
for i in range(nt):
145+
print("{0:4d} {1:5f}".format(i,C[i]),file=f)
146+
with open('Susceptibility.txt', 'w') as f:
147+
for i in range(nt):
148+
print("{0:4d} {1:5f}".format(i,X[i]),file=f)
149+
150+
f = plt.figure(figsize=(18, 10)); # plot the calculated values
151+
152+
sp = f.add_subplot(2, 2, 1 );
153+
plt.scatter(T, E, s=50, marker='o', color='IndianRed')
154+
plt.xlabel("Temperature (T)", fontsize=20);
155+
plt.ylabel("Energy ", fontsize=20); plt.axis('tight');
156+
157+
sp = f.add_subplot(2, 2, 2 );
158+
plt.scatter(T, abs(M), s=50, marker='o', color='RoyalBlue')
159+
plt.xlabel("Temperature (T)", fontsize=20);
160+
plt.ylabel("Magnetization ", fontsize=20); plt.axis('tight');
161+
162+
sp = f.add_subplot(2, 2, 3 );
163+
plt.scatter(T, C, s=50, marker='o', color='IndianRed')
164+
plt.xlabel("Temperature (T)", fontsize=20);
165+
plt.ylabel("Specific Heat ", fontsize=20); plt.axis('tight');
166+
167+
sp = f.add_subplot(2, 2, 4 );
168+
plt.scatter(T, X, s=50, marker='o', color='RoyalBlue')
169+
plt.xlabel("Temperature (T)", fontsize=20);
170+
plt.ylabel("Susceptibility", fontsize=20); plt.axis('tight');
171+
#plt.show()
172+
plt.savefig("MC.png")

0 commit comments

Comments
 (0)