Dear All

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)
  Inline image 1
The PRB.85.184412 result (please see fig 8_c of the paper):


Inline image 2

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.

W
hat 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