跳至内容
教程资源
蒙特卡洛模拟
MC-02 磁化率

MC-02 磁化率

磁化率 χ=∂⟨M⟩/∂h∣h=0\chi = \partial\langle M \rangle / \partial h \big|_{h=0} 衡量系统对外加小磁场的响应强度。 它的温度依赖关系是探测底层自旋关联的灵敏手段：弱相互作用自旋系统遵从居里定律 χ∝1/T\chi \propto 1/T，而相互作用和量子效应则会带来特征性的偏离。

本教程计算四个系统的 χ(T)\chi(T)——一维链和两腿梯子上的经典与量子海森堡模型——并把结果叠加在同一张图上。 这一比较突出了两个关键对照：量子涨落如何改变经典图像，以及把晶格几何从链改为梯子如何影响低温行为。

符号约定。 经典的 spinmc 程序和量子的 looper 程序都使用哈密顿量 H=J∑⟨i,j⟩S⃗i⋅S⃗jH = J \sum \langle i,j \rangle \vec{S}_i \cdot \vec{S}_j，其中 J<0J < 0 有利于铁磁排列，J>0J > 0 有利于反铁磁排列。下面的经典模拟使用 J=−1J = -1（铁磁体）；量子模拟使用 J=+1J = +1（反铁磁体）。这两种选择都给出在低温下增长的正磁化率，从而使经典行为与量子行为的比较更为直观。

经典海森堡模型

一维链

在命令行下设置并运行

参数文件 parm2a 设置了在一系列温度下、60 个格点的链上经典铁磁海森堡模型的模拟：

LATTICE="chain lattice"
L=60
J=-1
THERMALIZATION=15000
SWEEPS=500000
UPDATE="cluster"
MODEL="Heisenberg"
{T=0.05;}
{T=0.1;}
{T=0.2;}
{T=0.3;}
{T=0.4;}
{T=0.5;}
{T=0.6;}
{T=0.7;}
{T=0.8;}
{T=0.9;}
{T=1.0;}
{T=1.25;}
{T=1.5;}
{T=1.75;}
{T=2.0;}

用标准的命令序列运行模拟：

parameter2xml parm2a
spinmc --Tmin 10 --write-xml parm2a.in.xml

在 Python 中设置并运行

脚本 tutorial2a.py 设置并运行同样的模拟。把它放在与 parm2a 相同的目录下：

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

parms = []
for t in [0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0]:
    parms.append(
        {
            'LATTICE'        : "chain lattice",
            'T'              : t,
            'J'              : -1,
            'THERMALIZATION' : 10000,
            'SWEEPS'         : 500000,
            'UPDATE'         : "cluster",
            'MODEL'          : "Heisenberg",
            'L'              : 60
        }
    )

input_file = pyalps.writeInputFiles('parm2a', parms)
pyalps.runApplication('spinmc', input_file, Tmin=5, writexml=True)

计算与作图

从输出文件中载入磁化率，并画出它随温度的变化：

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

data = pyalps.loadMeasurements(pyalps.getResultFiles(prefix='parm2a'), 'Susceptibility')
susceptibility = pyalps.collectXY(data, x='T', y='Susceptibility')

plt.figure()
pyalps.plot.plot(susceptibility)
plt.xlabel('Temperature $T/J$')
plt.ylabel('Susceptibility $\chi J$')
plt.ylim(0, 0.22)
plt.title('Classical Heisenberg chain')
plt.show()

两腿梯子

梯子几何在腿上的耦合 J0J_0 之外，还引入了沿横档的第二个耦合 J1J_1。 除了晶格的改变和这两个耦合之外，模拟的设置与链的情形完全相同。

在命令行下设置并运行

下载 parm2b 并把它放在同一个目录下：

LATTICE="ladder"
L=60
J0=-1
J1=-1
THERMALIZATION=15000
SWEEPS=500000
UPDATE="cluster"
MODEL="Heisenberg"
{T=0.05;}
{T=0.1;}
{T=0.2;}
{T=0.3;}
{T=0.4;}
{T=0.5;}
{T=0.6;}
{T=0.7;}
{T=0.8;}
{T=0.9;}
{T=1.0;}
{T=1.25;}
{T=1.5;}
{T=1.75;}
{T=2.0;}

运行模拟：

parameter2xml parm2b
spinmc --Tmin 10 --write-xml parm2b.in.xml

在 Python 中设置并运行

脚本 tutorial2b.py 是 tutorial2a.py 的副本，只有三处改动：前缀改为 parm2b，LATTICE 设为 "ladder"，并把 J 换成 J0 和 J1（都取 -1）。

量子海森堡模型

量子模拟使用 ALPS 模型库来指定 S=1/2S = 1/2 的量子自旋模型，并用 looper 量子蒙特卡洛程序来运行模拟。 相对于经典情形，关键的参数改动是：

  • 用 MODEL="spin" 配合 local_S=1/2，而不是 MODEL="Heisenberg"
  • 用 ALGORITHM="loop" 来选择 looper 程序
  • 取 J=+1 作为反铁磁耦合（参见上面的符号约定说明）

一维链

在命令行下设置并运行

下载 parm2c：

LATTICE="chain lattice"
MODEL="spin"
local_S=1/2
L=60
J=1
THERMALIZATION=5000
SWEEPS=50000
ALGORITHM="loop"
{T=0.05;}
{T=0.1;}
{T=0.2;}
{T=0.3;}
{T=0.4;}
{T=0.5;}
{T=0.6;}
{T=0.7;}
{T=0.75;}
{T=0.8;}
{T=0.9;}
{T=1.0;}
{T=1.25;}
{T=1.5;}
{T=1.75;}
{T=2.0;}

用 loop 程序（而不是 spinmc）转换并运行：

parameter2xml parm2c
loop parm2c.in.xml

在 Python 中设置并运行

脚本 tutorial2c.py 把 tutorial2a.py 改写为量子参数，并调用 loop 而不是 spinmc：

input_file = pyalps.writeInputFiles('parm2c', parms)
pyalps.runApplication('loop', input_file)

两腿梯子

两腿梯子是准一维系统的一个简单例子：两条平行的海森堡链（腿），每条链上的耦合为 J0，它们由耦合为 J1 的横档键连接起来：

     J0      J0
  o-------o-------o--- ...     <- 腿（耦合 J0）
  |  J1   |  J1   |             <- 横档（耦合 J1）
  o-------o-------o--- ...     <- 腿（耦合 J0）

比值 α=J1/J0\alpha = J_1/J_0 控制两条腿的耦合强度，在解耦的链（α=0\alpha=0）与横档强耦合、近乎独立二聚体的系统（α≫1\alpha\gg 1）之间连续过渡。 与无能隙的链不同，两腿反铁磁海森堡梯子具有自旋能隙：它的基态是横档单重态的直积，χ\chi 在能隙以下被指数压低，在低温下陡峭地趋于零。

在命令行下设置并运行

下载 parm2d：

LATTICE="ladder"
MODEL="spin"
local_S=1/2
L=60
J0=1
J1=1
THERMALIZATION=5000
SWEEPS=50000
ALGORITHM="loop"
{T=0.1;}
{T=0.2;}
{T=0.3;}
{T=0.4;}
{T=0.5;}
{T=0.6;}
{T=0.7;}
{T=0.8;}
{T=1.0;}
{T=1.25;}
{T=1.5;}
{T=1.75;}
{T=2.0;}
parameter2xml parm2d
loop parm2d.in.xml

在 Python 中设置并运行

脚本 tutorial2d.py 是对 tutorial2c.py 的改写：把前缀改为 parm2d，把 LATTICE 改为 "ladder"，并把 J 换成 J0 和 J1（都取 1）。

汇总四次模拟

在同一个目录下运行完全部四次模拟之后，脚本 tutorial2full.py 会把所有结果一起载入，并叠加在同一张图上。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

data = pyalps.loadMeasurements(pyalps.getResultFiles(), 'Susceptibility')
data = pyalps.flatten(data)

不带参数调用 pyalps.getResultFiles() 会找出当前目录下所有的输出文件；pyalps.flatten 把来自不同文件的结果合并成一个列表，同时保留每个数据点对应的模拟参数。

我们使用 pyalps.collectXY 并配合 foreach 参数，为 MODEL 与 LATTICE 的每一种组合生成一条独立的曲线：

susceptibility = pyalps.collectXY(data, x='T', y='Susceptibility', foreach=['MODEL', 'LATTICE'])

利用每条曲线所保存的参数为它加上标签：

for s in susceptibility:
    if s.props['LATTICE'] == 'chain lattice':
        s.props['label'] = "chain"
    elif s.props['LATTICE'] == 'ladder':
        s.props['label'] = "ladder"
    if s.props['MODEL'] == 'spin':
        s.props['label'] = "quantum " + s.props['label']
    elif s.props['MODEL'] == 'Heisenberg':
        s.props['label'] = "classical " + s.props['label']

画出全部四条曲线：

plt.figure()
pyalps.plot.plot(susceptibility)
plt.xlabel('Temperature $T/J$')
plt.ylabel('Susceptibility $\chi J$')
plt.ylim(0, 0.25)
plt.legend()
plt.show()

所得的图应当与下图类似。 在高温下，四条曲线都趋于居里定律 χ∝1/T\chi \propto 1/T。 在低温下，量子梯子由于自旋能隙而陡然下落，而经典链、量子链以及经典梯子则保持有限值或继续增大。

思考题

  • 在高温下，四条曲线是否都趋于同样的 χ∝1/T\chi \propto 1/T 居里行为？是什么决定了前置因子？
  • 在低温下，经典链与量子链之间最明显的差别是什么？
  • 在低 TT 下，量子梯子的磁化率比经典梯子下降得快得多。是什么物理机制造成了这一点？
  • 试试更大的系统尺寸或不同的晶格（"cubic lattice"、"triangular lattice"——参见 lattices.xml）。结果会如何变化？
MC-01b 平衡化MC-03 磁化强度