diff --git a/docs/sphinx/en/tutorial/simple_tutorial-6.rst b/docs/sphinx/en/tutorial/simple_tutorial-6.rst index 26866c49..5b3f471a 100644 --- a/docs/sphinx/en/tutorial/simple_tutorial-6.rst +++ b/docs/sphinx/en/tutorial/simple_tutorial-6.rst @@ -1,7 +1,7 @@ .. highlight:: none Phase diagram of the hardcore boson model on a trianglar lattice ------------------------------------------------ +-------------------------------------------------------------------- Finally, let us consider a zero-temperature phase diagram of the hardcore boson model on a trianglar lattice. The Hamiltonian of this model is given as @@ -12,22 +12,28 @@ The Hamiltonian of this model is given as H = \sum_{\langle i,j \rangle} \Bigl[ -t (b_i^\dagger b_j + b_j^\dagger b_i) + V n_i n_j \Bigr] -\mu \sum_i n_i , \end{aligned} -where \ :math:`\langle i, j\rangle`\ indicates a pair of the nearest-neighbor sites, \ :math:`\mu`\ is a chemical potential, \ :math:`t`\ is a hopping energy, \ :math:`V`\ is a strength of the nearest-neighbor interaction. +where :math:`\langle i, j\rangle` indicates a pair of the nearest-neighbor sites, :math:`\mu` is a chemical potential, :math:`t` is a hopping energy, :math:`V` is a strength of the nearest-neighbor interaction. For a hardcore boson system, the maximum number of bosons at each site is restricted to 0 or 1. -It is known that several ordered phases characterized by two types of long-range order appear in this model :ref:`[Wessel] ` . -One is a solid-like order which exists at a 1/3 filling, where one of three sites is filled in a \ :math:`\sqrt{3}\times\sqrt{3}`\ ordering with wave vector \ :math:`\boldsymbol{Q}=(4\pi/3,0)`\ (see the inset of :numref:`fig_tutorial6_hardcore_boson` ). -This long-range order is characterized by the structure factor \ :math:`S(\boldsymbol{Q})`\. -The other is a superfluid order which is characterized by the offdiagonal order parameter \ :math:`|\langle b \rangle|`\. - -To perform calculation for this system, the user can use toml files named ``basic.toml`` , ``nn_obs.toml`` and a python script file ``run.py`` in the direction ``sample/06_hardcore_boson_triangular`` . Here, ``basic.toml`` specifies the model and its parameters. This file is almost the same as the triangular Heisenberg model described in the previous section and differs from it only in the section ``model`` in the last part, where the line ``type = "boson"`` specifies the hardcore boson model and ``t = 0.1`` , ``V = 1`` determines the strength of the hopping and nearest-neighbor interaction. - -The ``nn_obs.toml`` file describes the structure factor to be measured. - -.. literalinclude:: ../../../../sample/06_hardcore_boson_triangular/nn_obs.toml - -Using this file, the structure factor \ :math:`S(\boldsymbol{Q})`\ can be calculated. - -Let us execute calculation using the script ``run.py`` . +It is known that several ordered phases characterized by two types of long-range order appear in this model :ref:`[Wessel] `. +One is a superfluid order which is characterized by the offdiagonal order parameter :math:`|\langle b \rangle|`. +The other one is a solid-like order which exists at a 1/3 filling, where one of three sites is filled in a :math:`\sqrt{3}\times\sqrt{3}` ordering with wave vector :math:`\boldsymbol{Q}=(4\pi/3,0)` (see the inset of :numref:`fig_tutorial6_hardcore_boson`). +This long-range order is characterized by the structure factor :math:`S(\boldsymbol{Q}) = \sum_{ij}^{N_\text{sites}} \langle n_i n_j \rangle \exp[-i\boldsymbol{Q}\cdot(r_i-r_j)] /N_\text{sites}`. + +To perform calculation for this system, the user can use toml files named ``basic.toml``, ``nn_obs.toml`` and a python script file ``run.py`` in the direction ``sample/06_hardcore_boson_triangular``. +Here, ``basic.toml`` specifies the model and its parameters. +This file is almost the same as the triangular Heisenberg model described in the previous section and differs from it only in the section ``model`` in the last part, where the line ``type = "boson"`` specifies the hardcore boson model and ``t = 0.1``, ``V = 1`` determines the strength of the hopping and nearest-neighbor interaction. + +To calculate the structure factor :math:`S(\boldsymbol{Q})`, the correlations of densities at all the pairs of sites (in the unitcell) :math:`\langle n_i n_j \rangle` are required. +These observables are not defined by ``tenes_simple``, users should define them by themselves. +``nn_obs.toml`` specifies them for :math:`3 \times 3` unitcell, and the content of the file is appended to ``std__XXX_YYY.toml`` in ``run.py``. + +For larger unitcell, the computational cost increases, and hence the calculatation of the structure factor gets difficult. +On the other hand, TeNeS can calculate the correlation functions along the x and y directions (in the square lattice iTPS) at a low cost. +Structure factor can be calculated using these correlation functions. +Additionally, the Fourier transform of the density operator :math:`\langle n_i \rangle` is also available. +The ground state at 1/3 filling is three-fold degenerate, but one of them is selected in the finite bond dimension calculation, and calculated density has site dependence. + +Let us execute calculation using the script ``run.py``. After setting of the paths, execute calculation by typing the following command: :: @@ -41,18 +47,21 @@ After the calculation, start ``gnuplot`` and execute the following command: load 'plot.gp' -Then, we obtain a graph for the structure factor \ :math:`S(\boldsymbol{Q})`\ and the superfluid order parameter \ :math:`|\langle b \rangle|`\ as shown in :numref:`fig_tutorial6_hardcore_boson` (a). +Then, we obtain a graph as shown in :numref:`fig_tutorial6_hardcore_boson`. +:math:`S(\boldsymbol{Q})` is the structure factor calculated from the density correlations at all the pairs of sites in the unitcell, +:math:`S'(\boldsymbol{Q})` is the structure factor calculated from the density correlations along the x-axis, +:math:`n(\boldsymbol{Q})` is the Fourier transform of the densities, +and :math:`(|\langle b \rangle| + |\langle b^\dagger \rangle|)/2` is the superfluid order parameter. We note that this calculation is not so accurate because the bond dimension used in the calculation is small. -By commenting out the four lines in the beginning of the script ``run.py`` , we can perform more accurate calculation at the expense of execution time. -The result is shown in :numref:`fig_tutorial6_hardcore_boson` (b). -From these figures, we find that there exists three phases for the ground state, i.e., (a) a superfluid phase (\ :math:`-0.5 \lesssim \mu/V \lesssim -0.2`\), (b) a solid phase (\ :math:`-0.2 \lesssim \mu/V \lesssim 2.4`\), and (c) a supersolid phase (\ :math:`2.4 \lesssim \mu/V`\). +By increasing the bond dimensions specified in the beginning of the script ``run.py``, we can perform more accurate calculation at the expense of execution time. +From these figures, we find that there exists three phases for the ground state, i.e., (a) a superfluid phase (:math:`-0.5 \lesssim \mu/V \lesssim -0.2`), (b) a solid phase (:math:`-0.2 \lesssim \mu/V \lesssim 2.4`), and (c) a supersolid phase (:math:`2.4 \lesssim \mu/V`). This result is consistent with the previous work :ref:`[Wessel] ` . .. figure:: ../../img/tutorial_6_hardcore_boson.* :name: fig_tutorial6_hardcore_boson :width: 600px - Phase diagram of the hardcore boson model on a triangular lattice. (a) The bond dimension is 2. (b) The bond dimension is 5. Inset: a particle pattern for a solid phase. + Phase diagram of the hardcore boson model on a triangular lattice. Inset shows a particle pattern in a solid phase. .. rubric:: Reference diff --git a/docs/sphinx/img/tutorial_6_hardcore_boson.pdf b/docs/sphinx/img/tutorial_6_hardcore_boson.pdf index a00275e0..efaf5673 100644 Binary files a/docs/sphinx/img/tutorial_6_hardcore_boson.pdf and b/docs/sphinx/img/tutorial_6_hardcore_boson.pdf differ diff --git a/docs/sphinx/img/tutorial_6_hardcore_boson.png b/docs/sphinx/img/tutorial_6_hardcore_boson.png index 70846e69..1f3e1720 100644 Binary files a/docs/sphinx/img/tutorial_6_hardcore_boson.png and b/docs/sphinx/img/tutorial_6_hardcore_boson.png differ diff --git a/docs/sphinx/ja/tutorial/simple_tutorial-6.rst b/docs/sphinx/ja/tutorial/simple_tutorial-6.rst index 6134dc0d..6ea06174 100644 --- a/docs/sphinx/ja/tutorial/simple_tutorial-6.rst +++ b/docs/sphinx/ja/tutorial/simple_tutorial-6.rst @@ -15,16 +15,20 @@ ここで\ :math:`\langle i, j\rangle`\ は隣接サイトの組を表し、\ :math:`\mu`\ は化学ポテンシャル、\ :math:`t`\ はホッピングエネルギー、\ :math:`V`\ は隣接サイト間の相互作用の大きさを表します。 ハードコアボゾン模型では、各サイトのボソン数は0か1に制限されます。 この模型では、2つのタイプの長距離秩序によって特徴づけられるいくつかの秩序相が現れることが知られています :ref:`[Wessel] ` 。 -1/3フィリングでは、 :numref:`fig_tutorial6_hardcore_boson` の挿入図に示すような\ :math:`\sqrt{3}\times\sqrt{3}`\ 超格子構造をもつ固体相が出現します。これは波数\ :math:`\boldsymbol{Q}=(4\pi/3,0)`\の構造因子\ :math:`S(\boldsymbol{Q})`\を計算することで特徴づけることができます。 -もう一つの秩序は消滅演算子の期待値の大きさ\ :math:`|\langle b \rangle|`\で特徴づけられる超流動秩序です。 +まずは消滅演算子の期待値の大きさ\ :math:`|\langle b \rangle|`\で特徴づけられる超流動秩序です。 +もうひとつは固体秩序です。 1/3フィリングでは、 :numref:`fig_tutorial6_hardcore_boson` の挿入図に示すような\ :math:`\sqrt{3}\times\sqrt{3}`\ 超格子構造をもつ固体相が出現します。 +これは波数\ :math:`\boldsymbol{Q}=(4\pi/3,0)`\の構造因子\ :math:`S(\boldsymbol{Q}) = \sum_{ij}^{N_\text{sites}} \langle n_i n_j \rangle \exp[-i\boldsymbol{Q}\cdot(r_i - r_j)] / N_\text{site}`\ を計算することで特徴づけることができます。 -この模型の計算を行うには、 ``sample/06_hardcore_boson_triangular`` ディレクトリ中にある toml ファイル ``basic.toml`` , ``nn_obs.toml`` と、Pythonスクリプト ``run.py`` を利用します。 ``basic.toml`` ファイルには模型の設定やパラメータなどが記述されています。このファイルの記述は、前節の三角格子ハイゼンベルク模型とほぼ同じであるため、内容は割愛します。唯一の変更点は最後の ``model`` セクションだけです。ここで、 ``type = "boson"`` によってハードコアボゾン模型を指定しており、 ``t = 0.1`` , ``V = 1`` によってホッピングおよび隣接サイト間相互作用の大きさを指定しています。 +この模型の計算を行うには、 ``sample/06_hardcore_boson_triangular`` ディレクトリ中にある toml ファイル ``basic.toml``, ``nn_obs.toml`` と、Pythonスクリプト ``run.py`` を利用します。 ``basic.toml`` ファイルには模型の設定やパラメータなどが記述されています。このファイルの記述は、前節の三角格子ハイゼンベルク模型とほぼ同じであるため、内容は割愛します。唯一の変更点は最後の ``model`` セクションだけです。ここで、 ``type = "boson"`` によってハードコアボゾン模型を指定しており、 ``t = 0.1``, ``V = 1`` によってホッピングおよび隣接サイト間相互作用の大きさを指定しています。 -``nn_obs.toml`` ファイルには計算する構造因子の情報が記述されています。 +構造因子 :math:`S(\boldsymbol{Q})` を計算するためには、（ユニットセル中の）全サイト対における密度密度相関 :math:`\langle n_i n_j \rangle` を計算する必要があります。 +この演算子は ``tenes_simple`` では定義されないため、別途定義する必要があります。 +``nn_obs.toml`` ファイルに :math:`3 \times 3` ユニットセルにおける演算子が定義されており、 ``run.py`` で ``tenes_std`` の入力ファイル ``std_XXX_YYY.toml`` に追記しています。 -.. literalinclude:: ../../../../sample/06_hardcore_boson_triangular/nn_obs.toml - -このファイルによって、構造因子\ :math:`S(\boldsymbol{Q})`\を計算することができます。 +より大きなユニットセルを用いる場合には、密度密度相関を計算するコストが非常に大きくなってしまい、構造因子を計算することが現実的ではなくなります。 +一方、TeNeSは（正方格子iTPSの）x軸、y軸方向の相関関数を低コストで計算できます。これを利用して長距離秩序の有無を確認できます。 +また、1体の密度演算子の期待値 :math:`\langle n_i \rangle` のフーリエ変換 :math:`n(\boldsymbol{Q})` も使えます。 +1/3フィリングの基底状態は3重縮退していますが、有限ボンド次元計算ではそのうちのどれかが選ばれます。その結果、得られる密度はサイト依存性があります。 実際にスクリプト ``run.py`` を使って計算を実行してみましょう。 あらかじめ ``tenes`` などにパスを通した上で @@ -40,16 +44,17 @@ load 'plot.gp' -とすれば、 :numref:`fig_tutorial6_hardcore_boson` (a)のような構造因子\ :math:`S(\boldsymbol{Q})`\と超流動秩序パラメータ\ :math:`|\langle b \rangle|`\のグラフが得られます。 +とすれば、 :numref:`fig_tutorial6_hardcore_boson` のようなグラフが得られます。 +:math:`S(\boldsymbol{Q})` がユニットセル内の密度相関関数から計算した構造因子、 :math:`S'(\boldsymbol{Q})` がx方向の相関関数から計算した構造因子、 :math:`n(\boldsymbol{Q})` が密度のフーリエ変換、 :math:`(|\langle b \rangle| + |\langle b^\dagger \rangle|)/2` が超流動秩序パラメータです。 なお、ここで計算に用いたボンド次元は小さい(ボンド次元2)ため、計算精度はそれほど高くありません。 -スクリプト ``run.py`` の冒頭でコメントになっている4行をコメントアウトすることで、時間がかかりますが、より高精度の計算を行うこともできます。その結果を :numref:`fig_tutorial6_hardcore_boson` (b)に示します。 +スクリプト ``run.py`` の冒頭でボンド次元を増やすことで、時間がかかりますが、より高精度の計算を行うこともできます。 この図から、系の基底状態は3種類の相、つまり(a) 超流動相(\ :math:`-0.5 \lesssim \mu/V \lesssim -0.2`\), (b) 固体相(\ :math:`-0.2 \lesssim \mu/V \lesssim 2.4`\), および(c) 超流動固体相(\ :math:`2.4 \lesssim \mu/V`\)が現れることがわかります。これは先行研究の結果と整合します :ref:`[Wessel] ` 。 .. figure:: ../../img/tutorial_6_hardcore_boson.* :name: fig_tutorial6_hardcore_boson :width: 600px - 三角格子上のハードコアボゾン模型の基底状態相図. (a) ボンド次元が2のとき. (b) ボンド次元が5のとき. 挿入図: 固体相のときの粒子配置. + 三角格子上のハードコアボゾン模型の基底状態相図. .. rubric:: 参考文献 diff --git a/sample/06_hardcore_boson_triangular/basic.toml b/sample/06_hardcore_boson_triangular/basic.toml index 118ed06f..a0bf78e9 100644 --- a/sample/06_hardcore_boson_triangular/basic.toml +++ b/sample/06_hardcore_boson_triangular/basic.toml @@ -24,3 +24,7 @@ initial = "random" type = "boson" t = 0.1 V = 1.0 + +[correlation] +r_max = 30 +operators = [[0,0]] diff --git a/sample/06_hardcore_boson_triangular/nn_obs.toml b/sample/06_hardcore_boson_triangular/nn_obs.toml index 724d74d1..eb76cbf7 100644 --- a/sample/06_hardcore_boson_triangular/nn_obs.toml +++ b/sample/06_hardcore_boson_triangular/nn_obs.toml @@ -26,6 +26,54 @@ bonds = """ 2 0 2 2 1 2 2 2 2 +3 1 0 +3 2 0 +3 0 1 +3 1 1 +3 2 1 +3 0 2 +3 1 2 +3 2 2 +4 1 0 +4 2 0 +4 0 1 +4 1 1 +4 2 1 +4 0 2 +4 1 2 +4 2 2 +5 1 0 +5 2 0 +5 0 1 +5 1 1 +5 2 1 +5 0 2 +5 1 2 +5 2 2 +6 1 0 +6 2 0 +6 0 1 +6 1 1 +6 2 1 +6 0 2 +6 1 2 +6 2 2 +7 1 0 +7 2 0 +7 0 1 +7 1 1 +7 2 1 +7 0 2 +7 1 2 +7 2 2 +8 1 0 +8 2 0 +8 0 1 +8 1 1 +8 2 1 +8 0 2 +8 1 2 +8 2 2 """ ops = [0, 0] diff --git a/sample/06_hardcore_boson_triangular/plot.gp b/sample/06_hardcore_boson_triangular/plot.gp index 825629f2..ff9a3ebc 100644 --- a/sample/06_hardcore_boson_triangular/plot.gp +++ b/sample/06_hardcore_boson_triangular/plot.gp @@ -1,9 +1,13 @@ set key left top -set xlabel 'mu/V' -set ylabel 'S(Q)/N' -set y2label '|b|' +# set xlabel "mu/V" +# set ylabel "S(Q), N(Q)" +# set y2label "|b|+|b'|" +set ytics nomirror +set y2tics nomirror set style data lp plot \ - 'sq.dat' u 1:3 pt 4 t"S(Q)",\ - 'offdiag.dat' u 1:3 axes x1y2 pt 5 t"|b|" + 'nq.dat' u 1:(abs($3)) pt 4 t"|N(Q)|",\ + 'sq.dat' u 1:(abs($3)) pt 5 t"|S(Q)|",\ + 'sq_x.dat' u 1:(abs($3)) pt 6 t"|S'(Q)|",\ + 'offdiag.dat' u 1:(abs($3)) axes x1y2 pt 7 t"|b|+|b'|" diff --git a/sample/06_hardcore_boson_triangular/run.py b/sample/06_hardcore_boson_triangular/run.py index bdf81bfa..d9ff6920 100644 --- a/sample/06_hardcore_boson_triangular/run.py +++ b/sample/06_hardcore_boson_triangular/run.py @@ -1,31 +1,53 @@ +from __future__ import annotations import subprocess +from typing import TextIO from os.path import join import numpy as np import toml MPI_cmd = "" # mpiexec -np 4" +# calculate S(q) = \sum exp[-i q(r_i - r_j)] or not +# If True, only 3x3 unitcell is available +calculate_sq = True + +# if False, skip computing and only gather results into *.dat (e.g., nq.dat) do_calculation = True -decrease_mu = True # sweep mu from larger to smaller or not -num_h = 21 -min_h = -1.0 -max_h = 3.0 +# sweep mu from larger to smaller or not +decrease_mu = True + +num_mu = 21 +min_mu = -1.0 +max_mu = 3.0 num_step_table = [1000, 2000] D = 2 chi = 4 + # Please uncomment the following lines if you hope to obtain more accurate results -# num_h = 41 +# num_mu = 41 # num_step_table = [1500, 3000] # D = 5 # chi = 25 -fsq = open("sq.dat", "w") -fmag0 = open("offdiag.dat", "w") -fmag = open("density.dat", "w") -fene = open("energy.dat", "w") - -for f in (fsq, fmag, fmag0, fene): +# file objects for physical quantities +files: dict[str, TextIO] = {} + +# N(q) = (1/Nsite) \sum_i exp[-i q r_i] with q = (4pi/3, 0) +files["nq"] = open("nq.dat", "w") +if calculate_sq: + # S(q) = (1/Nsite) \sum_{ij} exp[-i q (r_i-r_j)] with q = (4pi/3, 0) + files["sq"] = open("sq.dat", "w") +# S(q)' = (1/Nsite) \sum_i \sum_dx exp[-i q dx] with q = (4pi/3, 0) +files["sq_x"] = open("sq_x.dat", "w") +# N0 = ( + ) / 2 +files["den0"] = open("offdiag.dat", "w") +# N = (1/Nsite) \sum_i +files["den"] = open("density.dat", "w") +# E = +files["ene"] = open("energy.dat", "w") + +for f in files.values(): f.write("# $1: mu\n") for i, num_step in enumerate(num_step_table, 2): f.write(f"# ${i}: num_step={num_step}\n") @@ -34,14 +56,12 @@ step_mu = -1 if decrease_mu else 1 -for idx, h in enumerate(np.linspace(min_h, max_h, num=num_h)[::step_mu]): - print(f"Calculation Process: {idx+1}/{num_h}") +for idx, mu in enumerate(np.linspace(min_mu, max_mu, num=num_mu)[::step_mu]): + print(f"Calculation Process: {idx+1}/{num_mu}") inum = 0 num_pre = 0 - fsq.write("{} ".format(h)) - fmag0.write("{} ".format(h)) - fmag.write("{} ".format(h)) - fene.write("{} ".format(h)) + for f in files.values(): + f.write("{} ".format(mu)) for num_step in num_step_table: ns = num_step - num_pre print(f"Steps: {num_step}") @@ -49,12 +69,21 @@ if do_calculation: with open("basic.toml") as f: dict_toml = toml.load(f) + + L = dict_toml["lattice"]["L"] + W = dict_toml["lattice"]["W"] + nsites = L * W + if calculate_sq and nsites != 9: + raise RuntimeError("The number of sites is not 9 (3x3)!") + dict_toml["parameter"]["general"]["output"] = output_dir dict_toml["parameter"]["general"]["tensor_save"] = "tensor_save" - dict_toml["model"]["mu"] = float(h) + dict_toml["model"]["mu"] = float(mu) dict_toml["lattice"]["virtual_dim"] = D dict_toml["parameter"]["ctm"]["dimension"] = chi dict_toml["parameter"]["simple_update"]["num_step"] = ns + + # load the previous optimized tensor if idx > 0 or inum > 0: dict_toml["parameter"]["general"]["tensor_load"] = "tensor_save" @@ -67,73 +96,120 @@ cmd = f"tenes_simple {simple_toml} -o {std_toml}" subprocess.call(cmd.split()) - with open(std_toml, "a") as fout: - with open("nn_obs.toml") as fin: - for line in fin: - fout.write(line) + if calculate_sq: + # append twosite observable to std.toml + with open(std_toml, "a") as f: + with open("nn_obs.toml") as fin: + for line in fin: + f.write(line) cmd = f"tenes_std {std_toml} -o {input_toml}" subprocess.call(cmd.split()) cmd = f"{MPI_cmd} tenes {input_toml}" subprocess.call(cmd.split()) - ene = "" - density = "" + with open(join(output_dir, "density.dat")) as f: for line in f: + if line.startswith("#"): + continue words = line.split("=") if words[0].strip() == "N": - density = words[1].strip() + density = float(words[1].strip().split()[0]) + files["den"].write(f"{density} ") if words[0].strip() == "Energy": - ene = words[1].strip() - fene.write(f"{ene} ") - fmag.write(f"{density} ") - fene.flush() - fmag.flush() - - n0 = 0.0 - N = np.zeros(9) + ene = float(words[1].strip().split()[0]) + files["ene"].write(f"{ene} ") + + coords_list = [] + with open("coordinates.dat") as f: + for line in f: + words = line.split() + if len(words) == 0: + continue + if words[0].startswith("#"): + continue + x = float(words[1]) + y = float(words[2]) + coords_list.append([x, y]) + coords = np.array(coords_list) + + nsites = coords.shape[0] + q = np.array([4 * np.pi / 3, 0]) + nq_coeffs = np.cos(coords @ q) + + n0 = 0.0 # ( + ) / 2 + nq = 0.0 # n(q), q = (4pi/a, 0) + N = np.zeros(nsites) with open(join(output_dir, "onesite_obs.dat")) as f: for line in f: + if line.startswith("#"): + continue words = line.split() if len(words) == 0: continue - if words[0].strip() == "0": + if words[0].strip() == "0": # N + i = int(words[1]) n = float(words[2]) - N[int(words[1])] = n - if words[0].strip() == "1": + nq += nq_coeffs[i] * n + N[i] = n + if words[0].strip() == "1": # Bdagger n = float(words[2]) - n0 += np.abs(n) / 18.0 - if words[0].strip() == "2": + n0 += 0.5 * np.abs(n) + if words[0].strip() == "2": # B n = float(words[2]) - n0 += np.abs(n) / 18.0 - fmag0.write("{} ".format(n0)) - fmag0.flush() + n0 += 0.5 * np.abs(n) + files["den0"].write(f"{n0 / nsites} ") + files["nq"].write(f"{nq / nsites} ") + + if calculate_sq: + sq = 0.0 # S(q), q = (4pi/a, 0) + for i in range(nsites): + sq += N[i] # \hat{n}_i^2 = \hat{n}_i for hardcore boson + with open(join(output_dir, "twosite_obs.dat")) as f: + for line in f: + words = line.split() + if len(words) == 0: + continue + if words[0].strip() == "4": # + dx = int(words[2]) + dy = int(words[3]) + nn = float(words[4]) + if dx == dy: + sq += nn + else: + sq -= 0.5 * nn + files["sq"].write(f"{sq / nsites} ") sq = 0.0 - for i in range(3): + for i in range(nsites): + # C(r=0) + # \hat{n}_i^2 = \hat{n}_i for hardcore boson sq += N[i] - - with open(join(output_dir, "twosite_obs.dat")) as f: + with open(join(output_dir, "correlation.dat")) as f: for line in f: + if line.startswith("#"): + continue words = line.split() if len(words) == 0: continue - if words[0].strip() == "4": - n = float(words[4]) - if words[2] != words[3]: - n *= -0.5 - sq += n - sq *= 3.0 / 9 / 9 - fsq.write(f"{sq} ") - fsq.flush() + dx = int(words[3]) + dy = int(words[4]) + nn = float(words[5]) + + # Consider only the C(r) along x-axis + if dy > 0: + continue + if dx % 3 == 0: + sq += nn + else: + sq += -0.5 * nn + files["sq_x"].write(f"{sq / nsites} ") inum = inum + 1 num_pre = num_step - fene.write("\n") - fmag.write("\n") - fmag0.write("\n") - fsq.write("\n") -fene.close() -fmag.close() -fmag0.close() -fsq.close() + for f in files.values(): + f.flush() + for f in files.values(): + f.write("\n") +for f in files.values(): + f.close()