Recently I tried to reproduce the result of a paper about two honeycomb bilayer (PRB.85.184412) .
The paper used simple MC simulations ( in fact authors used exchange MC technique) with classical Heisenberg model to derive
susceptibility of this system. Authors used three different exchange parameters:
1-J1 as NN interaction in honeycomb surfaces
2-J2 as NNN interactions in honeycomb surfaces
3-Jc as face-to face intra-bilayer interactions.
I used exactly the same parameters of the paper (I use minus sign for exchange). My setup is as follows:
******************************************
******************************************
Exchange parameters setting:
J1=30.2 K, J2/J1=0.1 and Jc/J1=0.1.
Landee g-factor:
g=2
steps for thermalization:
50000
steps for sampling:
1000000
lattice size:
L=24
******************************************
******************************************
All of my parameters are the same as the paper except that I used more steps for sampling and for thermalization because the system is frustrated.
The following plots are mine and from the paper ( I derived the data of fig(8,c) of the paper by using G3DATA):
The ALPS result (mine)
The PRB.85.184412 result (please see fig 8_c of the paper):
As you can see the plot of mine has a maximum around 40 K while the plot of PRB.85.184412 has a maximum around 60-70 K and also the plot shows a minimum around 20 K.
What do you think? what is my fault? Is it a bug or something else?
If somebody like to repeat the calculation, these are my input files:
the python script:
#=================================
import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
#prepare the input parameters
parms = []
for t in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 43, 46, 49, 52, 55, 58, 61, 64, 67, 70, 75, 80, 85, 90, 95, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400 ]:
parms.append(
{
'LATTICE_LIBRARY' : "graphene_bilayer.xml",
'LATTICE' : "graphene bilayer lattice",
'S' : 1.5,
'T' : t,
'J1' : -30.2,
'J2' : -3.02,
'J3' : -3.02,
'g' : 2.00,
'THERMALIZATION' : 50000,
'SWEEPS' : 1000000,
'UPDATE' : "local",
'MODEL' : "Heisenberg",
'L' : 24
}
)
#write the input file and run the simulation
input_file = pyalps.writeInputFiles('parm2b',parms)
pyalps.runApplication('spinmc',input_file,Tmin=300, MPI=4, mpirun='mpirun')
#get the list of result files
result_files = pyalps.getResultFiles(prefix='parm2b')
print "Loading results from the files: ", result_files
#print the observables stored in those files:
print "The files contain the following mesurements:",
print pyalps.loadObservableList(result_files)
#load a selection of measurements:
data = pyalps.loadMeasurements(result_files,['Susceptibility'])
#make a plot for the magnetization: collect Magnetziation as function of T
plotdata = pyalps.collectXY(data,'T','Susceptibility')
# convert the data to text file for plotting using another tool
print "The results in plain text are:"
print pyalps.plot.convertToText(plotdata)
print "Saving into file parm2b.txt"
f = open ('parm2b.txt','w')
f.write(pyalps.plot.convertToText(plotdata))
f.close()
#================================
the lattice file (xml format)
=====================================
<LATTICES>
<LATTICE name="hexagonal lattice" dimension="3">
<PARAMETER name="a" default="1"/>
<PARAMETER name="c" default="1"/>
<BASIS>
<VECTOR>a 0 0</VECTOR>
<VECTOR>a/2 a*sqrt(3)/2 0</VECTOR>
<VECTOR>0 0 c</VECTOR>
</BASIS>
<RECIPROCALBASIS>
<VECTOR>2*pi/a -2*pi/a/sqrt(3) 0 </VECTOR>
<VECTOR>0 4*pi/a/sqrt(3) 0 </VECTOR>
<VECTOR>0 0 2*pi/c</VECTOR>
</RECIPROCALBASIS>
</LATTICE>
<UNITCELL name="graphene bilayer" dimension="3" vertices="4">
<VERTEX><COORDINATE> 0.000000000000 0.0000000000 0 </COORDINATE></VERTEX>
<VERTEX><COORDINATE> 0.333333333333 0.3333333333 0</COORDINATE></VERTEX>
<VERTEX><COORDINATE> 0.000000000000 0.0000000000 0.5 </COORDINATE></VERTEX>
<VERTEX><COORDINATE> 0.333333333333 0.3333333333 0.5</COORDINATE></VERTEX>
<EDGE type="1"><SOURCE vertex="1"/><TARGET vertex="2" offset=" 0 0 0"/></EDGE>
<EDGE type="1"><SOURCE vertex="3"/><TARGET vertex="4" offset=" 0 0 0"/></EDGE>
<EDGE type="1"><SOURCE vertex="2"/><TARGET vertex="1" offset=" 1 0 0"/></EDGE>
<EDGE type="1"><SOURCE vertex="2"/><TARGET vertex="1" offset=" 0 1 0"/></EDGE>
<EDGE type="1"><SOURCE vertex="4"/><TARGET vertex="3" offset=" 1 0 0"/></EDGE>
<EDGE type="1"><SOURCE vertex="4"/><TARGET vertex="3" offset=" 0 1 0"/></EDGE>
<EDGE type="2"><SOURCE vertex="1"/><TARGET vertex="3" offset=" 0 0 0"/></EDGE>
<EDGE type="2"><SOURCE vertex="2"/><TARGET vertex="4" offset=" 0 0 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="1"/><TARGET vertex="1" offset=" 1 0 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="1"/><TARGET vertex="1" offset=" 0 1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="1"/><TARGET vertex="1" offset=" 1 -1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="2"/><TARGET vertex="2" offset=" 1 0 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="2"/><TARGET vertex="2" offset=" 0 1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="2"/><TARGET vertex="2" offset=" 1 -1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="3"/><TARGET vertex="3" offset=" 1 0 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="3"/><TARGET vertex="3" offset=" 0 1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="3"/><TARGET vertex="3" offset=" 1 -1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="4"/><TARGET vertex="4" offset=" 1 0 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="4"/><TARGET vertex="4" offset=" 0 1 0"/></EDGE>
<EDGE type="3"><SOURCE vertex="4"/><TARGET vertex="4" offset=" 1 -1 0"/></EDGE>
</UNITCELL>
<LATTICEGRAPH name = "graphene bilayer lattice">
<FINITELATTICE>
<LATTICE ref="hexagonal lattice"/>
<PARAMETER name="W" default="L"/>
<PARAMETER name="H" default="1"/>
<EXTENT dimension="1" size="L"/>
<EXTENT dimension="2" size="W"/>
<EXTENT dimension="3" size="H"/>
<BOUNDARY type="periodic"/>
</FINITELATTICE>
<UNITCELL ref="graphene bilayer"/>
</LATTICEGRAPH>
</LATTICES>
=======================================================================
Thanks in advanced
Mojtaba Alaei