diff --git a/doc/en/chap05_en.tex b/doc/en/chap05_en.tex index 59ef1b66..fb207921 100644 --- a/doc/en/chap05_en.tex +++ b/doc/en/chap05_en.tex @@ -2,77 +2,92 @@ %---------------------------------------------------------- \chapter{{Algorithm}} \label{Ch:algorithm} -\section{Variational MonteCalro Method} -In Variational MonteCalro (VMC) method, the importance sampling is performed in the system constructed by the Malkov chain toward the appropriate complete system. -Here, we choose the real spatial arrangement $\{| x\rangle\}$ at $S_z = 0$ as the complete system. +\section{Variational Monte Calro Method} + +The variational Monte Carlo (VMC) method is a method for calculating approximate wave functions +of a ground state and low-lying excited states by optimizing variational parameters included in a trial +wave function. +In calculating expectation values of physical quantities for the trial wave functions, +the Markov chain Monte Carlo method is applied for efficient important sampling. + +In the mVMC package, we choose a spatial configuration for electrons as a complete set of bases in sampling: \begin{equation} | x\rangle = \prod_{n=1}^{N/2} c_{r_{n\uparrow}}^{\dag} \prod_{n=1}^{N/2} c_{r_{n\downarrow}}^{\dag} |0 \rangle, \end{equation} -where we set $r_{n\sigma}$ as the position of $n$-th electron with $\sigma (=\uparrow \rm{or} \downarrow)$ spin. - -\subsection{Importance sampling} -When the importance of Malkov chain is defined by +where $r_{n\sigma}$ is a position of $n$-th electron with $\sigma (=\uparrow \rm{or} \downarrow)$ spin, +and $c_{r_{n\sigma}}^{\dag}$ is a creation operator of electrons. +By using this basis set, the expectation value of an operator $A$ is expressed as +\begin{equation} +\langle A \rangle =\frac{\langle \psi| A| \psi \rangle}{\langle \psi | \psi \rangle} +=\sum_x \frac{\langle \psi| A | x\rangle \langle x| \psi \rangle}{\langle \psi |\psi \rangle}. +\end{equation} +If we define a weight of the Markov chain Monte Carlo method as \begin{equation} -\rho(x)=\frac{|\langle x| \psi \rangle|^2}{\langle \psi | \psi \rangle} \ge 0, \sum{x} \rho(x)=1, +\rho(x)=\frac{|\langle x| \psi \rangle|^2}{\langle \psi | \psi \rangle} \ge 0, \quad \sum_{x} \rho(x)=1, \end{equation} -the expected value of the operator $A$ is given by +we can rewrite $\langle A \rangle$ in the following form: \begin{equation} -\langle A \rangle =\frac{\langle \psi| A| \psi \rangle}{\langle \psi | \psi \rangle} -=\sum_x \frac{\langle \psi| A | x\rangle \langle x| \psi \rangle}{\langle \psi |\psi \rangle} -=\sum_x \rho(x) \frac{\langle \psi| A | x\rangle }{\langle \psi |x \rangle} . +\langle A \rangle =\sum_x \rho(x) \frac{\langle \psi| A | x\rangle }{\langle \psi |x \rangle}. \end{equation} -In VMC, the summation in terms of $x$ is replaced by the importance sampling. -Then, local Green's function $G_{ij\sigma\sigma'}(x)$ is defined by +By using this form, the Markov chain Monte Carlo method is performed for sampling +with respect to $x$. The local Green's function $G_{ij\sigma\sigma'}(x)$, which is defined as \begin{equation} -G_{ij\sigma\sigma'}(x)=\frac{\langle \psi | c_{i\sigma}^{\dag} c_{j\sigma'} | \psi \rangle}{\langle \psi | x \rangle}. +G_{ij\sigma\sigma'}(x)=\frac{\langle \psi | c_{i\sigma}^{\dag} c_{j\sigma'} | \psi \rangle}{\langle \psi | x \rangle}, \end{equation} -In mVMC, the Mersenne twister method is used for sampling as the random number generator \cite{Mutsuo2008}. +is also evaluated by the same sampling method by taking $A = c_{i\sigma}^{\dag} c_{j\sigma'}$. +We adopt the Mersenne twister method as a random number generator for sampling~\cite{Mutsuo2008}. \section{Bogoliubov representation}\label{sec_bogoliubov_rep} -In the spin system, -the spin indices in input files of \verb|transfer|, \verb|InterAll|, -and correlation functions are specified as those of the Bogoliubov representation. -Spin operators are written by using creation/annihilation operators as follows: +In the VMC calculation for spin systems, we use the Bogoliubov representation. +In the input files defining the one-body term (\verb|transfer|) and the two-body term (\verb|InterAll|), +and the output files for correlation functions, the indices must be assigned by the Bogoliubov representation, +in which the spin operators are generally expressed by creation/annihilation operators of fermions as \begin{align} - S_{i z} &= \sum_{\sigma = -S}^{S} \sigma c_{i \sigma}^\dagger c_{i \sigma} + S_{i z} &= \sum_{\sigma = -S}^{S} \sigma c_{i \sigma}^\dagger c_{i \sigma}, \\ S_{i}^+ &= \sum_{\sigma = -S}^{S-1} \sqrt{S(S+1) - \sigma(\sigma+1)} - c_{i \sigma+1}^\dagger c_{i \sigma} + c_{i \sigma+1}^\dagger c_{i \sigma}, \\ S_{i}^- &= \sum_{\sigma = -S}^{S-1} \sqrt{S(S+1) - \sigma(\sigma+1)} - c_{i \sigma}^\dagger c_{i \sigma+1} + c_{i \sigma}^\dagger c_{i \sigma+1}. \end{align} +Since the present package support only $S=1/2$ spin systems, the Bogoliubov representation obtained +by substituting $S=1/2$ into the above equations is used. -\section{Relation between Pfaffian Slater determinant and single Slater determinant} +\section{Properties of the Pfaffian-Slater determinant} \label{sec:PuffAndSlater} -In this section, we show relation between Pfaffian Slater determinant and single Slater determinant. -We also discuss meaning of the singular value decomposition of coefficients $f_{ij}$. -\subsection{Relation between $f_{ij}$ and $\Phi_{in\sigma}$~(case of anti-parallel pairing)} -Pfaffian Slater determinant [one-body part of the many-variable variational Monte Carlo (mVMC) method] -is defined as + +In this section, we explain some properties of the Pfaffian-Slater determinant. +We derive the general relation between a Pfaffian-Slater determinant and a single Slater determinant +in Sec.~\ref{sec:PfaffianAP} and~\ref{sec:PfaffianP}. +We also discuss meaning of the singular value decomposition of coefficients $f_{ij}$ +in Sec.~\ref{sec:PfaffianSingular}. + +\subsection{Relation between $f_{ij}$ and $\Phi_{in\sigma}$~(the case of the anti-parallel pairing)} +\label{sec:PfaffianAP} + +In the many-variable variational Monte Carlo (mVMC) method, +the one-body part of the trial wave function is expressed by the Pfaffian Slater determinant defined as \begin{equation} |\phi_{\rm Pf}\rangle=\Big(\sum_{i,j=1}^{N_{s}}f_{ij}c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \end{equation} -where $N_{s}$ is number of sites, -$N_{e}$ is number of total particles, -and $f_{ij}$ are variational parameters. -For simplicity, we assume that $f_{ij}$ are real number. -Single Slater determinant is defined as +where $N_{s}$ is number of sites, $N_{e}$ is number of total particles, and $f_{ij}$ are variational parameters. +For simplicity, we assume that $f_{ij}$ are a real number. +The single Slater determinant is defined as \begin{align} |\phi_{\rm SL}\rangle&=\Big(\prod_{n=1}^{N_{e}/2}\psi_{n\uparrow}^{\dagger}\Big) \Big(\prod_{m=1}^{N_{e}/2}\psi_{m\downarrow}^{\dagger}\Big)|0\rangle, \\ -\psi_{n\sigma}^{\dagger}&=\sum_{i=1}^{N_{s}}\Phi_{in\sigma}c^{\dagger}_{i\sigma}. +\psi_{n\sigma}^{\dagger}&=\sum_{i=1}^{N_{s}}\Phi_{in\sigma}c^{\dagger}_{i\sigma}, \end{align} -We note that $\Phi$ is the normalized orthogonal basis, i.e, +Here, $\Phi_{in\sigma}$ is an orthonormal basis, i.e., satisfies \begin{equation} \sum_{i=1}^{N_{s}}\Phi_{in\sigma}\Phi_{im\sigma}=\delta_{nm}, \end{equation} where $\delta_{nm}$ is the Kronecker's delta. -Due to this normalized orthogonality, we obtain -following relation: +From this orthogonality, we can prove the relation \begin{align} [\psi^{\dagger}_{n\sigma},\psi_{m\sigma}]_{+}&=\delta_{nm},\\ G_{ij\sigma}=\langle c_{i\sigma}^{\dagger}c_{j\sigma}\rangle @@ -80,81 +95,76 @@ \subsection{Relation between $f_{ij}$ and $\Phi_{in\sigma}$~(case of anti-parall &=\sum_{n} \Phi_{in\sigma} \Phi_{jn\sigma}. \end{align} -Here, we rewrite $\phi_{\rm SL}$ and obtain explicit -relation between $f_{ij}$ and $\Phi_{in\sigma}$. -By using the commutation relation for $\psi^{\dagger}_{n\sigma}$, -we rewrite $\phi_{\rm SL}$ as +Next, let us prove the relation between $f_{ij}$ and $\Phi_{in\sigma}$ by modifying $|\phi_{\rm SL}\rangle $. +By the commutation relation for $\psi^{\dagger}_{n\sigma}$, $|\phi_{\rm SL}\rangle$ is rewritten as \begin{align} |\phi_{\rm SL}\rangle \propto \prod_{n=1}^{N_{e}/2}\Big(\psi_{n\uparrow}^{\dagger}\psi_{\mu(n)\downarrow}^{\dagger}\Big)|0\rangle, \end{align} -where $\mu(n)$ represents permutation of sequence of $n= 1, 2, \cdots, N_{e}/2$. -For simplicity, we take identity permutation and obtain the relation +where $\mu(n)$ represents permutation of a sequence of natural numbers, $n= 1, 2, \cdots, N_{e}/2$. +For simplicity, let us take identity permutation ($\mu(n) = n$). +By defining $K_{n}^{\dagger}=\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}$, and by +using the relation $K_{n}^{\dagger}K_{m}^{\dagger}=K_{m}^{\dagger}K_{n}^{\dagger}$, +we can derive the relation \begin{align} |\phi_{\rm SL}\rangle &\propto \prod_{n=1}^{N_{e}/2}\Big(\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}\Big)|0\rangle =\prod_{n=1}^{N_{e}/2} K_{n}^{\dagger}|0\rangle \\ &\propto\Big(\sum_{n=1}^{\frac{N_{e}}{2}}K_{n}^{\dagger}\Big)^{\frac{N_{e}}{2}} |0\rangle =\Big(\sum_{i,j=1}^{N_{s}}\Big[\sum_{n=1}^{\frac{N_{e}}{2}}\Phi_{in\uparrow}\Phi_{jn\downarrow}\Big] -c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)|0\rangle, +c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)|0\rangle. \end{align} -where $K_{n}^{\dagger}=\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}$ and -we use the relation $K_{n}^{\dagger}K_{m}^{\dagger}=K_{m}^{\dagger}K_{n}^{\dagger}$. -This result shows that $f_{ij}$ can be expressed by the -coefficients of the single Slater determinant as +This result indicates that $f_{ij}$ is expressed by the coefficients of the single Slater determinant as \begin{align} f_{ij}=\sum_{n=1}^{\frac{N_{e}}{2}}\Phi_{in\uparrow}\Phi_{jn\downarrow}. \end{align} -We note that this is one of expression of $f_{ij}$ for -single Slater determinant, i.e, $f_{ij}$ depend on -the pairing degrees of freedom (choices of $\mu(n)$) and -gauge degrees of freedom ( we can arbitrary change -the sign of $\Phi$ as $\Phi_{in\sigma}\rightarrow -\Phi_{in\sigma}$). -This large degrees of freedom is the origin of huge redundancy of -$f_{ij}$. +We note that this is one of a number of possible expressions of $f_{ij}$ derived from one single Slater determinant. +Since $f_{ij}$ depends not only on the choice of the pairing degrees of freedom (i.e., the choice of $\mu(n)$) but also on +the choice of the gauge degrees of freedom (i.e., the sign of $\Phi_{in\sigma}$), +the parameter $f_{ij}$ has huge redundancy. -\subsection{Relation between $F_{IJ}$ and $\Phi_{In}$~(case of general pairing including parallel pairing)} -Here, we consider the relation between the general Pfaffian -wave functions that include the parallel pairing and -the Slater determinant. +\subsection{Relation between $F_{IJ}$ and $\Phi_{In\sigma}$~(the case of the general pairing)} +\label{sec:PfaffianP} + +We extend the relation between the Pfaffian-Slater wave function and the single Slater wave function +into the general pairing case including the spin-parallel pairing. +We define the Pfaffian-Slater wave function and the single Slater wave function as \begin{align} -|\phi_{\rm Pf}\rangle&=\Big(\sum_{I,J=1}^{2N_{s}}f_{IJ}c_{I}^{\dagger}c_{J}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \\ +|\phi_{\rm Pf}\rangle&=\Big(\sum_{I,J=1}^{2N_{s}}F_{IJ}c_{I}^{\dagger}c_{J}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \\ |\phi_{\rm SL}\rangle&=\Big(\prod_{n=1}^{N_{e}}\psi_{n}^{\dagger}\Big)|0\rangle,~~\psi_{n}^{\dagger}=\sum_{I=1}^{2N_{s}}\Phi_{In}c^{\dagger}_{I}, \end{align} -where $I,J$ denote the site index including the spin degrees of freedom. - -By using the similar argument in the anti-parallel pairing case, -we obtain the following relation. +respectively, where $I$, $J$ denote the site index including the spin degrees of freedom. +By the similar argument as the anti-parallel pairing case, +we can derive the following relation: \begin{align} F_{IJ}=\sum_{n=1}^{\frac{N_{e}}{2}}\Big(\Phi_{In}\Phi_{Jn+1}-\Phi_{Jn}\Phi_{In+1}\Big). \end{align} -Because this relation hold for the case of anti-parallel pairing, -we employ this relation in mVMC ver 1.0. - - +Because this relation hold for the case of anti-parallel pairing, we employ this relation in mVMC ver 1.0 and later. \subsection{Singular value decomposition of $f_{ij}$} +\label{sec:PfaffianSingular} + We define matrices $F$, $\Phi_{\uparrow}$, $\Phi_{\downarrow}$, and $\Sigma$ as \begin{align} &(F)_{ij}=f_{ij},~~~ (\Phi_{\uparrow})_{in}=\Phi_{in\uparrow},~~~ (\Phi_{\downarrow})_{in}=\Phi_{in\downarrow}, \\ -&\Sigma={\rm diag}[1,\cdots,1,0,0,0]~~~\text{({\rm \# of 1} = $N_{e}/2$)}. +&\Sigma={\rm diag}[\underbrace{1,\cdots,1}_{N_e/2},0,0,0]. \end{align} -By using these notations, we can describe the -singular value decomposition of $f_{ij}$ (or equivalently $F$) as +When $f_{ij}$ (i.e., the matrix $F$) is related with a single Slater determinant +of the wave function, we can show that the singular value decomposition of $F$ becomes \begin{align} F=\Phi_{\uparrow}\Sigma\Phi_{\downarrow}^{t}. \end{align} -This result indicates that $f_{ij}$ can be -described by the mean-field solutions -if the number of -nonzero singular values are $N_{e}/2$ and -all the nonzero singular values of $F$ are one. -In other word, the singular values including their numbers -offers the quantitative criterion how the Pfaffian Slater determinant -deviates from the single Slate determinant. +This result indicates that when the number of nonzero singular values is $N_{e}/2$, +and when all the nonzero singular values of $F$ are one in the singular value decomposition of $F$, +the Pfaffian-Slater wave function parametrized by $f_{ij}$ coincides with a single Slater determinant +(i.e. a solution of the mean-field approximation). +In other words, the numbers of the nonzero singular values and their difference from one offer +a quantitative criterion how the Pfaffian-Slater determinant deviates from the single Slate determinant. \section{Power Lanczos method} -In this section, we show how to determine $\alpha$ in the power-Lanczos method. Within the finite- number sampling, it is important to use the decomposition that guarantees the positive definitive relation for variance. We also explain the calculation of physical quantities after the single-step Lanczos method. + +In this section, we show how to determine $\alpha$ in the power-Lanczos method. +We also explain the calculation of physical quantities after the single-step Lanczos method. \subsection{Determination of $\alpha$} First, we briefly explain the sampling procedure of the variational Monte Carlo (VMC) method. Physical properties $\hat{A}$ are calculated as follows: @@ -171,8 +181,9 @@ \subsection{Determination of $\alpha$} For example, we consider the expectation value of the variance, which is defined as $\sigma^2=\langle (\hat{H}-\langle \hat{H}\rangle)^2\rangle$. There are two ways to calculate the variance. \begin{align} -&\sigma^2=\sum_{x} \rho(x) F(x, (\hat{H}-\langle \hat{H}\rangle)^2) = \sum_{x} \rho(x) F(x, \hat{H}^2) - \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 ,\\ -&\sigma^2=\sum_{x} \rho(x) F^{\dag}(x, \hat{H}-\langle \hat{H}\rangle)F(x, \hat{H}-\langle \hat{H}\rangle) = \sum_{x} \rho(x) F^{\dag}(x, \hat{H}) F(x, \hat{H})- \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 +\sigma^2&=\sum_{x} \rho(x) F(x, (\hat{H}-\langle \hat{H}\rangle)^2) = \sum_{x} \rho(x) F(x, \hat{H}^2) - \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 ,\\ +\sigma^2&=\sum_{x} \rho(x) F^{\dag}(x, \hat{H}-\langle \hat{H}\rangle)F(x, \hat{H}-\langle \hat{H}\rangle) \nonumber \\ +&= \sum_{x} \rho(x) F^{\dag}(x, \hat{H}) F(x, \hat{H})- \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 \end{align} From its definition, the latter way gives the positive definitive variance even for the finite sampling while the former way does not guarantee the positive definitiveness of the variance. Here, we consider the expectation values of energy and variance for the (single-step) power Lanczos wave function $|\phi\rangle =(1+\alpha \hat{H}) |\psi \rangle$. The energy is calculated as \begin{align} diff --git a/doc/jp/chap05_jp.tex b/doc/jp/chap05_jp.tex index 09dca2e7..2772fab9 100644 --- a/doc/jp/chap05_jp.tex +++ b/doc/jp/chap05_jp.tex @@ -3,34 +3,46 @@ \chapter{アルゴリズム} \label{Ch:algorithm} \section{変分モンテカルロ法} -適当な完全系に対するマルコフ連鎖を構成して重みつきサンプリングを行います。 -ここでは完全系として$S_z = 0$ の実空間配置$\{| x\rangle\}$ を使用します。 + +変分モンテカルロ法では, 試行波動関数を用意して, その試行波動関数が含むパラメータを +変分原理に従って最適化することで量子多体系の基底状態(または低励起エネルギー状態)の波動関数を +近似的に求めます。 +試行波動関数に対する物理量の期待値を計算する部分で, マルコフ連鎖モンテカルロ法を利用し, +効率よく重み付きサンプリングを行います。 + +本パッケージでは, サンプリングに用いる完全系として電子の実空間配置$\{| x\rangle\}$をとっています: \begin{equation} | x\rangle = \prod_{n=1}^{N/2} c_{r_{n\uparrow}}^{\dag} \prod_{n=1}^{N/2} c_{r_{n\downarrow}}^{\dag} |0 \rangle \end{equation} -ここで、$n$番目の$\sigma$電子の位置を$r_{n\sigma}$としました。 -\subsection{Importance sampling} -マルコフ連鎖の重みを -\begin{equation} -\rho(x)=\frac{|\langle x| \psi \rangle|^2}{\langle \psi | \psi \rangle} \ge 0, \sum{x} \rho(x)=1 -\end{equation} -とすると演算子$A$の期待値は +ここで, $r_{n\sigma}$は$n$番目の電子(スピン$\sigma$)の位置, $c_{r_{n\sigma}}^{\dag}$はその位置での +電子(スピン$\sigma$)の生成演算子を表します。この基底を用いると, 演算子$A$の期待値は \begin{equation} \langle A \rangle =\frac{\langle \psi| A| \psi \rangle}{\langle \psi | \psi \rangle} =\sum_x \frac{\langle \psi| A | x\rangle \langle x| \psi \rangle}{\langle \psi |\psi \rangle} -=\sum_x \rho(x) \frac{\langle \psi| A | x\rangle }{\langle \psi |x \rangle} \end{equation} -となります。$x$に関する和を重み付きサンプリングに置き換えます。また、Local Green's function $G_{ij\sigma\sigma'}(x)$は +となるため, マルコフ連鎖の重みを +\begin{equation} +\rho(x)=\frac{|\langle x| \psi \rangle|^2}{\langle \psi | \psi \rangle} \ge 0, \quad \sum_{x} \rho(x)=1 +\end{equation} +と定義して, +\begin{equation} +\langle A \rangle =\sum_x \rho(x) \frac{\langle \psi| A | x\rangle }{\langle \psi |x \rangle} +\end{equation} +と書き直した後、$x$に関する和をマルコフ連鎖モンテカルロ法により +評価しています。Local Green's function $G_{ij\sigma\sigma'}(x)$は \begin{equation} G_{ij\sigma\sigma'}(x)=\frac{\langle \psi | c_{i\sigma}^{\dag} c_{j\sigma'} | \psi \rangle}{\langle \psi | x \rangle} \end{equation} -で定義されます。なお、サンプリングに使用する乱数生成については、メルセンヌツイスター法を使用しています\cite{Mutsuo2008}。 +と定義されますが, これも演算子$A$を$c_{i\sigma}^{\dag} c_{j\sigma'}$ととることで, +同じ方法により重み付きサンプリングを行うことができます。 +なお, サンプリングに使用する乱数生成については, メルセンヌツイスター法を使用しています\cite{Mutsuo2008}。 \section{Bogoliubov表現}\label{sec_bogoliubov_rep} -スピン系の計算において一体項(\verb|transfer|)、\verb|InterAll|形式での相互作用、 +スピン系の計算において一体項(\verb|transfer|), \verb|InterAll|形式での相互作用, 相関関数のインデックスの指定にはBogoliubov表現が使われています。 -スピンの演算子は次のように生成$\cdot$消滅演算子で書き換えられます。 +一般に、スピンの演算子は次のようにフェルミオンの生成$\cdot$消滅演算子$c_{i \sigma}$, +$c_{i \sigma}^\dagger$によって書き換えることができます: \begin{align} S_{i z} &= \sum_{\sigma = -S}^{S} \sigma c_{i \sigma}^\dagger c_{i \sigma} \\ @@ -42,27 +54,35 @@ \section{Bogoliubov表現}\label{sec_bogoliubov_rep} \sqrt{S(S+1) - \sigma(\sigma+1)} c_{i \sigma}^\dagger c_{i \sigma+1} \end{align} +本パッケージでは、$S=1/2$のスピン系のみ取り扱っており、上記の式で +$S=1/2$と置いたものを用いています。 -\section{{パフィアン行列式とスレータ行列式の関係}} +\section{{パフィアン-スレーター行列式の性質}} \label{sec:PuffAndSlater} -このセクションでは、パフィアン行列式とスレータ行列式の関係および$f_{ij}$の特異値分解の意味について説明します。 + +この節では, パフィアン-スレーター行列式のもつ性質について簡単にまとめます。 +\ref{sec:PfaffianAP}節と\ref{sec:PfaffianP}節でパフィアン-スレーター行列式と単一スレーター行列式の間の関係を導出し、 +\ref{sec:PfaffianSingular}節で$f_{ij}$の特異値分解の意味について説明します。 + \subsection{$f_{ij}$と$\Phi_{in\sigma}$の関係~(スピン反平行の場合)} -パフィアンスレーター行列式(多変数変分モンテカルロ法の一体部分)は +\label{sec:PfaffianAP} + +多変数変分モンテカルロ法で試行波動関数の一体部分として用いられるパフィアン-スレーター行列式は \begin{equation} |\phi_{\rm Pf}\rangle=\Big(\sum_{i,j=1}^{N_{s}}f_{ij}c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \end{equation} -のように定義されます。ここで、$N_{s}$はサイト数、$N_{e}$は全電子数、$f_{ij}$は変分パラメータとしました。 -簡単化のため、以降$f_{ij}$は実数と仮定します。また、シングルスレーター行列式として +のように定義されます。ここで, $N_{s}$はサイト数, $N_{e}$は全電子数, $f_{ij}$は変分パラメータです。 +簡単化のため, 以降$f_{ij}$は実数と仮定します。また, 単一スレーター行列式として \begin{align} |\phi_{\rm SL}\rangle&=\Big(\prod_{n=1}^{N_{e}/2}\psi_{n\uparrow}^{\dagger}\Big) \Big(\prod_{m=1}^{N_{e}/2}\psi_{m\downarrow}^{\dagger}\Big)|0\rangle, \\ \psi_{n\sigma}^{\dagger}&=\sum_{i=1}^{N_{s}}\Phi_{in\sigma}c^{\dagger}_{i\sigma}. \end{align} -を定義します。ただし、$\Phi$は正規直行基底であり、クロネッカーのデルタ$\delta_{nm}$を用い +を定義します。ただし, $\Phi$は正規直交基底であり, クロネッカーのデルタ$\delta_{nm}$を用い \begin{equation} \sum_{i=1}^{N_{s}}\Phi_{in\sigma}\Phi_{im\sigma}=\delta_{nm}, \end{equation} -で表されます。この直交性の関係から、以下の関係式 +で表されます。この直交性の関係から, 以下の関係式 \begin{align} [\psi^{\dagger}_{n\sigma},\psi_{m\sigma}]_{+}&=\delta_{nm},\\ G_{ij\sigma}=\langle c_{i\sigma}^{\dagger}c_{j\sigma}\rangle @@ -71,69 +91,76 @@ \subsection{$f_{ij}$と$\Phi_{in\sigma}$の関係~(スピン反平行の場合)} \end{align} が導かれます。 -次に、$\phi_{\rm SL}$を変形し、$f_{ij}$と$\Phi_{in\sigma}$の関係をあらわにします。 -$\psi^{\dagger}_{n\sigma}$の交換関係を用いると、$\phi_{\rm SL}$は +次に, $|\phi_{\rm SL}\rangle $を変形し, $f_{ij}$と$\Phi_{in\sigma}$の間に成り立つ関係式を示します。 +$\psi^{\dagger}_{n\sigma}$の交換関係を用いると, $|\phi_{\rm SL}\rangle $は \begin{align} |\phi_{\rm SL}\rangle \propto \prod_{n=1}^{N_{e}/2}\Big(\psi_{n\uparrow}^{\dagger}\psi_{\mu(n)\downarrow}^{\dagger}\Big)|0\rangle, \end{align} -と書き換えられます。ここで、$\mu(n)$は$n= 1, 2, \cdots, N_{e}/2$の置換を表します。 -ここで議論を簡単にするため、同一のペア$n=\mu(n)$を採用します。 -このとき、$K_{n}^{\dagger}=\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}$として、 -$K_{n}^{\dagger}K_{m}^{\dagger}=K_{m}^{\dagger}K_{n}^{\dagger}$の関係を用いることで、 +と書き換えられます。ここで, $\mu(n)$は$n= 1, 2, \cdots, N_{e}/2$の置換を表します。 +議論を簡単にするため, 同一のペア$n=\mu(n)$となる場合を考えましょう。 +このとき, $K_{n}^{\dagger}=\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}$として, +$K_{n}^{\dagger}K_{m}^{\dagger}=K_{m}^{\dagger}K_{n}^{\dagger}$の関係を用いることで, \begin{align} |\phi_{\rm SL}\rangle &\propto \prod_{n=1}^{N_{e}/2}\Big(\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}\Big)|0\rangle =\prod_{n=1}^{N_{e}/2} K_{n}^{\dagger}|0\rangle \\ &\propto\Big(\sum_{n=1}^{\frac{N_{e}}{2}}K_{n}^{\dagger}\Big)^{\frac{N_{e}}{2}} |0\rangle =\Big(\sum_{i,j=1}^{N_{s}}\Big[\sum_{n=1}^{\frac{N_{e}}{2}}\Phi_{in\uparrow}\Phi_{jn\downarrow}\Big] -c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)|0\rangle, +c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)^{N_e/2}|0\rangle, \end{align} -の関係が得られます。これより$f_{ij}$はシングルスレーター行列式の係数により +の関係が得られます。これより$f_{ij}$は単一スレーター行列式の係数により \begin{align} f_{ij}=\sum_{n=1}^{\frac{N_{e}}{2}}\Phi_{in\uparrow}\Phi_{jn\downarrow}. \end{align} -として表されることが分かります。なお、この形式はシングルスレーター行列式で与えられる$f_{ij}$の表式の一つであり、 -実際にはペアを組む自由度(どの$\mu(n)$を選ぶか)およびゲージの自由度(すなわち$\Phi$の符号の自由度)に依存します。 +として表されることが分かります。なお, この形式は単一スレーター行列式で与えられる$f_{ij}$の表式の一つであり, +実際にはペアを組む自由度(どの$\mu(n)$を選ぶか)およびゲージの自由度(すなわち$\Phi_{in\sigma}$の符号の自由度)に依存します。 この自由度の多さが$f_{ij}$の冗長性につながっています。 -\subsection{$F_{IJ}$と$\Phi_{In}$の関係~(スピン平行も含めた場合)} -以下の同種スピンのペアリングも考えたパフィアン波動関数と -スレーター波動関数を考えます(ここで$I,J$はスピン自由度も含めたサイトのインデックス). +\subsection{$F_{IJ}$と$\Phi_{In\sigma}$の関係~(スピン平行も含めた場合)} +\label{sec:PfaffianP} + +前節で考察したパフィアン-スレーター波動関数と単一スレーター波動関数の間の関係は、 +同種スピンのペアリングも考えた場合に拡張することができます。 +パフィアン-スレーター波動関数とスレーター波動関数をそれぞれ \begin{align} -|\phi_{\rm Pf}\rangle&=\Big(\sum_{I,J=1}^{2N_{s}}f_{IJ}c_{I}^{\dagger}c_{J}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \\ +|\phi_{\rm Pf}\rangle&=\Big(\sum_{I,J=1}^{2N_{s}}F_{IJ}c_{I}^{\dagger}c_{J}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \\ |\phi_{\rm SL}\rangle&=\Big(\prod_{n=1}^{N_{e}}\psi_{n}^{\dagger}\Big)|0\rangle,~~\psi_{n}^{\dagger}=\sum_{I=1}^{2N_{s}}\Phi_{In}c^{\dagger}_{I}. \end{align} -スピン反平行の場合とほぼ同様の議論を用いることで、 +と定義します。ここで$I,J$はスピン自由度も含めたサイトのインデックスです。 +スピン反平行の場合とほぼ同様の議論を用いることで, \begin{align} F_{IJ}=\sum_{n=1}^{\frac{N_{e}}{2}}\Big(\Phi_{In}\Phi_{Jn+1}-\Phi_{Jn}\Phi_{In+1}\Big). \end{align} -の関係を示すことができます。 -これはスピン反平行の場合にもそのまま適用できるので、 -mVMC ver1.0ではこの表式を使用しています。 +の関係を示すことができます。これはスピン反平行のペアリングにもそのまま適用できるので, +mVMC ver1.0以降ではこの表式を使用しています。 \subsection{$f_{ij}$の特異値分解} -行列$F$, $\Phi_{\uparrow}$, $\Phi_{\downarrow}$,$\Sigma$を +\label{sec:PfaffianSingular} + +行列$F$, $\Phi_{\uparrow}$, $\Phi_{\downarrow}$, $\Sigma$を \begin{align} &(F)_{ij}=f_{ij},~~~ (\Phi_{\uparrow})_{in}=\Phi_{in\uparrow},~~~ (\Phi_{\downarrow})_{in}=\Phi_{in\downarrow}, \\ -&\Sigma={\rm diag}[1,\cdots,1,0,0,0]~~~\text{({\rm \# of 1} = $N_{e}/2$)}. +&\Sigma={\rm diag}[\underbrace{1,\cdots,1}_{N_e/2},0,0,0], \end{align} -として定義します。これらの記法を用いると、$f_{ij}$(すなわち$F$)の特異値分解は +として定義します。前節のように$f_{ij}$(すなわち$F$)が単一スレーター行列と関係づけられて +いるとき、$F$の特異値分解は \begin{align} F=\Phi_{\uparrow}\Sigma\Phi_{\downarrow}^{t}. \end{align} -として記述することができます。 -この結果は、もし非ゼロの特異値が$N_{e}/2$個存在し、 -かつ全ての$F$の非ゼロの特異値が$1$であった場合、 -$f_{ij}$が平均場近似解として記述できることを示しています。 -言い換えると、特異値の非ゼロ成分の数とその値が、 -シングルスレータ行列式からパフィアンスレーター行列式がどのようにしてずれるのか、 +となることを示すことができます。 +この結果は、一般に$F$を特異値分解したとき、非ゼロの特異値が$N_{e}/2$個存在し, +かつ全ての$F$の非ゼロの特異値が$1$であった場合, $f_{ij}$が単一スレーター波動関数を +記述すること(つまり平均場近似解として記述できること)を表しています。 +言い換えると, 特異値の非ゼロ成分の数とその値が, +シングルスレータ行列式からパフィアンスレーター行列式がどのようにしてずれるのか, という点について定量的な基準を与えることを示しています。 \section{Power-Lanczos法} -このセクションでは、Power-Lanczos法での最適化パラメータ$\alpha$の決定方法について述べます。 -サンプリング数が有限のため、エネルギーの分散が常に正の値を取るよう計算を行うことが重要になります。 -また、ここではシングルステップのLanczos法を適用した後の物理量の計算についても説明します。 + +このセクションでは, Power-Lanczos法での最適化パラメータ$\alpha$の決定方法について述べます。 +また, ここではシングルステップのLanczos法を適用した後の物理量の計算についても説明します。 + \subsection{$\alpha$の決定} 最初に, 変分モンテカルロ法のサンプリングに関して簡単に説明します。 物理量$\hat{A}$は以下の手順で計算されます: @@ -141,23 +168,24 @@ \subsection{$\alpha$の決定} &\langle \hat{A}\rangle = \frac{\langle \phi| \hat{A}|\phi \rangle}{\langle \phi| \phi \rangle} = \sum_{x} \rho(x) F(x, {\hat{A}}),\\ & \rho(x)=\frac{|\langle \phi|x\rangle|^2}{\langle \phi | \phi \rangle}, ~~~~F(x, {\hat{A}}) = \frac{\langle x| \hat{A}|\phi \rangle}{\langle x| \phi \rangle}. \end{align} -演算子の積$\hat{A}\hat{B}$を計算する場合には、以下の二通りの方法があります。 +演算子の積$\hat{A}\hat{B}$を計算する場合には, 以下の二通りの方法があります。 \begin{align} &\langle \hat{A} \hat{B}\rangle = \sum_{x} \rho(x) F(x, {\hat{A}\hat{B}}),\\ &\langle \hat{A} \hat{B}\rangle = \sum_{x} \rho(x) F^{\dag}(x, {\hat{A})F(x, \hat{B}}). \end{align} -後述するように、後者の表記の方が数値的には安定します。 -例えば、エネルギーの期待値の分散$\sigma^2=\langle (\hat{H}-\langle \hat{H}\rangle)^2\rangle$を考えてみると、 +後述するように, 後者の表記の方が数値的には安定します。 +例えば, エネルギーの期待値の分散$\sigma^2=\langle (\hat{H}-\langle \hat{H}\rangle)^2\rangle$を考えてみると, 分散は以下の2通りの方法で計算できます。 \begin{align} -&\sigma^2=\sum_{x} \rho(x) F(x, (\hat{H}-\langle \hat{H}\rangle)^2) = \sum_{x} \rho(x) F(x, \hat{H}^2) - \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 ,\\ -&\sigma^2=\sum_{x} \rho(x) F^{\dag}(x, \hat{H}-\langle \hat{H}\rangle)F(x, \hat{H}-\langle \hat{H}\rangle) = \sum_{x} \rho(x) F^{\dag}(x, \hat{H}) F(x, \hat{H})- \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 +\sigma^2 &=\sum_{x} \rho(x) F(x, (\hat{H}-\langle \hat{H}\rangle)^2) = \sum_{x} \rho(x) F(x, \hat{H}^2) - \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 ,\\ +\sigma^2 &=\sum_{x} \rho(x) F^{\dag}(x, \hat{H}-\langle \hat{H}\rangle)F(x, \hat{H}-\langle \hat{H}\rangle) \nonumber \\ +&= \sum_{x} \rho(x) F^{\dag}(x, \hat{H}) F(x, \hat{H})- \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 \end{align} -この定義から、後者の方法では常に正の値となることが保証されているのに対して、前者の方法では分散が正の値になることが必ずしも保証されないことが分かります。次に, シングルステップでのpower-Lanczos波動関数$|\phi\rangle =(1+\alpha \hat{H}) |\psi \rangle$に対するエネルギーの期待値とその分散を考えます。エネルギーは以下の式で計算されます: +この定義から, 後者の方法では常に正の値となることが保証されているのに対して, 前者の方法では分散が正の値になることが必ずしも保証されないことが分かります。次に, シングルステップでのpower-Lanczos波動関数$|\phi\rangle =(1+\alpha \hat{H}) |\psi \rangle$に対するエネルギーの期待値とその分散を考えます。エネルギーは以下の式で計算されます: \begin{align} E_{LS}(\alpha) =\frac{\langle \phi| \hat{H} |\phi\rangle}{\langle \phi|\phi\rangle}=\frac{h_1 + \alpha(h_{2(20)} + h_{2(11)}) + \alpha^2 h_{3(12)}}{1 + 2\alpha h_1 + \alpha^2 h_{2(11)}}, \end{align} -ここで、$h_1$, $h_{2(11)},~h_{2(20)},$ and $h_{3(12)}$を以下のように定義しました: +ここで, $h_1$, $h_{2(11)},~h_{2(20)},$ and $h_{3(12)}$を以下のように定義しました: \begin{align} &h_1 =\sum_{x} \rho(x) F^{\dag}(x, \hat{H}),\\ &h_{2(11)}=\sum_{x} \rho(x) F^{\dag}(x, \hat{H}) F(x, \hat{H}),\\ @@ -171,12 +199,12 @@ \subsection{物理量の計算} \begin{align} A_{LS}(\alpha) =\frac{\langle \phi| \hat{A} |\phi\rangle}{\langle \phi|\phi\rangle}=\frac{A_0 + \alpha(A_{1(10)} + A_{1(01)}) + \alpha^2 A_{2(11)}}{1 + 2\alpha h_1 + \alpha^2 h_{2(11)}}, \end{align} -ここで、$A_0$, $A_{1(10)},~A_{1(01)},$ and $A_{2(11)}$は +ここで, $A_0$, $A_{1(10)},~A_{1(01)},$ and $A_{2(11)}$は \begin{align} &A_0 =\sum_{x} \rho(x) F(x, \hat{A}),\\ &A_{1(10)}=\sum_{x} \rho(x) F^{\dag}(x, \hat{H}) F(x, \hat{A}),\\ &A_{1(01)}=\sum_{x} \rho(x) F(x, \hat{A}\hat{H}),\\ &A_{2(11)}=\sum_{x} \rho(x) F^{\dag}(x, \hat{H})F(x, \hat{A}\hat{H}). \end{align} -として定義される変数を表します。プログラムでは、この表式に基づき一体グリーン関数および二体グリーン関数の計算を行っています。 +として定義される変数を表します。プログラムでは, この表式に基づき一体グリーン関数および二体グリーン関数の計算を行っています。 %---------------------------------------------------------- diff --git a/src/ComplexUHF/output.c b/src/ComplexUHF/output.c index 823293af..92171bff 100644 --- a/src/ComplexUHF/output.c +++ b/src/ComplexUHF/output.c @@ -27,6 +27,7 @@ along with this program. If not, see http://www.gnu.org/licenses/. int MakeOrbitalFile(struct BindStruct *X); void cal_cisajs(struct BindStruct *X); void OutputAntiParallel(struct BindStruct *X,double complex **UHF_Fij,double complex *ParamOrbital,int *CountOrbital); +void OutputAntiParallel_2(struct BindStruct *X,double complex **UHF_Fij,double complex *ParamOrbital,int *CountOrbital); void OutputParallel(struct BindStruct *X,double complex **UHF_Fij,double complex *ParamOrbital,int *CountOrbital); void OutputGeneral(struct BindStruct *X,double complex **UHF_Fij,double complex *ParamOrbital,int *CountOrbital); @@ -122,93 +123,152 @@ int MakeOrbitalFile(struct BindStruct *X){ int isite, jsite; double complex **UHF_Fij; double complex *ParamOrbital; + int *CountOrbital; - //int Orbitalidx; - //int int_i,int_j,int_k,int_l,xMsize; - //double complex **tmp_mat,**vec; - //double *r; - //char fileName[256]; - //int ini,fin,tmp_i; +/*this part only for anti-parallel*/ + int Orbitalidx; + int int_i,int_j,int_k,int_l,xMsize; + double complex **tmp_mat,**vec,**tmp_SLT_U,**tmp_SLT_D,**AP_UHF_fij; + double complex tmp; + double *r; + char fileName[256]; + int ini,fin,tmp_i; -/*this part only for anti-parallel +/*this part only for anti-parallel*/ //[s] for anti-pararell, rediag - xMsize = X->Def.Nsite; - c_malloc2(tmp_mat,xMsize,xMsize); - c_malloc2(vec,xMsize,xMsize); - d_malloc1(r,xMsize); - for(int_l = 0; int_l < 2*xMsize; int_l++){ - for(int_k = 0; int_k < 2*X->Def.Ne; int_k++){ - X->Large.R_SLT[int_l][int_k] = 0.0; - } - } -// for up - for(int_i = 0;int_i < xMsize; int_i++){ - for(int_j = 0;int_j < xMsize; int_j++){ - tmp_mat[int_i][int_j] = X->Large.Ham[int_i][int_j]; - } - } - ZHEEVall(xMsize,tmp_mat,r,vec); - for(int_k = 0; int_k < X->Def.Ne; int_k++){ - for(int_l = 0; int_l < xMsize; int_l++){ - X->Large.R_SLT[int_l][2*int_k] = (vec[int_k][int_l]); - } - } -// for down - for(int_i = 0;int_i < xMsize; int_i++){ - for(int_j = 0;int_j < xMsize; int_j++){ - tmp_mat[int_i][int_j] = X->Large.Ham[int_i+xMsize][int_j+xMsize]; - } - } - ZHEEVall(xMsize,tmp_mat,r,vec); - for(int_k = 0; int_k < X->Def.Ne; int_k++){ - for(int_l = 0; int_l < xMsize; int_l++){ - X->Large.R_SLT[int_l+xMsize][2*int_k+1] = (vec[int_k][int_l]); - } - } -//[e] for anti-pararell -*/ - if(X->Def.NOrbitalIdx>0){ - c_malloc2(UHF_Fij, X->Def.Nsite*2, X->Def.Nsite*2); - for(ispin=0; ispin<2; ispin++){ - for(jspin=0; jspin<2; jspin++){ - for(i=0;i< X->Def.Nsite;i++){ - for(j=0;j< X->Def.Nsite;j++){ - isite = i+ispin*X->Def.Nsite; - jsite = j+jspin*X->Def.Nsite; - UHF_Fij[isite][jsite]=0; - for(n=0;n< 2*X->Def.Ne;n+=2){ - UHF_Fij[isite][jsite] += conj(X->Large.R_SLT[isite][n])*conj(X->Large.R_SLT[jsite][n+1])- conj(X->Large.R_SLT[isite][n+1])*conj(X->Large.R_SLT[jsite][n]); - } + if(X->Def.NOrbitalIdx>0){/*[s]X->Def.NOrbitalIdx>0 */ + if(X->Def.OrbitalOutputMode==1){/*[s]X->Def.OrbitalOutputMode==1 */ + xMsize = X->Def.Nsite; + c_malloc2(tmp_mat,xMsize,xMsize); + c_malloc2(vec,xMsize,xMsize); + c_malloc2(tmp_SLT_U,xMsize,xMsize); + c_malloc2(tmp_SLT_D,xMsize,xMsize); + c_malloc2(AP_UHF_fij,xMsize,xMsize); + d_malloc1(r,xMsize); + for(int_l = 0; int_l < xMsize; int_l++){ + for(int_k = 0; int_k < xMsize; int_k++){ + tmp_SLT_U[int_l][int_k] = 0.0; + tmp_SLT_D[int_l][int_k] = 0.0; + } + } + // for up + for(int_i = 0;int_i < xMsize; int_i++){ + for(int_j = 0;int_j < xMsize; int_j++){ + tmp_mat[int_i][int_j] = X->Large.Ham[int_i][int_j]; + } + } + ZHEEVall(xMsize,tmp_mat,r,vec); + for(int_k = 0; int_k < xMsize; int_k++){ // int_k = n + for(int_l = 0; int_l < xMsize; int_l++){ + tmp_SLT_U[int_l][int_k] = vec[int_k][int_l]; + } + } + // for down + for(int_i = 0;int_i < xMsize; int_i++){ + for(int_j = 0;int_j < xMsize; int_j++){ + tmp_mat[int_i][int_j] = X->Large.Ham[int_i+xMsize][int_j+xMsize]; + } + } + ZHEEVall(xMsize,tmp_mat,r,vec); + for(int_k = 0; int_k < xMsize; int_k++){ // int_k: Ne + for(int_l = 0; int_l < xMsize; int_l++){ + tmp_SLT_D[int_l][int_k] = (vec[int_k][int_l]); + } + } + + for(int_i = 0; int_i < xMsize; int_i++){ // int_k: Ne + for(int_j = 0; int_j < xMsize; int_j++){ + tmp = 0.0; + for(n=0;n< X->Def.Ne;n++){ + tmp += tmp_SLT_U[int_i][n]*tmp_SLT_D[int_j][n]; } + //printf(" %d %d %lf %lf \n",int_i,int_j,creal(tmp),cimag(tmp)); + AP_UHF_fij[int_i][int_j] = tmp; } } - } - //printf(" %d %d %d \n",X->Def.NOrbitalAP,X->Def.NOrbitalP,X->Def.NOrbitalIdx); - if(X->Def.OrbitalOutputMode>0){ + c_malloc1(ParamOrbital, X->Def.NOrbitalIdx); i_malloc1(CountOrbital, X->Def.NOrbitalIdx); -// - OutputAntiParallel(X,UHF_Fij,ParamOrbital,CountOrbital); - if(X->Def.OrbitalOutputMode==2){ - OutputParallel(X,UHF_Fij,ParamOrbital,CountOrbital); - } -// + OutputAntiParallel_2(X,AP_UHF_fij,ParamOrbital,CountOrbital); c_free1(ParamOrbital, X->Def.NOrbitalIdx); i_free1(CountOrbital, X->Def.NOrbitalIdx); - }else if(X->Def.OrbitalOutputMode==0){ - c_malloc1(ParamOrbital, X->Def.NOrbitalIdx); - i_malloc1(CountOrbital, X->Def.NOrbitalIdx); + }/*[s]X->Def.OrbitalOutputMode==1 */ +//[e] for anti-pararell + else{ + c_malloc2(UHF_Fij, X->Def.Nsite*2, X->Def.Nsite*2); + for(ispin=0; ispin<2; ispin++){ + for(jspin=0; jspin<2; jspin++){ + for(i=0;i< X->Def.Nsite;i++){ + for(j=0;j< X->Def.Nsite;j++){ + isite = i+ispin*X->Def.Nsite; + jsite = j+jspin*X->Def.Nsite; + UHF_Fij[isite][jsite]=0; + for(n=0;n< 2*X->Def.Ne;n+=2){ + UHF_Fij[isite][jsite] += conj(X->Large.R_SLT[isite][n])*conj(X->Large.R_SLT[jsite][n+1])- conj(X->Large.R_SLT[isite][n+1])*conj(X->Large.R_SLT[jsite][n]); + } + } + } + } + } + //printf(" %d %d %d \n",X->Def.NOrbitalAP,X->Def.NOrbitalP,X->Def.NOrbitalIdx); + if(X->Def.OrbitalOutputMode==2){ // AP+P + c_malloc1(ParamOrbital, X->Def.NOrbitalIdx); + i_malloc1(CountOrbital, X->Def.NOrbitalIdx); // - OutputGeneral(X,UHF_Fij,ParamOrbital,CountOrbital); + OutputAntiParallel(X,UHF_Fij,ParamOrbital,CountOrbital); + if(X->Def.OrbitalOutputMode==2){ + OutputParallel(X,UHF_Fij,ParamOrbital,CountOrbital); + } // - c_free1(ParamOrbital, X->Def.NOrbitalIdx); - i_free1(CountOrbital, X->Def.NOrbitalIdx); + c_free1(ParamOrbital, X->Def.NOrbitalIdx); + i_free1(CountOrbital, X->Def.NOrbitalIdx); + }else if(X->Def.OrbitalOutputMode==0){ // Only General + c_malloc1(ParamOrbital, X->Def.NOrbitalIdx); + i_malloc1(CountOrbital, X->Def.NOrbitalIdx); +// + OutputGeneral(X,UHF_Fij,ParamOrbital,CountOrbital); +// + c_free1(ParamOrbital, X->Def.NOrbitalIdx); + i_free1(CountOrbital, X->Def.NOrbitalIdx); + } + c_free2(UHF_Fij, X->Def.Nsite * 2, X->Def.Nsite * 2); } - c_free2(UHF_Fij, X->Def.Nsite * 2, X->Def.Nsite * 2); - } + }/*[e]X->Def.NOrbitalIdx>0 */ return 0; } +void OutputAntiParallel_2(struct BindStruct *X,double complex **UHF_Fij,double complex *ParamOrbital,int *CountOrbital){ + int i,j,Orbitalidx,isite,jsite; + char fileName[256]; + + for (i = 0; i < X->Def.NOrbitalIdx; i++) { // all clear + ParamOrbital[i] = 0; + CountOrbital[i] = 0; + } + for(i = 0; i < X->Def.Nsite; i++) { + for(j = 0; j < X->Def.Nsite; j++) { + isite = i + 0 * X->Def.Nsite; + jsite = j + 1 * X->Def.Nsite; + Orbitalidx = X->Def.OrbitalIdx[isite][jsite]; + //printf(" %d %d %d \n", isite,jsite,Orbitalidx); + if(Orbitalidx != -1) { + ParamOrbital[Orbitalidx] += UHF_Fij[i][j]; + CountOrbital[Orbitalidx] += 1; + //printf(" %d %d %d %lf %lf \n", isite,jsite,Orbitalidx,creal(ParamOrbital[Orbitalidx]),cimag(ParamOrbital[Orbitalidx])); + } + } + } + for (i = 0; i < X->Def.NOrbitalAP; i++) { + ParamOrbital[i] /= (double) CountOrbital[i]; + ParamOrbital[i] += genrand_real2() * pow(10.0, -X->Def.eps_int_slater); + } + sprintf(fileName, "%s_APOrbital_opt.dat", X->Def.CParaFileHead); + Child_OutputOptData(fileName, "NOrbitalAP", ParamOrbital, X->Def.NOrbitalAP); + printf("fij for mVMC are outputted to %s.\n", fileName); +} + + + void OutputAntiParallel(struct BindStruct *X,double complex **UHF_Fij,double complex *ParamOrbital,int *CountOrbital){ int i,j,Orbitalidx,isite,jsite; char fileName[256]; diff --git a/src/StdFace/ChainLattice.c b/src/StdFace/ChainLattice.c index f02ee95e..b6b75d67 100644 --- a/src/StdFace/ChainLattice.c +++ b/src/StdFace/ChainLattice.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the chain lattice +*/ #include "StdFace_vals.h" #include #include @@ -24,24 +27,23 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the Hubbard model on a Chain lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Chain(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Hubbard model on a Chain lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Chain( + struct StdIntList *StdI//!<[inout] +) { FILE *fp; - int isite, jsite; + int isite, jsite, ntransMax, nintrMax; int iL; double complex Cphase; - /**/ fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.gp", "w"); /**/ @@ -67,7 +69,9 @@ void StdFace_Chain(struct StdIntList *StdI, char *model) StdI->W = 1; StdFace_InitSite(StdI, fp, 2); StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; StdI->tau[0][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_J("J1", StdI->J1All, StdI->J1); StdFace_NotUsed_J("J2", StdI->J2All, StdI->J2); @@ -122,8 +126,9 @@ void StdFace_Chain(struct StdIntList *StdI, char *model) } }/*if (strcmp(StdI->model, "spin") != 0 )*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->L; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -138,30 +143,28 @@ void StdFace_Chain(struct StdIntList *StdI, char *model) StdI->locspinflag[isite] = StdI->S2; StdI->locspinflag[isite + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->L * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->L * (StdI->NsiteUC/*D*/ + 1/*J*/ + 1/*J'*/) + ntransMax = StdI->L * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->L * (StdI->NsiteUC/*D*/ + 1/*J*/ + 1/*J'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->L * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 2/*t*/ + 2/*t'*/); - StdI->nintr = StdI->L * (StdI->NsiteUC/*U*/ + 4 * (1/*V*/ + 1/*V'*/)); + ntransMax = StdI->L * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 2/*t*/ + 2/*t'*/); + nintrMax = StdI->L * (StdI->NsiteUC/*U*/ + 4 * (1/*V*/ + 1/*V'*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->L * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * 1 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->L * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * 1 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (iL = 0; iL < StdI->L; iL++){ isite = iL; @@ -171,7 +174,7 @@ void StdFace_Chain(struct StdIntList *StdI, char *model) */ if (strcmp(StdI->model, "spin") == 0 ) { StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite, jsite); + StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite, isite); }/*if (strcmp(StdI->model, "spin") == 0 )*/ else { StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite); @@ -214,10 +217,8 @@ void StdFace_Chain(struct StdIntList *StdI, char *model) #if defined(_HPhi) /** -* -* Setup a Hamiltonian for the generalized Heisenberg model on a Chain lattice -* -* @author Mitsuaki Kawamura (The University of Tokyo) +@brief Setup a Hamiltonian for the generalized Heisenberg model on a Chain lattice +@author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_Chain_Boost(struct StdIntList *StdI) { diff --git a/src/StdFace/FCOrtho.c b/src/StdFace/FCOrtho.c index 2b9a6627..fcdd14db 100644 --- a/src/StdFace/FCOrtho.c +++ b/src/StdFace/FCOrtho.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the face centered orthorhombic lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,14 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the Face-Centered Orthorhombic lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_FCOrtho(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Face-Centered Orthorhombic lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_FCOrtho( + struct StdIntList *StdI//!<[inout] +) { - int isite, jsite; + int isite, jsite, ntransMax, nintrMax; int iL, iW, iH, kCell; FILE *fp; double complex Cphase; @@ -39,8 +42,8 @@ void StdFace_FCOrtho(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.xsf", "w"); /**/ @@ -68,7 +71,9 @@ void StdFace_FCOrtho(struct StdIntList *StdI, char *model) /**/ StdFace_InitSite(StdI, fp, 3); StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; ; StdI->tau[0][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_d("K", StdI->K); StdFace_PrintVal_d("h", &StdI->h, 0.0); @@ -137,8 +142,9 @@ void StdFace_FCOrtho(struct StdIntList *StdI, char *model) }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -153,30 +159,28 @@ void StdFace_FCOrtho(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 6/*J*/ + 3/*J'*/ + 0/*J''*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 6/*J*/ + 3/*J'*/ + 0/*J''*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 12/*t*/ + 6/*t'*/ + 0/*t''*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (6/*V*/ + 3/*V'*/ + 0/*V''*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 12/*t*/ + 6/*t'*/ + 0/*t''*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (6/*V*/ + 3/*V'*/ + 0/*V''*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++){ /**/ iW = StdI->Cell[kCell][0]; diff --git a/src/StdFace/HoneycombLattice.c b/src/StdFace/HoneycombLattice.c index b73e8fcb..ded6f609 100644 --- a/src/StdFace/HoneycombLattice.c +++ b/src/StdFace/HoneycombLattice.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the honeycomb lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,12 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the Hubbard model on a Honeycomb lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Honeycomb(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Hubbard model on a Honeycomb lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Honeycomb(struct StdIntList *StdI) { - int isite, jsite, kCell; + int isite, jsite, kCell, ntransMax, nintrMax; int iL, iW; FILE *fp; double complex Cphase; @@ -39,8 +40,8 @@ void StdFace_Honeycomb(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.gp", "w"); /**/ @@ -62,7 +63,9 @@ void StdFace_Honeycomb(struct StdIntList *StdI, char *model) StdFace_InitSite(StdI, fp, 2); StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; StdI->tau[0][2] = 0.0; StdI->tau[1][0] = 1.0 / 3.0; StdI->tau[1][1] = 1.0 / 3.0; StdI->tau[1][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_J("J1'", StdI->J1pAll, StdI->J1p); StdFace_NotUsed_J("J2'", StdI->J2pAll, StdI->J2p); @@ -125,8 +128,9 @@ void StdFace_Honeycomb(struct StdIntList *StdI, char *model) }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -141,30 +145,28 @@ void StdFace_Honeycomb(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 3/*J*/ + 6/*J'*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 3/*J*/ + 6/*J'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 6/*t*/ + 12/*t'*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (3/*V*/ + 6/*V'*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 6/*t*/ + 12/*t'*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (3/*V*/ + 6/*V'*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++) { /**/ iW = StdI->Cell[kCell][0]; diff --git a/src/StdFace/Kagome.c b/src/StdFace/Kagome.c index 338cc512..dc43ad73 100644 --- a/src/StdFace/Kagome.c +++ b/src/StdFace/Kagome.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the kagome lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,14 @@ along with this program. If not, see . #include /** -* -* Setup a Hamiltonian for the Kagome lattice -* -* @author Mitsuaki Kawamura (The University of Tokyo) +@brief Setup a Hamiltonian for the Kagome lattice +@author Mitsuaki Kawamura (The University of Tokyo) */ -void StdFace_Kagome(struct StdIntList *StdI, char *model) +void StdFace_Kagome( + struct StdIntList *StdI//!<[inout] +) { - int isite, jsite, kCell; + int isite, jsite, isiteUC, kCell, ntransMax, nintrMax; int iL, iW; FILE *fp; double complex Cphase; @@ -39,8 +42,8 @@ void StdFace_Kagome(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.gp", "w"); /**/ @@ -63,7 +66,9 @@ void StdFace_Kagome(struct StdIntList *StdI, char *model) StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; StdI->tau[0][2] = 0.0; StdI->tau[1][0] = 0.5; StdI->tau[1][1] = 0.0; StdI->tau[1][2] = 0.0; StdI->tau[2][0] = 0.0; StdI->tau[2][1] = 0.5; StdI->tau[2][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); /**/ StdFace_NotUsed_J("J1'", StdI->J1pAll, StdI->J1p); @@ -124,12 +129,13 @@ void StdFace_Kagome(struct StdIntList *StdI, char *model) StdFace_InputSpin(StdI, StdI->J, StdI->JAll, "J"); }/*if (model != "hubbard")*/ - }/*if (model != "spin")*/ + }/*if (model != "spin")@@*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ - StdI->nsite = 3 * StdI->NCell; + StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; StdI->locspinflag = (int *)malloc(sizeof(int) * StdI->nsite); /**/ @@ -142,64 +148,57 @@ void StdFace_Kagome(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ - if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 6/*J*/ + 6/*J'*/) + if (strcmp(StdI->model, "spin") == 0 ) {//>> + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 6/*J*/ + 6/*J'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 12/*t*/ + 12/*t'*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (6/*V*/ + 6/*V'*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 12/*t*/ + 12/*t'*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (6/*V*/ + 6/*V'*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ - } + }//<< /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++) { /**/ iW = StdI->Cell[kCell][0]; iL = StdI->Cell[kCell][1]; - /* + /*>> Local term */ isite = StdI->NsiteUC * kCell; if (strcmp(StdI->model, "kondo") == 0 ) isite += StdI->nsite / 2; /**/ if (strcmp(StdI->model, "spin") == 0 ) { - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + 1); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + 2); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite, isite); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + 1, isite + 1); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + 2, isite + 2); + for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) { + StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + isiteUC); + StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + isiteUC, isite + isiteUC); + }/*for (jsite = 0; jsite < StdI->NsiteUC; jsite++)*/ }/*if (strcmp(StdI->model, "spin") == 0 )*/ else { - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite); - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + 1); - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + 2); + for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) + StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + isiteUC); /**/ if (strcmp(StdI->model, "kondo") == 0 ) { jsite = StdI->NsiteUC * kCell; - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite, jsite); - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + 1, jsite + 1); - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + 2, jsite + 2); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + 1); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + 2); + for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) { + StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + isiteUC, jsite + isiteUC); + StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + isiteUC); + }/*for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)*/ }/*if (strcmp(StdI->model, "kondo") == 0 )*/ - }/*if (strcmp(StdI->model, "spin") != 0 )*/ - /* + }/*if (strcmp(StdI->model, "spin") != 0 )<<*/ + /*>> Nearest neighbor intra cell 0 -> 1 */ StdFace_SetLabel(StdI, fp, iW, iL, 0, 0, 0, 1, &isite, &jsite, 1, &Cphase); @@ -210,7 +209,7 @@ void StdFace_Kagome(struct StdIntList *StdI, char *model) else { StdFace_Hopping(StdI, Cphase * StdI->t2, isite, jsite); StdFace_Coulomb(StdI, StdI->V2, isite, jsite); - } + }//<< /* Nearest neighbor intra cell 0 -> 2 */ diff --git a/src/StdFace/Ladder.c b/src/StdFace/Ladder.c index 391014d8..ddbb23f5 100644 --- a/src/StdFace/Ladder.c +++ b/src/StdFace/Ladder.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the Ladder lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,23 +27,23 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the generalized Heisenberg model on a square lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Ladder(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the generalized Heisenberg model on a square lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Ladder( + struct StdIntList *StdI//!<[inout] +) { FILE *fp; - int isite, jsite; + int isite, jsite, ntransMax, nintrMax; int iL, isiteUC; double complex Cphase; fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.gp", "w"); /**/ @@ -74,7 +77,9 @@ void StdFace_Ladder(struct StdIntList *StdI, char *model) StdI->tau[isite][0] = (double)isite / (double)StdI->NsiteUC; StdI->tau[isite][1] = 0.0; StdI->tau[isite][2] = 0.0; } - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_J("J", StdI->JAll, StdI->J); StdFace_NotUsed_J("J'", StdI->JpAll, StdI->Jp); @@ -139,8 +144,9 @@ void StdFace_Ladder(struct StdIntList *StdI, char *model) } }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->L * StdI->NsiteUC; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -155,35 +161,33 @@ void StdFace_Ladder(struct StdIntList *StdI, char *model) StdI->locspinflag[isite] = StdI->S2; StdI->locspinflag[isite + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->L * StdI->NsiteUC * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->L * StdI->NsiteUC * (1/*D*/ + 1/*J1*/ + 1/*J1'*/) + ntransMax = StdI->L * StdI->NsiteUC * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->L * StdI->NsiteUC * (1/*D*/ + 1/*J1*/ + 1/*J1'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1) + StdI->L * (StdI->NsiteUC - 1) * (1/*J0*/ + 1/*J2*/ + 1/*J2'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "spin") == 0 )*/ else { - StdI->ntrans = StdI->L*StdI->NsiteUC * (2/*mu+h+Gamma*/ + 2/*t1*/ + 2/*t1'*/) + ntransMax = StdI->L*StdI->NsiteUC * (2/*mu+h+Gamma*/ + 2/*t1*/ + 2/*t1'*/) + StdI->L*(StdI->NsiteUC - 1) * (2/*t0*/ + 2/*t2*/ + 2/*t2'*/); - StdI->nintr = StdI->L*StdI->NsiteUC * 1/*U*/ + nintrMax = StdI->L*StdI->NsiteUC * 1/*U*/ + StdI->L*StdI->NsiteUC * 4 * (1/*V1*/ + 1/*V1'*/) + StdI->L*(StdI->NsiteUC - 1) * 4 * (1/*V0*/ + 1/*V2*/ + 1/*V2'*/); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->L * StdI->NsiteUC * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * 1 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->L * StdI->NsiteUC * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * 1 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (iL = 0; iL < StdI->L; iL++) { for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) { @@ -194,7 +198,7 @@ void StdFace_Ladder(struct StdIntList *StdI, char *model) */ if (strcmp(StdI->model, "spin") == 0 ) { StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite, jsite); + StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite, isite); }/*if (strcmp(StdI->model, "spin") == 0 )*/ else { StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite); diff --git a/src/StdFace/Orthorhombic.c b/src/StdFace/Orthorhombic.c index e6f7cf7f..63aa60ab 100644 --- a/src/StdFace/Orthorhombic.c +++ b/src/StdFace/Orthorhombic.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the orthorhombic lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,14 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the Simple Orthorhombic lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Orthorhombic(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Simple Orthorhombic lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Orthorhombic( + struct StdIntList *StdI//!<[inout] +) { - int isite, jsite; + int isite, jsite, ntransMax, nintrMax; int iL, iW, iH, kCell; FILE *fp; double complex Cphase; @@ -39,8 +42,8 @@ void StdFace_Orthorhombic(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.xsf", "w"); /**/ @@ -68,7 +71,9 @@ void StdFace_Orthorhombic(struct StdIntList *StdI, char *model) /**/ StdFace_InitSite(StdI, fp, 3); StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; ; StdI->tau[0][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_d("K", StdI->K); StdFace_PrintVal_d("h", &StdI->h, 0.0); @@ -137,8 +142,9 @@ void StdFace_Orthorhombic(struct StdIntList *StdI, char *model) }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -153,30 +159,28 @@ void StdFace_Orthorhombic(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 3/*J*/ + 6/*J'*/ + 4/*J''*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 3/*J*/ + 6/*J'*/ + 4/*J''*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 6/*t*/ + 12/*t'*/ + 8/*t''*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (3/*V*/ + 6/*V'*/ + 4/*V''*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 6/*t*/ + 12/*t'*/ + 8/*t''*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (3/*V*/ + 6/*V'*/ + 4/*V''*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++){ /**/ iW = StdI->Cell[kCell][0]; diff --git a/src/StdFace/Pyrochlore.c b/src/StdFace/Pyrochlore.c index 6e5fe2a0..24f7b17e 100644 --- a/src/StdFace/Pyrochlore.c +++ b/src/StdFace/Pyrochlore.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the pyrochlore lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,14 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the Pyrochlore structure - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Pyrochlore(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Pyrochlore structure +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Pyrochlore( + struct StdIntList *StdI//!<[inout] +) { - int isite, jsite; + int isite, jsite, isiteUC, ntransMax, nintrMax; int iL, iW, iH, kCell; FILE *fp; double complex Cphase; @@ -39,8 +42,8 @@ void StdFace_Pyrochlore(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.xsf", "w"); /**/ @@ -71,7 +74,9 @@ void StdFace_Pyrochlore(struct StdIntList *StdI, char *model) StdI->tau[1][0] = 0.5; StdI->tau[1][1] = 0.0; ; StdI->tau[1][2] = 0.0; StdI->tau[2][0] = 0.0; StdI->tau[2][1] = 0.5; ; StdI->tau[2][2] = 0.0; StdI->tau[3][0] = 0.0; StdI->tau[3][1] = 0.0; ; StdI->tau[3][2] = 0.5; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_d("K", StdI->K); StdFace_PrintVal_d("h", &StdI->h, 0.0); @@ -140,8 +145,9 @@ void StdFace_Pyrochlore(struct StdIntList *StdI, char *model) }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -156,30 +162,28 @@ void StdFace_Pyrochlore(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 12/*J*/ + 0/*J'*/ + 0/*J''*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 12/*J*/ + 0/*J'*/ + 0/*J''*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 24/*t*/ + 0/*t'*/ + 0/*t''*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (12/*V*/ + 0/*V'*/ + 0/*V''*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 24/*t*/ + 0/*t'*/ + 0/*t''*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (12/*V*/ + 0/*V'*/ + 0/*V''*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++){ /**/ iW = StdI->Cell[kCell][0]; @@ -192,29 +196,22 @@ void StdFace_Pyrochlore(struct StdIntList *StdI, char *model) if (strcmp(StdI->model, "kondo") == 0) isite += StdI->nsite / 2; /**/ if (strcmp(StdI->model, "spin") == 0) { - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + 1); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + 2); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite, isite); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + 1, isite + 1); - StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + 2, isite + 2); + for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) { + StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + isiteUC); + StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + isiteUC, isite + isiteUC); + }/*for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)*/ }/*if (strcmp(StdI->model, "spin") == 0 )*/ else { - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite); - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + 1); - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + 2); - StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + 3); + for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) { + StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + isiteUC); + }/*for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)*/ /**/ if (strcmp(StdI->model, "kondo") == 0) { jsite = StdI->NsiteUC * kCell; - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite, jsite); - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + 1, jsite + 1); - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + 2, jsite + 2); - StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + 3, jsite + 3); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + 1); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + 2); - StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + 3); + for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) { + StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + 3, jsite + isiteUC); + StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + isiteUC); + }/*for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)*/ }/*if (strcmp(StdI->model, "kondo") == 0 )*/ }/*if (strcmp(StdI->model, "spin") != 0 )*/ /* diff --git a/src/StdFace/SquareLattice.c b/src/StdFace/SquareLattice.c index 2d864c7b..92ea2c6d 100644 --- a/src/StdFace/SquareLattice.c +++ b/src/StdFace/SquareLattice.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the tetragonal lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,12 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the square lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Tetragonal(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the square lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Tetragonal(struct StdIntList *StdI) { - int isite, jsite; + int isite, jsite, ntransMax, nintrMax; int iL, iW, kCell; FILE *fp; double complex Cphase; @@ -39,8 +40,8 @@ void StdFace_Tetragonal(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.gp", "w"); /**/ @@ -61,7 +62,9 @@ void StdFace_Tetragonal(struct StdIntList *StdI, char *model) /**/ StdFace_InitSite(StdI, fp, 2); StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; StdI->tau[0][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_J("J2", StdI->J2All, StdI->J2); StdFace_NotUsed_J("J1'", StdI->J1pAll, StdI->J1p); @@ -120,8 +123,9 @@ void StdFace_Tetragonal(struct StdIntList *StdI, char *model) }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -136,30 +140,28 @@ void StdFace_Tetragonal(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 2/*J*/ + 2/*J'*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 2/*J*/ + 2/*J'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 4/*t*/ + 4/*t'*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (2/*V*/ + 2/*V'*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 4/*t*/ + 4/*t'*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (2/*V*/ + 2/*V'*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++){ /**/ iW = StdI->Cell[kCell][0]; diff --git a/src/StdFace/StdFace_ModelUtil.c b/src/StdFace/StdFace_ModelUtil.c index 91ca210a..dc2a50f0 100644 --- a/src/StdFace/StdFace_ModelUtil.c +++ b/src/StdFace/StdFace_ModelUtil.c @@ -53,7 +53,7 @@ increment StdIntList::ntrans @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_trans( -struct StdIntList *StdI, + struct StdIntList *StdI,//!<[inout] double complex trans0,//!<[in] Hopping integral @f$t, mu@f$, etc. int isite,//!<[in] @f$i@f$ for @f$c_{i \sigma}^\dagger@f$ int ispin,//!<[in] @f$\sigma@f$ for @f$c_{i \sigma}^\dagger@f$ @@ -73,7 +73,7 @@ struct StdIntList *StdI, @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_Hopping( -struct StdIntList *StdI, +struct StdIntList *StdI,//!<[inout] double complex trans0,//!<[in] Hopping integral @f$t@f$ int isite,//!<[in] @f$i@f$ for @f$c_{i \sigma}^\dagger@f$ int jsite//!<[in] @f$j@f$ for @f$c_{j \sigma}@f$ @@ -81,8 +81,8 @@ struct StdIntList *StdI, { int ispin; /**@brief - Both @f$c_{i \sigma}^\daggerc_{j \sigma}@f$ and - @f$c_{j \sigma}^\daggerc_{i \sigma}@f$ for every spin channel + Both @f$c_{i \sigma}^\dagger c_{j \sigma}@f$ and + @f$c_{j \sigma}^\dagger c_{i \sigma}@f$ for every spin channel (@f$\sigma@f$) is specified */ for (ispin = 0; ispin < 2; ispin++) { @@ -96,7 +96,7 @@ itenerant electron @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_HubbardLocal( - struct StdIntList *StdI, + struct StdIntList *StdI,//!<[inout] double mu0,//!<[in] Chemical potential double h0,//!<[in] Longitudinal magnetic feild double Gamma0,//!<[in] Transvers magnetic feild @@ -122,7 +122,7 @@ void StdFace_HubbardLocal( @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_MagField( -struct StdIntList *StdI, + struct StdIntList *StdI,//!<[inout] int S2,//!<[in] Spin moment in @f$i@f$ site double h,//!<[in] Longitudinal magnetic field @f$h@f$ double Gamma,//!<[in] Transvars magnetic field @f$h@f$ @@ -170,7 +170,7 @@ increase the number of that (StdIntList::nintr). @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_intr( -struct StdIntList *StdI, + struct StdIntList *StdI,//!<[inout] double complex intr0,//!<[in] Interaction @f$U, V, J@f$, etc. int site1,//!<[in] @f$i_1@f$ for @f$c_{i_1 \sigma_1}^\dagger@f$ int spin1,//!<[in] @f$sigma1_1@f$ for @f$c_{i_1 \sigma_1}^\dagger@f$ @@ -194,7 +194,7 @@ struct StdIntList *StdI, @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_GeneralJ( -struct StdIntList *StdI, +struct StdIntList *StdI,//!<[inout] double J[3][3],//!<[in] The Spin interaction @f$J_x, J_{xy}@f$, ... int Si2,//!<[in] Spin moment in @f$i@f$ site int Sj2,//!<[in] Spin moment in @f$j@f$ site @@ -303,7 +303,7 @@ struct StdIntList *StdI, StdFace_intr(StdI, conj(intr0), isite, ispin, isite, ispin + 1, jsite, jspin + 1, jsite, jspin); } - /**@brief + /**@brief (3) @f[ I S_i^+ S_j^+ + I^* S_j^- S_i^- = \sum_{\sigma, \sigma' = -S}^{S-1} @@ -323,13 +323,15 @@ struct StdIntList *StdI, StdFace_intr(StdI, conj(intr0), isite, ispin, isite, ispin + 1, jsite, jspin, jsite, jspin + 1); } - /**@brief + /**@brief (4) + @f[ I S_i^+ S_{j z} + I^* S_{j z} S_i^-= \sum_{\sigma=-S}^{S-1} \sum_{\sigma' = -S}^{S} \sqrt{S_i(S_i+1) - \sigma(\sigma+1)} (I c_{i\sigma+1}^\dagger c_{i\sigma} c_{j\sigma'}^\dagger c_{j\sigma'} + I^* c_{j\sigma'}^\dagger c_{j\sigma'} c_{i\sigma}^\dagger c_{i\sigma+1}) \\ I \equiv \frac{J_{xz} - i J_{yz}}{2} + @f] */ if (ispin < Si2) { intr0 = 0.5 * (J[0][2] - I * J[1][2]) * sqrt(Si * (Si + 1.0) - Siz * (Siz + 1.0)) * Sjz; @@ -338,14 +340,16 @@ struct StdIntList *StdI, StdFace_intr(StdI, conj(intr0), jsite, jspin, jsite, jspin, isite, ispin, isite, ispin + 1); }/*if (ispin < Si2)*/ - /**@brief - I S_{i z} S_j^+ + I^* S_j^- S_{i z}= + /**@brief (5) + @f[ + I S_{i z} S_j^+ + I^* S_j^- S_{i z} = \sum_{\sigma=-S}^{S} \sum_{\sigma' = -S}^{S-1} \sqrt{S_j(S_j+1) - \sigma'(\sigma'+1)} (I c_{i\sigma}^\dagger c_{i\sigma} c_{j\sigma'+1}^\dagger c_{j\sigma'} + I^* c_{j\sigma'}^\dagger c_{j\sigma'+1} c_{i\sigma}^\dagger c_{i\sigma}) \\ I \equiv \frac{J_{zx} - i J_{zy}}{2} - */ + @f] + */ if (jspin < Sj2) { intr0 = 0.5 * (J[2][0] - I * J[2][1]) * Siz * sqrt(Sj * (Sj + 1.0) - Sjz * (Sjz + 1.0)); StdFace_intr(StdI, intr0, @@ -363,7 +367,7 @@ and increase the number of them (StdIntList::NCinter). @author Mitsuaki Kawamura (The University of Tokyo) */ void StdFace_Coulomb( -struct StdIntList *StdI, +struct StdIntList *StdI,//!<[inout] double V,//!<[in] Coulomb integral U, V, etc. int isite,//!<[in] i of n_i int jsite//!<[in] j of n_j @@ -561,10 +565,11 @@ void StdFace_RequiredVal_i( original supercell. @author Mitsuaki Kawamura (The University of Tokyo) */ -static void StdFace_FoldSite(struct StdIntList *StdI, +static void StdFace_FoldSite( + struct StdIntList *StdI,//!<[inout] int iCellV[3],//!<[in] The fractional coordinate of a site int nBox[3], //!<[out] the index of supercell - int iCellV_fold[3]/**<@brief [out] The fractional coordinate of a site + int iCellV_fold[3]/**<[out] The fractional coordinate of a site which is moved into the original cell*/ ) { @@ -597,7 +602,8 @@ static void StdFace_FoldSite(struct StdIntList *StdI, @brief Initialize the super-cell where simulation is performed. @author Mitsuaki Kawamura (The University of Tokyo) */ -void StdFace_InitSite(struct StdIntList *StdI, +void StdFace_InitSite( + struct StdIntList *StdI,//!<[inout] FILE *fp,//!<[in] File pointer to lattice.gp int dim//!<[in] dimension of system, if = 2, print lattice.gp ) @@ -787,7 +793,8 @@ void StdFace_InitSite(struct StdIntList *StdI, /** @brief Find the index of transfer and interaction */ -void StdFace_FindSite(struct StdIntList *StdI, +void StdFace_FindSite( + struct StdIntList *StdI,//!<[inout] int iW,//!<[in] position of initial site int iL,//!<[in] position of initial site int iH,//!<[in] position of initial site @@ -835,8 +842,9 @@ void StdFace_FindSite(struct StdIntList *StdI, /** @brief Set Label in the gnuplot display (Only used in 2D system) */ -void StdFace_SetLabel(struct StdIntList *StdI, - FILE *fp, +void StdFace_SetLabel( + struct StdIntList *StdI,//!<[inout] + FILE *fp,//!<[in] File pointer to lattice.gp int iW,//!<[in] position of initial site int iL,//!<[in] position of initial site int diW,//!<[in] Translation from the initial site @@ -929,7 +937,8 @@ void StdFace_PrintXSF(struct StdIntList *StdI) { /** @brief Input nearest-neighbor spin-spin interaction */ -void StdFace_InputSpinNN(struct StdIntList *StdI, +void StdFace_InputSpinNN( + struct StdIntList *StdI,//!<[inout] double J0[3][3],//!<[in] The anisotropic spin interaction double J0All,//!<[in] The isotropic interaction char *J0name//!<[in] The name of this spin interaction (e.g. J1) @@ -1015,7 +1024,8 @@ void StdFace_InputSpinNN(struct StdIntList *StdI, /** @brief Input spin-spin interaction other than nearest-neighbor */ -void StdFace_InputSpin(struct StdIntList *StdI, +void StdFace_InputSpin( + struct StdIntList *StdI,//!<[inout] double Jp[3][3],//!<[in] Fully anisotropic spin interaction double JpAll,//!<[in] The isotropic interaction char *Jpname//!transindx = (int **)malloc(sizeof(int*) * StdI->ntrans); - StdI->trans = (double complex *)malloc(sizeof(double complex) * StdI->ntrans); - for (ii = 0; ii < StdI->ntrans; ii++) { + StdI->transindx = (int **)malloc(sizeof(int*) * ntransMax); + StdI->trans = (double complex *)malloc(sizeof(double complex) * ntransMax); + for (ii = 0; ii < ntransMax; ii++) { StdI->transindx[ii] = (int *)malloc(sizeof(int) * 4); } + StdI->ntrans = 0; /**@brief (2) InterAll StdIntList::intr, StdIntList::intrindx */ - StdI->intrindx = (int **)malloc(sizeof(int*) * StdI->nintr); - StdI->intr = (double complex *)malloc(sizeof(double complex) * StdI->nintr); - for (ii = 0; ii < StdI->nintr; ii++) { + StdI->intrindx = (int **)malloc(sizeof(int*) * nintrMax); + StdI->intr = (double complex *)malloc(sizeof(double complex) * nintrMax); + for (ii = 0; ii < nintrMax; ii++) { StdI->intrindx[ii] = (int *)malloc(sizeof(int) * 8); } + StdI->nintr = 0; /**@brief (3) Coulomb intra StdIntList::Cintra, StdIntList::CintraIndx */ - StdI->CintraIndx = (int **)malloc(sizeof(int*) * StdI->nintr); - StdI->Cintra = (double *)malloc(sizeof(double) * StdI->nintr); - for (ii = 0; ii < StdI->nintr; ii++) { + StdI->CintraIndx = (int **)malloc(sizeof(int*) * nintrMax); + StdI->Cintra = (double *)malloc(sizeof(double) * nintrMax); + for (ii = 0; ii < nintrMax; ii++) { StdI->CintraIndx[ii] = (int *)malloc(sizeof(int) * 1); } + StdI->NCintra = 0; /**@brief (4) Coulomb inter StdIntList::Cinter, StdIntList::CinterIndx */ - StdI->CinterIndx = (int **)malloc(sizeof(int*) * StdI->nintr); - StdI->Cinter = (double *)malloc(sizeof(double) * StdI->nintr); - for (ii = 0; ii < StdI->nintr; ii++) { + StdI->CinterIndx = (int **)malloc(sizeof(int*) * nintrMax); + StdI->Cinter = (double *)malloc(sizeof(double) * nintrMax); + for (ii = 0; ii < nintrMax; ii++) { StdI->CinterIndx[ii] = (int *)malloc(sizeof(int) * 2); } + StdI->NCinter = 0; /**@brief (5) Hund StdIntList::Hund, StdIntList::HundIndx */ - StdI->HundIndx = (int **)malloc(sizeof(int*) * StdI->nintr); - StdI->Hund = (double *)malloc(sizeof(double) * StdI->nintr); - for (ii = 0; ii < StdI->nintr; ii++) { + StdI->HundIndx = (int **)malloc(sizeof(int*) * nintrMax); + StdI->Hund = (double *)malloc(sizeof(double) * nintrMax); + for (ii = 0; ii < nintrMax; ii++) { StdI->HundIndx[ii] = (int *)malloc(sizeof(int) * 2); } + StdI->NHund = 0; /**@brief (6) Excahnge StdIntList::Ex, StdIntList::ExIndx */ - StdI->ExIndx = (int **)malloc(sizeof(int*) * StdI->nintr); - StdI->Ex = (double *)malloc(sizeof(double) * StdI->nintr); - for (ii = 0; ii < StdI->nintr; ii++) { + StdI->ExIndx = (int **)malloc(sizeof(int*) * nintrMax); + StdI->Ex = (double *)malloc(sizeof(double) * nintrMax); + for (ii = 0; ii < nintrMax; ii++) { StdI->ExIndx[ii] = (int *)malloc(sizeof(int) * 2); } + StdI->NEx = 0; /**@brief (7) PairLift StdIntList::PairLift, StdIntList::PLIndx */ - StdI->PLIndx = (int **)malloc(sizeof(int*) * StdI->nintr); - StdI->PairLift = (double *)malloc(sizeof(double) * StdI->nintr); - for (ii = 0; ii < StdI->nintr; ii++) { + StdI->PLIndx = (int **)malloc(sizeof(int*) * nintrMax); + StdI->PairLift = (double *)malloc(sizeof(double) * nintrMax); + for (ii = 0; ii < nintrMax; ii++) { StdI->PLIndx[ii] = (int *)malloc(sizeof(int) * 2); } - - StdI->NCintra = 0; - StdI->NCinter = 0; - StdI->NHund = 0; - StdI->NEx = 0; StdI->NPairLift = 0; }/*void StdFace_MallocInteractions*/ #if defined(_mVMC) @@ -1217,7 +1234,8 @@ void StdFace_MallocInteractions(struct StdIntList *StdI) { @brief Define whether the specified site is in the unit cell or not. @author Mitsuaki Kawamura (The University of Tokyo) */ -static void StdFace_FoldSiteSub(struct StdIntList *StdI, +static void StdFace_FoldSiteSub( + struct StdIntList *StdI,//!<[inout] int iCellV[3],//!<[in] int nBox[3], //!<[out] int iCellV_fold[3]//!<[out] @@ -1576,60 +1594,60 @@ void PrintJastrow(struct StdIntList *StdI) { Jastrow[isite] = (int *)malloc(sizeof(int) * StdI->nsite); if (abs(StdI->NMPTrans) == 1 || StdI->NMPTrans == StdI->NaN_i) { - /**@brief - (1) Copy Orbital index - */ - for (isite = 0; isite < StdI->nsite; isite++) { - for (jsite = 0; jsite < StdI->nsite; jsite++) { - Jastrow[isite][jsite] = StdI->Orb[isite][jsite]; - }/*for (jsite = 0; jsite < isite; jsite++)*/ - }/*for (isite = 0; isite < StdI->nsite; isite++)*/ - /**@brief - (2) Symmetrize - */ - for (iorb = 0; iorb < StdI->NOrb; iorb++) { - for (isite = 0; isite < StdI->nsite; isite++) { - for (jsite = 0; jsite < StdI->nsite; jsite++) { - if (Jastrow[isite][jsite] == iorb) { - Jastrow[jsite][isite] = Jastrow[isite][jsite]; - } - }/*for (jsite = 0; jsite < isite; jsite++)*/ - }/*for (isite = 0; isite < StdI->nsite; isite++)*/ - }/*for (iorb = 0; iorb < StdI->NOrb; iorb++)*/ - /**/ - if (strcmp(StdI->model, "hubbard") == 0) NJastrow = 0; - else NJastrow = -1; - for (isite = 0; isite < StdI->nsite; isite++) { - /* - For Local spin - */ - if (StdI->locspinflag[isite] != 0) { - for (jsite = 0; jsite < StdI->nsite; jsite++) { - Jastrow[isite][jsite] = -1; - Jastrow[jsite][isite] = -1; - } - continue; - } - /**/ - for (jsite = 0; jsite < isite; jsite++) { - if (Jastrow[isite][jsite] >= 0) { - iJastrow = Jastrow[isite][jsite]; - NJastrow -= 1; - for (isite1 = 0; isite1 < StdI->nsite; isite1++) { - for (jsite1 = 0; jsite1 < StdI->nsite; jsite1++) { - if (Jastrow[isite1][jsite1] == iJastrow) - Jastrow[isite1][jsite1] = NJastrow; - }/*for (jsite1 = 0; jsite1 < StdI->nsite; jsite1++)*/ - }/*for (isite1 = 0; isite1 < StdI->nsite; isite1++)*/ - }/*if (Jastrow[isite][jsite] >= 0)*/ - }/*for (jsite = 0; jsite < isite; jsite++)*/ - }/*for (isite = 0; isite < StdI->nsite; isite++)*/ - /**/ - NJastrow = -NJastrow; - for (isite = 0; isite < StdI->nsite; isite++) { - for (jsite = 0; jsite < StdI->nsite; jsite++) { - Jastrow[isite][jsite] = -1 - Jastrow[isite][jsite]; - }/*for (jsite = 0; jsite < isite; jsite++)*/ + /**@brief + (1) Copy Orbital index + */ + for (isite = 0; isite < StdI->nsite; isite++) { + for (jsite = 0; jsite < StdI->nsite; jsite++) { + Jastrow[isite][jsite] = StdI->Orb[isite][jsite]; + }/*for (jsite = 0; jsite < isite; jsite++)*/ + }/*for (isite = 0; isite < StdI->nsite; isite++)*/ + /**@brief + (2) Symmetrize + */ + for (iorb = 0; iorb < StdI->NOrb; iorb++) { + for (isite = 0; isite < StdI->nsite; isite++) { + for (jsite = 0; jsite < StdI->nsite; jsite++) { + if (Jastrow[isite][jsite] == iorb) { + Jastrow[jsite][isite] = Jastrow[isite][jsite]; + } + }/*for (jsite = 0; jsite < isite; jsite++)*/ + }/*for (isite = 0; isite < StdI->nsite; isite++)*/ + }/*for (iorb = 0; iorb < StdI->NOrb; iorb++)*/ + /**/ + if (strcmp(StdI->model, "hubbard") == 0) NJastrow = 0; + else NJastrow = -1; + for (isite = 0; isite < StdI->nsite; isite++) { + /* + For Local spin + */ + if (StdI->locspinflag[isite] != 0) { + for (jsite = 0; jsite < StdI->nsite; jsite++) { + Jastrow[isite][jsite] = -1; + Jastrow[jsite][isite] = -1; + } + continue; + } + /**/ + for (jsite = 0; jsite < isite; jsite++) { + if (Jastrow[isite][jsite] >= 0) { + iJastrow = Jastrow[isite][jsite]; + NJastrow -= 1; + for (isite1 = 0; isite1 < StdI->nsite; isite1++) { + for (jsite1 = 0; jsite1 < StdI->nsite; jsite1++) { + if (Jastrow[isite1][jsite1] == iJastrow) + Jastrow[isite1][jsite1] = NJastrow; + }/*for (jsite1 = 0; jsite1 < StdI->nsite; jsite1++)*/ + }/*for (isite1 = 0; isite1 < StdI->nsite; isite1++)*/ + }/*if (Jastrow[isite][jsite] >= 0)*/ + }/*for (jsite = 0; jsite < isite; jsite++)*/ + }/*for (isite = 0; isite < StdI->nsite; isite++)*/ + /**/ + NJastrow = -NJastrow; + for (isite = 0; isite < StdI->nsite; isite++) { + for (jsite = 0; jsite < StdI->nsite; jsite++) { + Jastrow[isite][jsite] = -1 - Jastrow[isite][jsite]; + }/*for (jsite = 0; jsite < isite; jsite++)*/ }/*for (isite = 0; isite < StdI->nsite; isite++)*/ }/*if (abs(StdI->NMPTrans) == 1)*/ else { diff --git a/src/StdFace/StdFace_ModelUtil.h b/src/StdFace/StdFace_ModelUtil.h index 0b2a63ad..e063b6bf 100644 --- a/src/StdFace/StdFace_ModelUtil.h +++ b/src/StdFace/StdFace_ModelUtil.h @@ -56,23 +56,23 @@ void StdFace_SetLabel(struct StdIntList *StdI, FILE *fp, int iW, int iL, int diW, int diL, int isiteUC, int jsiteUC, int *isite, int *jsite, int connect, double complex *Cphase); void StdFace_PrintGeometry(struct StdIntList *StdI); -void StdFace_MallocInteractions(struct StdIntList *StdI); +void StdFace_MallocInteractions(struct StdIntList *StdI, int ntransMax, int nintrMax); void StdFace_FindSite(struct StdIntList *StdI, int iW, int iL, int iH, int diW, int diL, int diH, int isiteUC, int jsiteUC, int *isite, int *jsite, double complex *Cphase); void StdFace_PrintXSF(struct StdIntList *StdI); -void StdFace_Tetragonal(struct StdIntList *StdI, char *model); -void StdFace_Chain(struct StdIntList *StdI, char *model); -void StdFace_Ladder(struct StdIntList *StdI, char *model); -void StdFace_Triangular(struct StdIntList *StdI, char *model); -void StdFace_Honeycomb(struct StdIntList *StdI, char *model); -void StdFace_Kagome(struct StdIntList *StdI, char *model); -void StdFace_Orthorhombic(struct StdIntList *StdI, char *model); -void StdFace_FCOrtho(struct StdIntList *StdI, char *model); -void StdFace_Pyrochlore(struct StdIntList *StdI, char *model); -void StdFace_Wannier90(struct StdIntList *StdI, char *model); +void StdFace_Tetragonal(struct StdIntList *StdI); +void StdFace_Chain(struct StdIntList *StdI); +void StdFace_Ladder(struct StdIntList *StdI); +void StdFace_Triangular(struct StdIntList *StdI); +void StdFace_Honeycomb(struct StdIntList *StdI); +void StdFace_Kagome(struct StdIntList *StdI); +void StdFace_Orthorhombic(struct StdIntList *StdI); +void StdFace_FCOrtho(struct StdIntList *StdI); +void StdFace_Pyrochlore(struct StdIntList *StdI); +void StdFace_Wannier90(struct StdIntList *StdI); #if defined(_HPhi) void StdFace_Chain_Boost(struct StdIntList *StdI); diff --git a/src/StdFace/StdFace_main.c b/src/StdFace/StdFace_main.c index 0a5e85c7..4fee32bd 100644 --- a/src/StdFace/StdFace_main.c +++ b/src/StdFace/StdFace_main.c @@ -19,6 +19,18 @@ along with this program. If not, see . @brief Read Input file and write files for Expert mode. Initialize variables. Check parameters. + +The following lattices are supported: +- 1D Chain : StdFace_Chain() +- 1D Ladder : StdFace_Ladder() +- 2D Tetragonal : StdFace_Tetragonal() +- 2D Triangular : StdFace_Triangular() +- 2D Honeycomb : StdFace_Honeycomb() +- 2D Kagome : StdFace_Kagome() +- 3D Simple Orthorhombic : StdFace_Orthorhombic() +- 3D Face Centered Orthorhombic : StdFace_FCOrtho() +- 3D Pyrochlore : StdFace_Pyrochlore() + */ #include #include @@ -485,33 +497,33 @@ static void PrintGutzwiller(struct StdIntList *StdI) Gutz = (int *)malloc(sizeof(int) * StdI->nsite); if (abs(StdI->NMPTrans) == 1 || StdI->NMPTrans == StdI->NaN_i) { - if (strcmp(StdI->model, "hubbard") == 0) NGutzwiller = 0; - else NGutzwiller = -1; - - for (isite = 0; isite < StdI->nsite; isite++) Gutz[isite] = StdI->Orb[isite][isite]; - - for (isite = 0; isite < StdI->nsite; isite++) { - /* - For Local spin - */ - if (StdI->locspinflag[isite] != 0) { - Gutz[isite] = -1; - continue; - } - /**/ - if (Gutz[isite] >= 0) { - iGutz = Gutz[isite]; - NGutzwiller -= 1; - for (jsite = 0; jsite < StdI->nsite; jsite++) { - if (Gutz[jsite] == iGutz) - Gutz[jsite] = NGutzwiller; - }/*for (jsite = 0; jsite < StdI->nsite; jsite++)*/ - }/*if (Gutz[isite] >= 0)*/ - }/*for (isite = 0; isite < StdI->nsite; isite++)*/ - /**/ - NGutzwiller = -NGutzwiller; - for (isite = 0; isite < StdI->nsite; isite++) { - Gutz[isite] = -1 - Gutz[isite]; + if (strcmp(StdI->model, "hubbard") == 0) NGutzwiller = 0; + else NGutzwiller = -1; + + for (isite = 0; isite < StdI->nsite; isite++) Gutz[isite] = StdI->Orb[isite][isite]; + + for (isite = 0; isite < StdI->nsite; isite++) { + /* + For Local spin + */ + if (StdI->locspinflag[isite] != 0) { + Gutz[isite] = -1; + continue; + } + /**/ + if (Gutz[isite] >= 0) { + iGutz = Gutz[isite]; + NGutzwiller -= 1; + for (jsite = 0; jsite < StdI->nsite; jsite++) { + if (Gutz[jsite] == iGutz) + Gutz[jsite] = NGutzwiller; + }/*for (jsite = 0; jsite < StdI->nsite; jsite++)*/ + }/*if (Gutz[isite] >= 0)*/ + }/*for (isite = 0; isite < StdI->nsite; isite++)*/ + /**/ + NGutzwiller = -NGutzwiller; + for (isite = 0; isite < StdI->nsite; isite++) { + Gutz[isite] = -1 - Gutz[isite]; }/*for (isite = 0; isite < StdI->nsite; isite++)*/ }/*if (abs(StdI->NMPTrans) == 1)*/ else { @@ -1788,8 +1800,7 @@ void StdFace_main( char *fname//!<[in] Input file name for the standard mode ) { - struct StdIntList StdI; - + struct StdIntList StdI[1]; FILE *fp; int ktrans, kintr; char ctmpline[256]; @@ -1804,7 +1815,7 @@ void StdFace_main( fprintf(stdout, "\n Open Standard-Mode Inputfile %s \n\n", fname); } - StdFace_ResetVals(&StdI); + StdFace_ResetVals(StdI); while (fgets(ctmpline, 256, fp) != NULL) { @@ -1826,217 +1837,217 @@ void StdFace_main( Text2Lower(keyword); fprintf(stdout, " KEYWORD : %-20s | VALUE : %s \n", keyword, value); - if (strcmp(keyword, "a") == 0) StoreWithCheckDup_d(keyword, value, &StdI.a); - else if (strcmp(keyword, "a0h") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[0][2]); - else if (strcmp(keyword, "a0l") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[0][1]); - else if (strcmp(keyword, "a0w") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[0][0]); - else if (strcmp(keyword, "a1h") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[1][2]); - else if (strcmp(keyword, "a1l") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[1][1]); - else if (strcmp(keyword, "a1w") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[1][0]); - else if (strcmp(keyword, "a2h") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[2][2]); - else if (strcmp(keyword, "a2l") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[2][1]); - else if (strcmp(keyword, "a2w") == 0) StoreWithCheckDup_i(keyword, value, &StdI.box[2][0]); - else if (strcmp(keyword, "d") == 0) StoreWithCheckDup_d(keyword, value, &StdI.D[2][2]); - else if (strcmp(keyword, "gamma") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Gamma); - else if (strcmp(keyword, "h") == 0) StoreWithCheckDup_d(keyword, value, &StdI.h); - else if (strcmp(keyword, "height") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Height); - else if (strcmp(keyword, "hlength") == 0) StoreWithCheckDup_d(keyword, value, &StdI.length[2]); - else if (strcmp(keyword, "hx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[2][0]); - else if (strcmp(keyword, "hy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[2][1]); - else if (strcmp(keyword, "hz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[2][2]); - else if (strcmp(keyword, "j") == 0) StoreWithCheckDup_d(keyword, value, &StdI.JAll); - else if (strcmp(keyword, "jx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[0][0]); - else if (strcmp(keyword, "jxy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[0][1]); - else if (strcmp(keyword, "jxz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[0][2]); - else if (strcmp(keyword, "jy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[1][1]); - else if (strcmp(keyword, "jyx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[1][0]); - else if (strcmp(keyword, "jyz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[1][2]); - else if (strcmp(keyword, "jz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[2][2]); - else if (strcmp(keyword, "jzx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[2][0]); - else if (strcmp(keyword, "jzy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J[2][1]); - else if (strcmp(keyword, "j0") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0All); - else if (strcmp(keyword, "j0x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[0][0]); - else if (strcmp(keyword, "j0xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[0][1]); - else if (strcmp(keyword, "j0xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[0][2]); - else if (strcmp(keyword, "j0y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[1][1]); - else if (strcmp(keyword, "j0yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[1][0]); - else if (strcmp(keyword, "j0yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[1][2]); - else if (strcmp(keyword, "j0z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[2][2]); - else if (strcmp(keyword, "j0zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[2][0]); - else if (strcmp(keyword, "j0zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0[2][1]); - else if (strcmp(keyword, "j0'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0pAll); - else if (strcmp(keyword, "j0'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[0][0]); - else if (strcmp(keyword, "j0'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[0][1]); - else if (strcmp(keyword, "j0'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[0][2]); - else if (strcmp(keyword, "j0'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[1][1]); - else if (strcmp(keyword, "j0'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[1][0]); - else if (strcmp(keyword, "j0'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[1][2]); - else if (strcmp(keyword, "j0'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[2][2]); - else if (strcmp(keyword, "j0'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[2][0]); - else if (strcmp(keyword, "j0'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J0p[2][1]); - else if (strcmp(keyword, "j1") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1All); - else if (strcmp(keyword, "j1x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[0][0]); - else if (strcmp(keyword, "j1xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[0][1]); - else if (strcmp(keyword, "j1xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[0][2]); - else if (strcmp(keyword, "j1y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[1][1]); - else if (strcmp(keyword, "j1yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[1][0]); - else if (strcmp(keyword, "j1yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[1][2]); - else if (strcmp(keyword, "j1z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[2][2]); - else if (strcmp(keyword, "j1zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[2][0]); - else if (strcmp(keyword, "j1zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1[2][1]); - else if (strcmp(keyword, "j1'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1pAll); - else if (strcmp(keyword, "j1'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[0][0]); - else if (strcmp(keyword, "j1'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[0][1]); - else if (strcmp(keyword, "j1'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[0][2]); - else if (strcmp(keyword, "j1'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[1][1]); - else if (strcmp(keyword, "j1'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[1][0]); - else if (strcmp(keyword, "j1'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[1][2]); - else if (strcmp(keyword, "j1'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[2][2]); - else if (strcmp(keyword, "j1'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[2][0]); - else if (strcmp(keyword, "j1'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J1p[2][1]); - else if (strcmp(keyword, "j2") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2All); - else if (strcmp(keyword, "j2x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[0][0]); - else if (strcmp(keyword, "j2xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[0][1]); - else if (strcmp(keyword, "j2xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[0][2]); - else if (strcmp(keyword, "j2y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[1][1]); - else if (strcmp(keyword, "j2yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[1][0]); - else if (strcmp(keyword, "j2yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[1][2]); - else if (strcmp(keyword, "j2z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[2][2]); - else if (strcmp(keyword, "j2zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[2][0]); - else if (strcmp(keyword, "j2zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2[2][1]); - else if (strcmp(keyword, "j2'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2pAll); - else if (strcmp(keyword, "j2'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[0][0]); - else if (strcmp(keyword, "j2'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[0][1]); - else if (strcmp(keyword, "j2'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[0][2]); - else if (strcmp(keyword, "j2'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[1][1]); - else if (strcmp(keyword, "j2'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[1][0]); - else if (strcmp(keyword, "j2'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[1][2]); - else if (strcmp(keyword, "j2'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[2][2]); - else if (strcmp(keyword, "j2'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[2][0]); - else if (strcmp(keyword, "j2'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.J2p[2][1]); - else if (strcmp(keyword, "j'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.JpAll); - else if (strcmp(keyword, "j'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[0][0]); - else if (strcmp(keyword, "j'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[0][1]); - else if (strcmp(keyword, "j'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[0][2]); - else if (strcmp(keyword, "j'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[1][1]); - else if (strcmp(keyword, "j'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[1][0]); - else if (strcmp(keyword, "j'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[1][2]); - else if (strcmp(keyword, "j'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[2][2]); - else if (strcmp(keyword, "j'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[2][0]); - else if (strcmp(keyword, "j'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jp[2][1]); - else if (strcmp(keyword, "j''") == 0) StoreWithCheckDup_d(keyword, value, &StdI.JppAll); - else if (strcmp(keyword, "j''x") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[0][0]); - else if (strcmp(keyword, "j''xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[0][1]); - else if (strcmp(keyword, "j''xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[0][2]); - else if (strcmp(keyword, "j''y") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[1][1]); - else if (strcmp(keyword, "j''yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[1][0]); - else if (strcmp(keyword, "j''yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[1][2]); - else if (strcmp(keyword, "j''z") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[2][2]); - else if (strcmp(keyword, "j''zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[2][0]); - else if (strcmp(keyword, "j''zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Jpp[2][1]); - else if (strcmp(keyword, "k") == 0) StoreWithCheckDup_d(keyword, value, &StdI.K); - else if (strcmp(keyword, "l") == 0) StoreWithCheckDup_i(keyword, value, &StdI.L); - else if (strcmp(keyword, "lattice") == 0) StoreWithCheckDup_sl(keyword, value, StdI.lattice); - else if (strcmp(keyword, "llength") == 0) StoreWithCheckDup_d(keyword, value, &StdI.length[1]); - else if (strcmp(keyword, "lx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[1][0]); - else if (strcmp(keyword, "ly") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[1][1]); - else if (strcmp(keyword, "lz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[1][2]); - else if (strcmp(keyword, "model") == 0) StoreWithCheckDup_sl(keyword, value, StdI.model); - else if (strcmp(keyword, "mu") == 0) StoreWithCheckDup_d(keyword, value, &StdI.mu); - else if (strcmp(keyword, "nelec") == 0) StoreWithCheckDup_i(keyword, value, &StdI.nelec); - else if (strcmp(keyword, "outputmode") == 0) StoreWithCheckDup_sl(keyword, value, StdI.outputmode); - else if (strcmp(keyword, "phase0") == 0) StoreWithCheckDup_d(keyword, value, &StdI.phase[0]); - else if (strcmp(keyword, "phase1") == 0) StoreWithCheckDup_d(keyword, value, &StdI.phase[1]); - else if (strcmp(keyword, "phase2") == 0) StoreWithCheckDup_d(keyword, value, &StdI.phase[2]); - else if (strcmp(keyword, "t") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t); - else if (strcmp(keyword, "t0") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t0); - else if (strcmp(keyword, "t0'") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t0p); - else if (strcmp(keyword, "t1") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t1); - else if (strcmp(keyword, "t1'") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t1p); - else if (strcmp(keyword, "t2") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t2); - else if (strcmp(keyword, "t2'") == 0) StoreWithCheckDup_c(keyword, value, &StdI.t2p); - else if (strcmp(keyword, "t'") == 0) StoreWithCheckDup_c(keyword, value, &StdI.tp); - else if (strcmp(keyword, "t''") == 0) StoreWithCheckDup_c(keyword, value, &StdI.tpp); - else if (strcmp(keyword, "u") == 0) StoreWithCheckDup_d(keyword, value, &StdI.U); - else if (strcmp(keyword, "v") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V); - else if (strcmp(keyword, "v0") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V0); - else if (strcmp(keyword, "v0'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V0p); - else if (strcmp(keyword, "v1") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V1); - else if (strcmp(keyword, "v1'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V1p); - else if (strcmp(keyword, "v2") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V2); - else if (strcmp(keyword, "v2p") == 0) StoreWithCheckDup_d(keyword, value, &StdI.V2); - else if (strcmp(keyword, "v'") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Vp); - else if (strcmp(keyword, "v''") == 0) StoreWithCheckDup_d(keyword, value, &StdI.Vpp); - else if (strcmp(keyword, "w") == 0) StoreWithCheckDup_i(keyword, value, &StdI.W); - else if (strcmp(keyword, "wlength") == 0) StoreWithCheckDup_d(keyword, value, &StdI.length[0]); - else if (strcmp(keyword, "wx") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[0][0]); - else if (strcmp(keyword, "wy") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[0][1]); - else if (strcmp(keyword, "wz") == 0) StoreWithCheckDup_d(keyword, value, &StdI.direct[0][2]); - else if (strcmp(keyword, "w90_cutoff") == 0) StoreWithCheckDup_d(keyword, value, &StdI.W90_cutoff); - else if (strcmp(keyword, "w90_geom") == 0) StoreWithCheckDup_s(keyword, value, StdI.W90_geom); - else if (strcmp(keyword, "w90_hr") == 0) StoreWithCheckDup_s(keyword, value, StdI.W90_hr); - else if (strcmp(keyword, "2sz") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Sz2); + if (strcmp(keyword, "a") == 0) StoreWithCheckDup_d(keyword, value, &StdI->a); + else if (strcmp(keyword, "a0h") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[0][2]); + else if (strcmp(keyword, "a0l") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[0][1]); + else if (strcmp(keyword, "a0w") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[0][0]); + else if (strcmp(keyword, "a1h") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[1][2]); + else if (strcmp(keyword, "a1l") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[1][1]); + else if (strcmp(keyword, "a1w") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[1][0]); + else if (strcmp(keyword, "a2h") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[2][2]); + else if (strcmp(keyword, "a2l") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[2][1]); + else if (strcmp(keyword, "a2w") == 0) StoreWithCheckDup_i(keyword, value, &StdI->box[2][0]); + else if (strcmp(keyword, "d") == 0) StoreWithCheckDup_d(keyword, value, &StdI->D[2][2]); + else if (strcmp(keyword, "gamma") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Gamma); + else if (strcmp(keyword, "h") == 0) StoreWithCheckDup_d(keyword, value, &StdI->h); + else if (strcmp(keyword, "height") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Height); + else if (strcmp(keyword, "hlength") == 0) StoreWithCheckDup_d(keyword, value, &StdI->length[2]); + else if (strcmp(keyword, "hx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[2][0]); + else if (strcmp(keyword, "hy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[2][1]); + else if (strcmp(keyword, "hz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[2][2]); + else if (strcmp(keyword, "j") == 0) StoreWithCheckDup_d(keyword, value, &StdI->JAll); + else if (strcmp(keyword, "jx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[0][0]); + else if (strcmp(keyword, "jxy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[0][1]); + else if (strcmp(keyword, "jxz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[0][2]); + else if (strcmp(keyword, "jy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[1][1]); + else if (strcmp(keyword, "jyx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[1][0]); + else if (strcmp(keyword, "jyz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[1][2]); + else if (strcmp(keyword, "jz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[2][2]); + else if (strcmp(keyword, "jzx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[2][0]); + else if (strcmp(keyword, "jzy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J[2][1]); + else if (strcmp(keyword, "j0") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0All); + else if (strcmp(keyword, "j0x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[0][0]); + else if (strcmp(keyword, "j0xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[0][1]); + else if (strcmp(keyword, "j0xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[0][2]); + else if (strcmp(keyword, "j0y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[1][1]); + else if (strcmp(keyword, "j0yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[1][0]); + else if (strcmp(keyword, "j0yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[1][2]); + else if (strcmp(keyword, "j0z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[2][2]); + else if (strcmp(keyword, "j0zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[2][0]); + else if (strcmp(keyword, "j0zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0[2][1]); + else if (strcmp(keyword, "j0'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0pAll); + else if (strcmp(keyword, "j0'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[0][0]); + else if (strcmp(keyword, "j0'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[0][1]); + else if (strcmp(keyword, "j0'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[0][2]); + else if (strcmp(keyword, "j0'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[1][1]); + else if (strcmp(keyword, "j0'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[1][0]); + else if (strcmp(keyword, "j0'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[1][2]); + else if (strcmp(keyword, "j0'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[2][2]); + else if (strcmp(keyword, "j0'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[2][0]); + else if (strcmp(keyword, "j0'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J0p[2][1]); + else if (strcmp(keyword, "j1") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1All); + else if (strcmp(keyword, "j1x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[0][0]); + else if (strcmp(keyword, "j1xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[0][1]); + else if (strcmp(keyword, "j1xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[0][2]); + else if (strcmp(keyword, "j1y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[1][1]); + else if (strcmp(keyword, "j1yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[1][0]); + else if (strcmp(keyword, "j1yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[1][2]); + else if (strcmp(keyword, "j1z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[2][2]); + else if (strcmp(keyword, "j1zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[2][0]); + else if (strcmp(keyword, "j1zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1[2][1]); + else if (strcmp(keyword, "j1'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1pAll); + else if (strcmp(keyword, "j1'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[0][0]); + else if (strcmp(keyword, "j1'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[0][1]); + else if (strcmp(keyword, "j1'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[0][2]); + else if (strcmp(keyword, "j1'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[1][1]); + else if (strcmp(keyword, "j1'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[1][0]); + else if (strcmp(keyword, "j1'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[1][2]); + else if (strcmp(keyword, "j1'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[2][2]); + else if (strcmp(keyword, "j1'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[2][0]); + else if (strcmp(keyword, "j1'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J1p[2][1]); + else if (strcmp(keyword, "j2") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2All); + else if (strcmp(keyword, "j2x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[0][0]); + else if (strcmp(keyword, "j2xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[0][1]); + else if (strcmp(keyword, "j2xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[0][2]); + else if (strcmp(keyword, "j2y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[1][1]); + else if (strcmp(keyword, "j2yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[1][0]); + else if (strcmp(keyword, "j2yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[1][2]); + else if (strcmp(keyword, "j2z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[2][2]); + else if (strcmp(keyword, "j2zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[2][0]); + else if (strcmp(keyword, "j2zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2[2][1]); + else if (strcmp(keyword, "j2'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2pAll); + else if (strcmp(keyword, "j2'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[0][0]); + else if (strcmp(keyword, "j2'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[0][1]); + else if (strcmp(keyword, "j2'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[0][2]); + else if (strcmp(keyword, "j2'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[1][1]); + else if (strcmp(keyword, "j2'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[1][0]); + else if (strcmp(keyword, "j2'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[1][2]); + else if (strcmp(keyword, "j2'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[2][2]); + else if (strcmp(keyword, "j2'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[2][0]); + else if (strcmp(keyword, "j2'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->J2p[2][1]); + else if (strcmp(keyword, "j'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->JpAll); + else if (strcmp(keyword, "j'x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[0][0]); + else if (strcmp(keyword, "j'xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[0][1]); + else if (strcmp(keyword, "j'xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[0][2]); + else if (strcmp(keyword, "j'y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[1][1]); + else if (strcmp(keyword, "j'yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[1][0]); + else if (strcmp(keyword, "j'yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[1][2]); + else if (strcmp(keyword, "j'z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[2][2]); + else if (strcmp(keyword, "j'zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[2][0]); + else if (strcmp(keyword, "j'zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jp[2][1]); + else if (strcmp(keyword, "j''") == 0) StoreWithCheckDup_d(keyword, value, &StdI->JppAll); + else if (strcmp(keyword, "j''x") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[0][0]); + else if (strcmp(keyword, "j''xy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[0][1]); + else if (strcmp(keyword, "j''xz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[0][2]); + else if (strcmp(keyword, "j''y") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[1][1]); + else if (strcmp(keyword, "j''yx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[1][0]); + else if (strcmp(keyword, "j''yz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[1][2]); + else if (strcmp(keyword, "j''z") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[2][2]); + else if (strcmp(keyword, "j''zx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[2][0]); + else if (strcmp(keyword, "j''zy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Jpp[2][1]); + else if (strcmp(keyword, "k") == 0) StoreWithCheckDup_d(keyword, value, &StdI->K); + else if (strcmp(keyword, "l") == 0) StoreWithCheckDup_i(keyword, value, &StdI->L); + else if (strcmp(keyword, "lattice") == 0) StoreWithCheckDup_sl(keyword, value, StdI->lattice); + else if (strcmp(keyword, "llength") == 0) StoreWithCheckDup_d(keyword, value, &StdI->length[1]); + else if (strcmp(keyword, "lx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[1][0]); + else if (strcmp(keyword, "ly") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[1][1]); + else if (strcmp(keyword, "lz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[1][2]); + else if (strcmp(keyword, "model") == 0) StoreWithCheckDup_sl(keyword, value, StdI->model); + else if (strcmp(keyword, "mu") == 0) StoreWithCheckDup_d(keyword, value, &StdI->mu); + else if (strcmp(keyword, "nelec") == 0) StoreWithCheckDup_i(keyword, value, &StdI->nelec); + else if (strcmp(keyword, "outputmode") == 0) StoreWithCheckDup_sl(keyword, value, StdI->outputmode); + else if (strcmp(keyword, "phase0") == 0) StoreWithCheckDup_d(keyword, value, &StdI->phase[0]); + else if (strcmp(keyword, "phase1") == 0) StoreWithCheckDup_d(keyword, value, &StdI->phase[1]); + else if (strcmp(keyword, "phase2") == 0) StoreWithCheckDup_d(keyword, value, &StdI->phase[2]); + else if (strcmp(keyword, "t") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t); + else if (strcmp(keyword, "t0") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t0); + else if (strcmp(keyword, "t0'") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t0p); + else if (strcmp(keyword, "t1") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t1); + else if (strcmp(keyword, "t1'") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t1p); + else if (strcmp(keyword, "t2") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t2); + else if (strcmp(keyword, "t2'") == 0) StoreWithCheckDup_c(keyword, value, &StdI->t2p); + else if (strcmp(keyword, "t'") == 0) StoreWithCheckDup_c(keyword, value, &StdI->tp); + else if (strcmp(keyword, "t''") == 0) StoreWithCheckDup_c(keyword, value, &StdI->tpp); + else if (strcmp(keyword, "u") == 0) StoreWithCheckDup_d(keyword, value, &StdI->U); + else if (strcmp(keyword, "v") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V); + else if (strcmp(keyword, "v0") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V0); + else if (strcmp(keyword, "v0'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V0p); + else if (strcmp(keyword, "v1") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V1); + else if (strcmp(keyword, "v1'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V1p); + else if (strcmp(keyword, "v2") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V2); + else if (strcmp(keyword, "v2p") == 0) StoreWithCheckDup_d(keyword, value, &StdI->V2); + else if (strcmp(keyword, "v'") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Vp); + else if (strcmp(keyword, "v''") == 0) StoreWithCheckDup_d(keyword, value, &StdI->Vpp); + else if (strcmp(keyword, "w") == 0) StoreWithCheckDup_i(keyword, value, &StdI->W); + else if (strcmp(keyword, "wlength") == 0) StoreWithCheckDup_d(keyword, value, &StdI->length[0]); + else if (strcmp(keyword, "wx") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[0][0]); + else if (strcmp(keyword, "wy") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[0][1]); + else if (strcmp(keyword, "wz") == 0) StoreWithCheckDup_d(keyword, value, &StdI->direct[0][2]); + else if (strcmp(keyword, "w90_cutoff") == 0) StoreWithCheckDup_d(keyword, value, &StdI->W90_cutoff); + else if (strcmp(keyword, "w90_geom") == 0) StoreWithCheckDup_s(keyword, value, StdI->W90_geom); + else if (strcmp(keyword, "w90_hr") == 0) StoreWithCheckDup_s(keyword, value, StdI->W90_hr); + else if (strcmp(keyword, "2sz") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Sz2); #if defined(_HPhi) - else if (strcmp(keyword, "calcspec") == 0) StoreWithCheckDup_sl(keyword, value, StdI.CalcSpec); - else if (strcmp(keyword, "exct") == 0) StoreWithCheckDup_i(keyword, value, &StdI.exct); - else if (strcmp(keyword, "eigenvecio") == 0) StoreWithCheckDup_sl(keyword, value, StdI.EigenVecIO); - else if (strcmp(keyword, "expecinterval") == 0) StoreWithCheckDup_i(keyword, value, &StdI.ExpecInterval); - else if (strcmp(keyword, "cdatafilehead") == 0) StoreWithCheckDup_s(keyword, value, StdI.CDataFileHead); - else if (strcmp(keyword, "flgtemp") == 0) StoreWithCheckDup_i(keyword, value, &StdI.FlgTemp); - else if (strcmp(keyword, "initialvectype") == 0) StoreWithCheckDup_sl(keyword, value, StdI.InitialVecType); - else if (strcmp(keyword, "initial_iv") == 0) StoreWithCheckDup_i(keyword, value, &StdI.initial_iv); - else if (strcmp(keyword, "lanczoseps") == 0) StoreWithCheckDup_i(keyword, value, &StdI.LanczosEps); - else if (strcmp(keyword, "lanczostarget") == 0) StoreWithCheckDup_i(keyword, value, &StdI.LanczosTarget); - else if (strcmp(keyword, "lanczos_max") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Lanczos_max); - else if (strcmp(keyword, "largevalue") == 0) StoreWithCheckDup_d(keyword, value, &StdI.LargeValue); - else if (strcmp(keyword, "method") == 0) StoreWithCheckDup_sl(keyword, value, StdI.method); - else if (strcmp(keyword, "nomega") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Nomega); - else if (strcmp(keyword, "numave") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NumAve); - else if (strcmp(keyword, "nvec") == 0) StoreWithCheckDup_i(keyword, value, &StdI.nvec); - else if (strcmp(keyword, "omegamax") == 0) StoreWithCheckDup_d(keyword, value, &StdI.OmegaMax); - else if (strcmp(keyword, "omegamin") == 0) StoreWithCheckDup_d(keyword, value, &StdI.OmegaMin); - else if (strcmp(keyword, "omegaim") == 0) StoreWithCheckDup_d(keyword, value, &StdI.OmegaIm); - else if (strcmp(keyword, "restart") == 0) StoreWithCheckDup_sl(keyword, value, StdI.Restart); - else if (strcmp(keyword, "spectrumqh") == 0) StoreWithCheckDup_d(keyword, value, &StdI.SpectrumQH); - else if (strcmp(keyword, "spectrumql") == 0) StoreWithCheckDup_d(keyword, value, &StdI.SpectrumQL); - else if (strcmp(keyword, "spectrumqw") == 0) StoreWithCheckDup_d(keyword, value, &StdI.SpectrumQW); - else if (strcmp(keyword, "spectrumtype") == 0) StoreWithCheckDup_sl(keyword, value, StdI.SpectrumType); - else if (strcmp(keyword, "2s") == 0) StoreWithCheckDup_i(keyword, value, &StdI.S2); + else if (strcmp(keyword, "calcspec") == 0) StoreWithCheckDup_sl(keyword, value, StdI->CalcSpec); + else if (strcmp(keyword, "exct") == 0) StoreWithCheckDup_i(keyword, value, &StdI->exct); + else if (strcmp(keyword, "eigenvecio") == 0) StoreWithCheckDup_sl(keyword, value, StdI->EigenVecIO); + else if (strcmp(keyword, "expecinterval") == 0) StoreWithCheckDup_i(keyword, value, &StdI->ExpecInterval); + else if (strcmp(keyword, "cdatafilehead") == 0) StoreWithCheckDup_s(keyword, value, StdI->CDataFileHead); + else if (strcmp(keyword, "flgtemp") == 0) StoreWithCheckDup_i(keyword, value, &StdI->FlgTemp); + else if (strcmp(keyword, "initialvectype") == 0) StoreWithCheckDup_sl(keyword, value, StdI->InitialVecType); + else if (strcmp(keyword, "initial_iv") == 0) StoreWithCheckDup_i(keyword, value, &StdI->initial_iv); + else if (strcmp(keyword, "lanczoseps") == 0) StoreWithCheckDup_i(keyword, value, &StdI->LanczosEps); + else if (strcmp(keyword, "lanczostarget") == 0) StoreWithCheckDup_i(keyword, value, &StdI->LanczosTarget); + else if (strcmp(keyword, "lanczos_max") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Lanczos_max); + else if (strcmp(keyword, "largevalue") == 0) StoreWithCheckDup_d(keyword, value, &StdI->LargeValue); + else if (strcmp(keyword, "method") == 0) StoreWithCheckDup_sl(keyword, value, StdI->method); + else if (strcmp(keyword, "nomega") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Nomega); + else if (strcmp(keyword, "numave") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NumAve); + else if (strcmp(keyword, "nvec") == 0) StoreWithCheckDup_i(keyword, value, &StdI->nvec); + else if (strcmp(keyword, "omegamax") == 0) StoreWithCheckDup_d(keyword, value, &StdI->OmegaMax); + else if (strcmp(keyword, "omegamin") == 0) StoreWithCheckDup_d(keyword, value, &StdI->OmegaMin); + else if (strcmp(keyword, "omegaim") == 0) StoreWithCheckDup_d(keyword, value, &StdI->OmegaIm); + else if (strcmp(keyword, "restart") == 0) StoreWithCheckDup_sl(keyword, value, StdI->Restart); + else if (strcmp(keyword, "spectrumqh") == 0) StoreWithCheckDup_d(keyword, value, &StdI->SpectrumQH); + else if (strcmp(keyword, "spectrumql") == 0) StoreWithCheckDup_d(keyword, value, &StdI->SpectrumQL); + else if (strcmp(keyword, "spectrumqw") == 0) StoreWithCheckDup_d(keyword, value, &StdI->SpectrumQW); + else if (strcmp(keyword, "spectrumtype") == 0) StoreWithCheckDup_sl(keyword, value, StdI->SpectrumType); + else if (strcmp(keyword, "2s") == 0) StoreWithCheckDup_i(keyword, value, &StdI->S2); #elif defined(_mVMC) - else if (strcmp(keyword, "a0hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[0][2]); - else if (strcmp(keyword, "a0lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[0][1]); - else if (strcmp(keyword, "a0wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[0][0]); - else if (strcmp(keyword, "a1hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[1][2]); - else if (strcmp(keyword, "a1lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[1][1]); - else if (strcmp(keyword, "a1wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[1][0]); - else if (strcmp(keyword, "a2hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[2][2]); - else if (strcmp(keyword, "a2lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[2][1]); - else if (strcmp(keyword, "a2wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.boxsub[2][0]); - else if (strcmp(keyword, "complextype") == 0) StoreWithCheckDup_i(keyword, value, &StdI.ComplexType); - else if (strcmp(keyword, "cparafilehead") == 0) StoreWithCheckDup_s(keyword, value, StdI.CParaFileHead); - else if (strcmp(keyword, "dsroptredcut") == 0) StoreWithCheckDup_d(keyword, value, &StdI.DSROptRedCut); - else if (strcmp(keyword, "dsroptstadel") == 0) StoreWithCheckDup_d(keyword, value, &StdI.DSROptStaDel); - else if (strcmp(keyword, "dsroptstepdt") == 0) StoreWithCheckDup_d(keyword, value, &StdI.DSROptStepDt); - else if (strcmp(keyword, "hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Hsub); - else if (strcmp(keyword, "lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Lsub); - else if (strcmp(keyword, "nvmccalmode") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NVMCCalMode); - else if (strcmp(keyword, "ndataidxstart") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NDataIdxStart); - else if (strcmp(keyword, "ndataqtysmp") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NDataQtySmp); - else if (strcmp(keyword, "nlanczosmode") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NLanczosMode); - else if (strcmp(keyword, "nmptrans") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NMPTrans); - else if (strcmp(keyword, "nspgaussleg") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NSPGaussLeg); - else if (strcmp(keyword, "nsplitsize") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NSplitSize); - else if (strcmp(keyword, "nspstot") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NSPStot); - else if (strcmp(keyword, "nsroptitrsmp") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NSROptItrSmp); - else if (strcmp(keyword, "nsroptitrstep") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NSROptItrStep); - else if (strcmp(keyword, "nstore") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NStore); - else if (strcmp(keyword, "nsrcg") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NSRCG); - else if (strcmp(keyword, "nvmcinterval") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NVMCInterval); - else if (strcmp(keyword, "nvmcsample") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NVMCSample); - else if (strcmp(keyword, "nvmcwarmup") == 0) StoreWithCheckDup_i(keyword, value, &StdI.NVMCWarmUp); - else if (strcmp(keyword, "rndseed") == 0) StoreWithCheckDup_i(keyword, value, &StdI.RndSeed); - else if (strcmp(keyword, "wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI.Wsub); + else if (strcmp(keyword, "a0hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[0][2]); + else if (strcmp(keyword, "a0lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[0][1]); + else if (strcmp(keyword, "a0wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[0][0]); + else if (strcmp(keyword, "a1hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[1][2]); + else if (strcmp(keyword, "a1lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[1][1]); + else if (strcmp(keyword, "a1wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[1][0]); + else if (strcmp(keyword, "a2hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[2][2]); + else if (strcmp(keyword, "a2lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[2][1]); + else if (strcmp(keyword, "a2wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->boxsub[2][0]); + else if (strcmp(keyword, "complextype") == 0) StoreWithCheckDup_i(keyword, value, &StdI->ComplexType); + else if (strcmp(keyword, "cparafilehead") == 0) StoreWithCheckDup_s(keyword, value, StdI->CParaFileHead); + else if (strcmp(keyword, "dsroptredcut") == 0) StoreWithCheckDup_d(keyword, value, &StdI->DSROptRedCut); + else if (strcmp(keyword, "dsroptstadel") == 0) StoreWithCheckDup_d(keyword, value, &StdI->DSROptStaDel); + else if (strcmp(keyword, "dsroptstepdt") == 0) StoreWithCheckDup_d(keyword, value, &StdI->DSROptStepDt); + else if (strcmp(keyword, "hsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Hsub); + else if (strcmp(keyword, "lsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Lsub); + else if (strcmp(keyword, "nvmccalmode") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NVMCCalMode); + else if (strcmp(keyword, "ndataidxstart") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NDataIdxStart); + else if (strcmp(keyword, "ndataqtysmp") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NDataQtySmp); + else if (strcmp(keyword, "nlanczosmode") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NLanczosMode); + else if (strcmp(keyword, "nmptrans") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NMPTrans); + else if (strcmp(keyword, "nspgaussleg") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NSPGaussLeg); + else if (strcmp(keyword, "nsplitsize") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NSplitSize); + else if (strcmp(keyword, "nspstot") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NSPStot); + else if (strcmp(keyword, "nsroptitrsmp") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NSROptItrSmp); + else if (strcmp(keyword, "nsroptitrstep") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NSROptItrStep); + else if (strcmp(keyword, "nstore") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NStore); + else if (strcmp(keyword, "nsrcg") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NSRCG); + else if (strcmp(keyword, "nvmcinterval") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NVMCInterval); + else if (strcmp(keyword, "nvmcsample") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NVMCSample); + else if (strcmp(keyword, "nvmcwarmup") == 0) StoreWithCheckDup_i(keyword, value, &StdI->NVMCWarmUp); + else if (strcmp(keyword, "rndseed") == 0) StoreWithCheckDup_i(keyword, value, &StdI->RndSeed); + else if (strcmp(keyword, "wsub") == 0) StoreWithCheckDup_i(keyword, value, &StdI->Wsub); #endif else { fprintf(stdout, "ERROR ! Unsupported Keyword !\n"); @@ -2048,129 +2059,190 @@ void StdFace_main( /* Check the model */ - StdI.lGC = 0; - StdI.lBoost = 0; - if (strcmp(StdI.model, "fermionhubbard") == 0 - || strcmp(StdI.model, "hubbard") == 0) - strcpy(StdI.model, "hubbard\0"); - else if(strcmp(StdI.model, "fermionhubbardgc") == 0 - || strcmp(StdI.model, "hubbardgc") == 0) { - strcpy(StdI.model, "hubbard\0"); - StdI.lGC = 1; + StdI->lGC = 0; + StdI->lBoost = 0; + if (strcmp(StdI->model, "fermionhubbard") == 0 + || strcmp(StdI->model, "hubbard") == 0) + strcpy(StdI->model, "hubbard\0"); + else if(strcmp(StdI->model, "fermionhubbardgc") == 0 + || strcmp(StdI->model, "hubbardgc") == 0) { + strcpy(StdI->model, "hubbard\0"); + StdI->lGC = 1; } - else if (strcmp(StdI.model, "spin") == 0) - strcpy(StdI.model, "spin\0"); - else if (strcmp(StdI.model, "spingc") == 0) { - strcpy(StdI.model, "spin\0"); - StdI.lGC = 1; + else if (strcmp(StdI->model, "spin") == 0) + strcpy(StdI->model, "spin\0"); + else if (strcmp(StdI->model, "spingc") == 0) { + strcpy(StdI->model, "spin\0"); + StdI->lGC = 1; } #if defined(_HPhi) - else if(strcmp(StdI.model, "spingcboost") == 0 || - strcmp(StdI.model, "spingccma") == 0) { - strcpy(StdI.model, "spin\0"); - StdI.lGC = 1; - StdI.lBoost = 1; + else if(strcmp(StdI->model, "spingcboost") == 0 || + strcmp(StdI->model, "spingccma") == 0) { + strcpy(StdI->model, "spin\0"); + StdI->lGC = 1; + StdI->lBoost = 1; } #endif - else if (strcmp(StdI.model, "kondolattice") == 0 - || strcmp(StdI.model, "kondo") == 0) { - strcpy(StdI.model, "kondo\0"); + else if (strcmp(StdI->model, "kondolattice") == 0 + || strcmp(StdI->model, "kondo") == 0) { + strcpy(StdI->model, "kondo\0"); } - else if(strcmp(StdI.model, "kondolatticegc") == 0 - || strcmp(StdI.model, "kondogc") == 0) { - strcpy(StdI.model, "kondo\0"); - StdI.lGC = 1; + else if(strcmp(StdI->model, "kondolatticegc") == 0 + || strcmp(StdI->model, "kondogc") == 0) { + strcpy(StdI->model, "kondo\0"); + StdI->lGC = 1; } - else UnsupportedSystem(StdI.model, StdI.lattice); + else UnsupportedSystem(StdI->model, StdI->lattice); - /* + /*>> Generate Hamiltonian definition files */ - if (strcmp(StdI.lattice, "chain") == 0 - || strcmp(StdI.lattice, "chainlattice") == 0) StdFace_Chain(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "facecenteredorthorhombic") == 0 - || strcmp(StdI.lattice, "fcorthorhombic") == 0 - || strcmp(StdI.lattice, "fco") == 0) StdFace_FCOrtho(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "honeycomb") == 0 - || strcmp(StdI.lattice, "honeycomblattice") == 0) StdFace_Honeycomb(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "kagome") == 0 - || strcmp(StdI.lattice, "kagomelattice") == 0) StdFace_Kagome(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "ladder") == 0 - || strcmp(StdI.lattice, "ladderlattice") == 0) StdFace_Ladder(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "orthorhombic") == 0 - || strcmp(StdI.lattice, "simpleorthorhombic") == 0) StdFace_Orthorhombic(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "pyrochlore") == 0) StdFace_Pyrochlore(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "tetragonal") == 0 - || strcmp(StdI.lattice, "tetragonallattice") == 0 - || strcmp(StdI.lattice, "square") == 0 - || strcmp(StdI.lattice, "squarelattice") == 0) StdFace_Tetragonal(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "triangular") == 0 - || strcmp(StdI.lattice, "triangularlattice") == 0) StdFace_Triangular(&StdI, StdI.model); - else if (strcmp(StdI.lattice, "wannier90") == 0) StdFace_Wannier90(&StdI, StdI.model); - else UnsupportedSystem(StdI.model, StdI.lattice); + if (strcmp(StdI->lattice, "chain") == 0 + || strcmp(StdI->lattice, "chainlattice") == 0) StdFace_Chain(StdI); + else if (strcmp(StdI->lattice, "facecenteredorthorhombic") == 0 + || strcmp(StdI->lattice, "fcorthorhombic") == 0 + || strcmp(StdI->lattice, "fco") == 0) StdFace_FCOrtho(StdI); + else if (strcmp(StdI->lattice, "honeycomb") == 0 + || strcmp(StdI->lattice, "honeycomblattice") == 0) StdFace_Honeycomb(StdI); + else if (strcmp(StdI->lattice, "kagome") == 0 + || strcmp(StdI->lattice, "kagomelattice") == 0) StdFace_Kagome(StdI); + else if (strcmp(StdI->lattice, "ladder") == 0 + || strcmp(StdI->lattice, "ladderlattice") == 0) StdFace_Ladder(StdI); + else if (strcmp(StdI->lattice, "orthorhombic") == 0 + || strcmp(StdI->lattice, "simpleorthorhombic") == 0) StdFace_Orthorhombic(StdI); + else if (strcmp(StdI->lattice, "pyrochlore") == 0) StdFace_Pyrochlore(StdI); + else if (strcmp(StdI->lattice, "tetragonal") == 0 + || strcmp(StdI->lattice, "tetragonallattice") == 0 + || strcmp(StdI->lattice, "square") == 0 + || strcmp(StdI->lattice, "squarelattice") == 0) StdFace_Tetragonal(StdI); + else if (strcmp(StdI->lattice, "triangular") == 0 + || strcmp(StdI->lattice, "triangularlattice") == 0) StdFace_Triangular(StdI); + else if (strcmp(StdI->lattice, "wannier90") == 0) StdFace_Wannier90(StdI); + else UnsupportedSystem(StdI->model, StdI->lattice);//<< /**/ #if defined(_HPhi) - StdFace_LargeValue(&StdI); + StdFace_LargeValue(StdI); /* Generate Hamiltonian for Boost */ - if (StdI.lBoost == 1) { - if (strcmp(StdI.lattice, "chain") == 0 - || strcmp(StdI.lattice, "chainlattice") == 0) StdFace_Chain_Boost(&StdI); - else if (strcmp(StdI.lattice, "honeycomb") == 0 - || strcmp(StdI.lattice, "honeycomblattice") == 0) StdFace_Honeycomb_Boost(&StdI); - else if (strcmp(StdI.lattice, "kagome") == 0 - || strcmp(StdI.lattice, "kagomelattice") == 0) StdFace_Kagome_Boost(&StdI); - else if (strcmp(StdI.lattice, "ladder") == 0 - || strcmp(StdI.lattice, "ladderlattice") == 0) StdFace_Ladder_Boost(&StdI); - else UnsupportedSystem(StdI.model, StdI.lattice); + if (StdI->lBoost == 1) { + if (strcmp(StdI->lattice, "chain") == 0 + || strcmp(StdI->lattice, "chainlattice") == 0) StdFace_Chain_Boost(StdI); + else if (strcmp(StdI->lattice, "honeycomb") == 0 + || strcmp(StdI->lattice, "honeycomblattice") == 0) StdFace_Honeycomb_Boost(StdI); + else if (strcmp(StdI->lattice, "kagome") == 0 + || strcmp(StdI->lattice, "kagomelattice") == 0) StdFace_Kagome_Boost(StdI); + else if (strcmp(StdI->lattice, "ladder") == 0 + || strcmp(StdI->lattice, "ladderlattice") == 0) StdFace_Ladder_Boost(StdI); + else UnsupportedSystem(StdI->model, StdI->lattice); } #endif /**/ fprintf(stdout, "\n"); fprintf(stdout, "###### Print Expert input files ######\n"); fprintf(stdout, "\n"); - PrintLocSpin(&StdI); - PrintTrans(&StdI); - PrintInteractions(&StdI); + PrintLocSpin(StdI); + PrintTrans(StdI); + PrintInteractions(StdI); #if defined(_HPhi) - PrintExcitation(&StdI); - PrintCalcMod(&StdI); + PrintExcitation(StdI); + PrintCalcMod(StdI); #elif defined(_mVMC) - if(StdI.lGC == 0 && (StdI.Sz2 == 0 || StdI.Sz2 == StdI.NaN_i)) - StdFace_PrintVal_i("ComplexType", &StdI.ComplexType, 0); - else StdFace_PrintVal_i("ComplexType", &StdI.ComplexType, 1); - - StdFace_generate_orb(&StdI); - StdFace_Proj(&StdI); - PrintJastrow(&StdI); - if(StdI.lGC == 1 || (StdI.Sz2 != 0 && StdI.Sz2 != StdI.NaN_i) ) - PrintOrbPara(&StdI); - PrintGutzwiller(&StdI); - PrintOrb(&StdI); + if(StdI->lGC == 0 && (StdI->Sz2 == 0 || StdI->Sz2 == StdI->NaN_i)) + StdFace_PrintVal_i("ComplexType", &StdI->ComplexType, 0); + else StdFace_PrintVal_i("ComplexType", &StdI->ComplexType, 1); + + StdFace_generate_orb(StdI); + StdFace_Proj(StdI); + PrintJastrow(StdI); + if(StdI->lGC == 1 || (StdI->Sz2 != 0 && StdI->Sz2 != StdI->NaN_i) ) + PrintOrbPara(StdI); + PrintGutzwiller(StdI); + PrintOrb(StdI); #endif - CheckModPara(&StdI); - PrintModPara(&StdI); - CheckOutputMode(&StdI); - Print1Green(&StdI); - Print2Green(&StdI); - PrintNamelist(&StdI); + CheckModPara(StdI); + PrintModPara(StdI); + CheckOutputMode(StdI); + Print1Green(StdI); + Print2Green(StdI); + PrintNamelist(StdI); /* Finalize All */ - free(StdI.locspinflag); - for (ktrans = 0; ktrans < StdI.ntrans; ktrans++) { - free(StdI.transindx[ktrans]); + free(StdI->locspinflag); + for (ktrans = 0; ktrans < StdI->ntrans; ktrans++) { + free(StdI->transindx[ktrans]); } - free(StdI.transindx); - free(StdI.trans); - for (kintr = 0; kintr < StdI.nintr; kintr++) { - free(StdI.intrindx[kintr]); + free(StdI->transindx); + free(StdI->trans); + for (kintr = 0; kintr < StdI->nintr; kintr++) { + free(StdI->intrindx[kintr]); } - free(StdI.intrindx); - free(StdI.intr); + free(StdI->intrindx); + free(StdI->intr); fprintf(stdout, "\n###### Input files are generated. ######\n\n"); }/*void StdFace_main*/ +/** +@page page_addstandard Add new lattice model into Standard mode + +@section sec_stan_proc Overall procedure + +If you want to create new lattice file, do as these files. + +-# Copy one of laattice files such as Kagome.c + (Probably the most similar one) and rename it. +-# @ref sec_lattice +-# Add that function in the header file, StdFace_ModelUtil.h +-# Add entry at + @dontinclude StdFace_main.c + @skip StdFace\_main + @until StdIntList + : + @skip >> + @until << +. +
+@section sec_lattice Modify lattice model file + +To create new lattice file, please modify the following part +(Kagome.c as an example): + +@dontinclude Kagome.c +Define function as +@skip StdFace\_Kagome( +@until { +Lattice parameter used only in geometry.dat and lattice.gp +@skip StdFace\_PrintVal\_d +@until Ly +these are unit lattice vectors.\n +Just call this function to initialize all lattice related parameters +@skipline StdFace\_InitSite +where "2" indicates 2D +@skip tau +@until tau\[2\]\[0\] +These are the fractional coordinate of internal sites. +Then set parameters of Hamiltonian +@skip StdFace\_NotUsed\_J +@until @@ +to determine the default value of them and unused parameters. +For more details, please see the description of each functions. +Then Compute the upper limit of the number of Transfer & Interaction and malloc them. +@skip >> +@until << +Please estimate the number of bonds per site. +@skipline kCell +In this loop, the parameters of Hamiltonian are computed & stored. +The local term is computed as follows: +@skip >> +@until << +Probably, it is not necessary to modify this part. +The non-local term is as follows: +@skip >> +@until << +For more details, please see each functions. + +StdFace_Kagome_Boost()? Forget!! +*/ \ No newline at end of file diff --git a/src/StdFace/StdFace_vals.h b/src/StdFace/StdFace_vals.h index 0b2a995e..ae5dd08a 100644 --- a/src/StdFace/StdFace_vals.h +++ b/src/StdFace/StdFace_vals.h @@ -136,176 +136,176 @@ struct StdIntList { /* Transfer, Interaction, Locspin */ - int nsite;/**@brief Number of sites, set in the each lattice file.*/ - int *locspinflag;/**@brief [StdIntList::nsite] LocSpin in Expert mode, + int nsite;/**<@brief Number of sites, set in the each lattice file.*/ + int *locspinflag;/**<@brief [StdIntList::nsite] LocSpin in Expert mode, malloc and set in each lattice file.*/ - int ntrans;/**@brief Number of transfer, counted in each lattice file.*/ - int Ltrans;/**@brief Print trans.def or not, set in PrintTrans().*/ - int **transindx;/**@brief [StdIntList::ntrans][4] Site/spin indices of + int ntrans;/**<@brief Number of transfer, counted in each lattice file.*/ + int Ltrans;/**<@brief Print trans.def or not, set in PrintTrans().*/ + int **transindx;/**<@brief [StdIntList::ntrans][4] Site/spin indices of one-body term, malloc in StdFace_MallocInteractions() and set in StdFace_trans().*/ - double complex *trans;/**@brief [StdIntList::ntrans] Coefficient of + double complex *trans;/**<@brief [StdIntList::ntrans] Coefficient of one-body term, malloc in StdFace_MallocInteractions() and set in StdFace_trans().*/ - int nintr;/**@brief Number of InterAll, counted in each lattice file.*/ - int Lintr;/**@brief Print interall.def or not, set in PrintInteractions().*/ - int **intrindx;/**@brief [StdIntList::nintr][8] Site/spin indices of + int nintr;/**<@brief Number of InterAll, counted in each lattice file.*/ + int Lintr;/**<@brief Print interall.def or not, set in PrintInteractions().*/ + int **intrindx;/**<@brief [StdIntList::nintr][8] Site/spin indices of two-body term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - double complex *intr;/**@brief [StdIntList::nintr] Coefficient of general + double complex *intr;/**<@brief [StdIntList::nintr] Coefficient of general two-body term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - int NCintra;/**@brief Number of intra-site Coulomb interaction, + int NCintra;/**<@brief Number of intra-site Coulomb interaction, counted in each lattice file.*/ - int LCintra;/**@brief Print coulombintra.def or not, set in PrintInteractions().*/ - int **CintraIndx;/**@brief [StdIntList::NCintra][1] Site indices of + int LCintra;/**<@brief Print coulombintra.def or not, set in PrintInteractions().*/ + int **CintraIndx;/**<@brief [StdIntList::NCintra][1] Site indices of intra-site Coulomb term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - double *Cintra;/**@brief [StdIntList::NCintra] Coefficient of intra-site + double *Cintra;/**<@brief [StdIntList::NCintra] Coefficient of intra-site Coulomb term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - int NCinter;/**@brief Number of inter-site Coulomb interaction, + int NCinter;/**<@brief Number of inter-site Coulomb interaction, counted in each lattice file.*/ - int LCinter;/**@brief Print coulombinter.def or not, set in PrintInteractions().*/ - int **CinterIndx;/**@brief [StdIntList::NCinter][2] Site indices of + int LCinter;/**<@brief Print coulombinter.def or not, set in PrintInteractions().*/ + int **CinterIndx;/**<@brief [StdIntList::NCinter][2] Site indices of inter-site Coulomb term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - double *Cinter;/**@brief [StdIntList::NCinter] Coefficient of inter-site + double *Cinter;/**<@brief [StdIntList::NCinter] Coefficient of inter-site Coulomb term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - int NHund;/**@brief Number of Hund term, counted in each lattice file.*/ - int LHund;/**@brief Print hund.def or not, set in PrintInteractions().*/ - int **HundIndx;/**@brief [StdIntList::NHund][2] Site indices of + int NHund;/**<@brief Number of Hund term, counted in each lattice file.*/ + int LHund;/**<@brief Print hund.def or not, set in PrintInteractions().*/ + int **HundIndx;/**<@brief [StdIntList::NHund][2] Site indices of Hund term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - double *Hund;/**@brief [StdIntList::NHund] Coefficient of Hund term, + double *Hund;/**<@brief [StdIntList::NHund] Coefficient of Hund term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - int NEx;/**@brief Number of exchange term, counted in each lattice file.*/ - int LEx;/**@brief Print exchange.def or not, set in PrintInteractions().*/ - int **ExIndx;/**@brief [StdIntList::NEx][2] Site indices of + int NEx;/**<@brief Number of exchange term, counted in each lattice file.*/ + int LEx;/**<@brief Print exchange.def or not, set in PrintInteractions().*/ + int **ExIndx;/**<@brief [StdIntList::NEx][2] Site indices of exchange term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - double *Ex;/**@brief [StdIntList::NEx] Coefficient of exchange term, + double *Ex;/**<@brief [StdIntList::NEx] Coefficient of exchange term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - int NPairLift;/**@brief Number of pair-lift term, counted in each lattice file.*/ - int LPairLift;/**@brief Print pairlift.def or not, set in PrintInteractions().*/ - int **PLIndx;/**@brief [StdIntList::NPairLift][2] Site indices of + int NPairLift;/**<@brief Number of pair-lift term, counted in each lattice file.*/ + int LPairLift;/**<@brief Print pairlift.def or not, set in PrintInteractions().*/ + int **PLIndx;/**<@brief [StdIntList::NPairLift][2] Site indices of pair-lift term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ - double *PairLift;/**@brief [StdIntList::NPairLift] Coefficient of + double *PairLift;/**<@brief [StdIntList::NPairLift] Coefficient of pair-lift term, malloc in StdFace_MallocInteractions() and set in StdFace_intr().*/ int lBoost; /* Calculation conditions */ - int lGC;/**@brief Switch for computing Grandcanonical ensemble(== 1). + int lGC;/**<@brief Switch for computing Grandcanonical ensemble(== 1). Setted in StdFace_main() after all keywords are read.*/ - int nelec;/**@brief Number of electrons, input from file.*/ - int S2;/**@brief Total spin |S| of a local spin, input from file.*/ - char outputmode[256];/**@brief Select amount of correlation function, + int nelec;/**<@brief Number of electrons, input from file.*/ + int S2;/**<@brief Total spin |S| of a local spin, input from file.*/ + char outputmode[256];/**<@brief Select amount of correlation function, input from file.*/ - char CDataFileHead[256];/**@brief String fron tof the output files, + char CDataFileHead[256];/**<@brief String fron tof the output files, input from file*/ - int Sz2;/**@brief Total Sz, input from file.*/ - int ioutputmode;/**@brief Switch associated to StdIntList::outputmode*/ + int Sz2;/**<@brief Total Sz, input from file.*/ + int ioutputmode;/**<@brief Switch associated to StdIntList::outputmode*/ /* Wannier90 mode */ - char W90_hr[256];/**@brief Name of hopping parameter file from wannier90, + char W90_hr[256];/**<@brief Name of hopping parameter file from wannier90, input from file.*/ - char W90_geom[256];/**@brief Name of geometry file from wannier90 converter, + char W90_geom[256];/**<@brief Name of geometry file from wannier90 converter, input from file.*/ - int W90_nt;/**@brief Number of transfer in wannier90 HR file.*/ - int **W90_indx;/**@brief [StdIntList::W90_nt][5] Hopping index, + int W90_nt;/**<@brief Number of transfer in wannier90 HR file.*/ + int **W90_indx;/**<@brief [StdIntList::W90_nt][5] Hopping index, malloc in read_W90().*/ - double complex *W90_t;/**@brief [StdIntList::W90_nt] Hopping parameter, + double complex *W90_t;/**<@brief [StdIntList::W90_nt] Hopping parameter, malloc in read_W90().*/ - double W90_cutoff;/**@brief Cutoof for the hopping in wannier90, input from file*/ + double W90_cutoff;/**<@brief Cutoof for the hopping in wannier90, input from file*/ #if defined(_HPhi) /* HPhi modpara */ - char method[256];/**@brief The name of method, input from file.*/ - char Restart[256];/**@brief The name of restart mode, input from file.*/ - char InitialVecType[256];/**@brief The name of initialguess-type, input from file.*/ - char EigenVecIO[256];/**@brief The name of I/O mode for eigenvector, input from file*/ - int FlgTemp;/**@brief */ - int Lanczos_max;/**@brief The maxixmum number of iterations, input from file*/ - int initial_iv; /**@brief the number for generating random number, input from file.*/ - int nvec;/**@brief */ - int exct;/**@brief The number of eigenvectors to be computed. input from file*/ - int LanczosEps;/**@brief Convergence threshold for the Lanczos method.*/ - int LanczosTarget;/**@brief Which eigenvector is used for the convergence check.*/ - int NumAve;/**@brief Number of trials for TPQ calculation.*/ - int ExpecInterval;/**@brief Interval for the iteration when the expectation + char method[256];/**<@brief The name of method, input from file.*/ + char Restart[256];/**<@brief The name of restart mode, input from file.*/ + char InitialVecType[256];/**<@brief The name of initialguess-type, input from file.*/ + char EigenVecIO[256];/**<@brief The name of I/O mode for eigenvector, input from file*/ + int FlgTemp;/**<@brief */ + int Lanczos_max;/**<@brief The maxixmum number of iterations, input from file*/ + int initial_iv; /**<@brief the number for generating random number, input from file.*/ + int nvec;/**<@brief */ + int exct;/**<@brief The number of eigenvectors to be computed. input from file*/ + int LanczosEps;/**<@brief Convergence threshold for the Lanczos method.*/ + int LanczosTarget;/**<@brief Which eigenvector is used for the convergence check.*/ + int NumAve;/**<@brief Number of trials for TPQ calculation.*/ + int ExpecInterval;/**<@brief Interval for the iteration when the expectation value is computed.*/ - double LargeValue;/**@brief The shift parameter for the TPQ calculation.*/ + double LargeValue;/**<@brief The shift parameter for the TPQ calculation.*/ /* Boost */ - int ***list_6spin_pair;/**@brief */ - int **list_6spin_star;/**@brief */ - int num_pivot;/**@brief */ - int ishift_nspin;/**@brief */ + int ***list_6spin_pair;/**<@brief */ + int **list_6spin_star;/**<@brief */ + int num_pivot;/**<@brief */ + int ishift_nspin;/**<@brief */ /*Spectrum*/ - char CalcSpec[256];/**@brief The name of mode for spectrum, input from file.*/ - char SpectrumType[256];/**@brief The type of mode for spectrum, input from file.*/ - int Nomega;/**@brief Number of frequencies, input from file.*/ - double OmegaMax;/**@brief Maximum of frequency for spectrum, input from file.*/ - double OmegaMin;/**@brief Minimum of frequency for spectrum, input from file.*/ - double OmegaIm;/**@brief Imaginary part of frequency.*/ - double SpectrumQL;/**@brief wavenumver (q-vector) in fractional coordinate*/ - double SpectrumQW;/**@brief wavenumver (q-vector) in fractional coordinate*/ - double SpectrumQH;/**@brief wavenumver (q-vector) in fractional coordinate*/ - int SpectrumBody;/**@brief one- or two-body excitation, defined from + char CalcSpec[256];/**<@brief The name of mode for spectrum, input from file.*/ + char SpectrumType[256];/**<@brief The type of mode for spectrum, input from file.*/ + int Nomega;/**<@brief Number of frequencies, input from file.*/ + double OmegaMax;/**<@brief Maximum of frequency for spectrum, input from file.*/ + double OmegaMin;/**<@brief Minimum of frequency for spectrum, input from file.*/ + double OmegaIm;/**<@brief Imaginary part of frequency.*/ + double SpectrumQL;/**<@brief wavenumver (q-vector) in fractional coordinate*/ + double SpectrumQW;/**<@brief wavenumver (q-vector) in fractional coordinate*/ + double SpectrumQH;/**<@brief wavenumver (q-vector) in fractional coordinate*/ + int SpectrumBody;/**<@brief one- or two-body excitation, defined from StdIntList::SpectrumType*/ #elif defined(_mVMC) /*mVMC modpara*/ - char CParaFileHead[256];/**@brief Header of the optimized wavefunction, + char CParaFileHead[256];/**<@brief Header of the optimized wavefunction, input from file*/ - int NVMCCalMode;/**@brief Optimization(=0) or compute correlation + int NVMCCalMode;/**<@brief Optimization(=0) or compute correlation function(=1), input from file.*/ - int NLanczosMode;/**@brief Power Lanczos(=1), input from file*/ - int NDataIdxStart;/**@brief Start index of trials, input from file.*/ - int NDataQtySmp;/**@brief Number of trials, input from file.*/ - int NSPGaussLeg;/**@brief Number of Gauss-Legendre points for spin projection, + int NLanczosMode;/**<@brief Power Lanczos(=1), input from file*/ + int NDataIdxStart;/**<@brief Start index of trials, input from file.*/ + int NDataQtySmp;/**<@brief Number of trials, input from file.*/ + int NSPGaussLeg;/**<@brief Number of Gauss-Legendre points for spin projection, input from file.*/ - int NMPTrans;/**@brief Number of translation symmetry*/ - int NSROptItrStep;/**@brief Number of iterations for stocastic reconfiguration*/ - int NSROptItrSmp;/**@brief Number of steps for sampling*/ - int NSROptFixSmp;/**@brief */ - double DSROptRedCut;/**@brief Stocastic reconfiguration parameter, input from file.*/ - double DSROptStaDel;/**@brief Stocastic reconfiguration parameter, input from file.*/ - double DSROptStepDt;/**@brief Stocastic reconfiguration parameter, input from file.*/ - int NVMCWarmUp;/**@brief */ - int NVMCInterval;/**@brief */ - int NVMCSample;/**@brief */ - int NExUpdatePath;/**@brief */ - int RndSeed;/**@brief */ - int NSplitSize;/**@brief */ - int NSPStot;/**@brief */ - int NStore;/**@brief */ - int NSRCG;/**@brief */ - int ComplexType;/**@brief */ + int NMPTrans;/**<@brief Number of translation symmetry*/ + int NSROptItrStep;/**<@brief Number of iterations for stocastic reconfiguration*/ + int NSROptItrSmp;/**<@brief Number of steps for sampling*/ + int NSROptFixSmp;/**<@brief */ + double DSROptRedCut;/**<@brief Stocastic reconfiguration parameter, input from file.*/ + double DSROptStaDel;/**<@brief Stocastic reconfiguration parameter, input from file.*/ + double DSROptStepDt;/**<@brief Stocastic reconfiguration parameter, input from file.*/ + int NVMCWarmUp;/**<@brief */ + int NVMCInterval;/**<@brief */ + int NVMCSample;/**<@brief */ + int NExUpdatePath;/**<@brief */ + int RndSeed;/**<@brief */ + int NSplitSize;/**<@brief */ + int NSPStot;/**<@brief */ + int NStore;/**<@brief */ + int NSRCG;/**<@brief */ + int ComplexType;/**<@brief */ /* Sub-lattice */ - int Lsub;/**@brief Sublattice*/ - int Wsub;/**@brief Sublattice*/ - int Hsub;/**@brief Sublattice*/ - int NCellsub;/**@brief Number of cells in a sublattice*/ - int boxsub[3][3];/**@brief Sublattice*/ - int rboxsub[3][3];/**@brief Sublattice*/ + int Lsub;/**<@brief Sublattice*/ + int Wsub;/**<@brief Sublattice*/ + int Hsub;/**<@brief Sublattice*/ + int NCellsub;/**<@brief Number of cells in a sublattice*/ + int boxsub[3][3];/**<@brief Sublattice*/ + int rboxsub[3][3];/**<@brief Sublattice*/ /* 2-body part of the trial wavefunction */ - int **Orb;/**@brief [StdIntList::nsite][StdIntList::nsite] Orbital index*/ - int **AntiOrb;/**@brief [StdIntList::nsite][StdIntList::nsite] Anti-periodic switch*/ - int NOrb;/**@brief Number of independent orbital index*/ - int NSym;/**@brief Number of translation symmetries, + int **Orb;/**<@brief [StdIntList::nsite][StdIntList::nsite] Orbital index*/ + int **AntiOrb;/**<@brief [StdIntList::nsite][StdIntList::nsite] Anti-periodic switch*/ + int NOrb;/**<@brief Number of independent orbital index*/ + int NSym;/**<@brief Number of translation symmetries, Defined from the number of cells in the sub-lattice.*/ #endif }; diff --git a/src/StdFace/TriangularLattice.c b/src/StdFace/TriangularLattice.c index 94a25e2e..3c621080 100644 --- a/src/StdFace/TriangularLattice.c +++ b/src/StdFace/TriangularLattice.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for the triangular lattice +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -24,14 +27,12 @@ along with this program. If not, see . #include /** - * - * Setup a Hamiltonian for the Triangular lattice - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Triangular(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Triangular lattice +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Triangular(struct StdIntList *StdI) { - int isite, jsite, kCell; + int isite, jsite, kCell, ntransMax, nintrMax; int iL, iW; FILE *fp; double complex Cphase; @@ -39,8 +40,8 @@ void StdFace_Triangular(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.gp", "w"); /**/ @@ -60,7 +61,9 @@ void StdFace_Triangular(struct StdIntList *StdI, char *model) /**/ StdFace_InitSite(StdI, fp, 2); StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; StdI->tau[0][2] = 0.0; - /**/ + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_J("J1'", StdI->J1pAll, StdI->J1p); StdFace_NotUsed_J("J2'", StdI->J2pAll, StdI->J2p); @@ -123,8 +126,9 @@ void StdFace_Triangular(struct StdIntList *StdI, char *model) }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; if (strcmp(StdI->model, "kondo") == 0 ) StdI->nsite *= 2; @@ -139,30 +143,28 @@ void StdFace_Triangular(struct StdIntList *StdI, char *model) StdI->locspinflag[iL] = StdI->S2; StdI->locspinflag[iL + StdI->nsite / 2] = 0; } - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + 3/*J*/ + 3/*J'*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 3/*J*/ + 3/*J'*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 6/*t*/ + 6/*t'*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (3/*V*/ + 3/*V'*/)); + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 6/*t*/ + 6/*t'*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (3/*V*/ + 3/*V'*/)); if (strcmp(StdI->model, "kondo") == 0) { - StdI->ntrans += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); + ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); }/*if (strcmp(StdI->model, "kondo") == 0)*/ } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++) { /**/ iW = StdI->Cell[kCell][0]; diff --git a/src/StdFace/Wannier90.c b/src/StdFace/Wannier90.c index 63203f09..60e465a6 100644 --- a/src/StdFace/Wannier90.c +++ b/src/StdFace/Wannier90.c @@ -15,6 +15,9 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ +/**@file +@brief Standard mode for wannier90 +*/ #include "StdFace_vals.h" #include "StdFace_ModelUtil.h" #include @@ -22,27 +25,28 @@ along with this program. If not, see . #include #include #include - -/* -* Read Geometry file for wannier90 -* -* @author Mitsuaki Kawamura (The University of Tokyo) +/** +@brief Read Geometry file for wannier90 +@author Mitsuaki Kawamura (The University of Tokyo) */ -static void geometry_W90(struct StdIntList *StdI, int *wan_num) { +static void geometry_W90( + struct StdIntList *StdI,//!<[inout] + int *wan_num//!<[out] +) { int isite, ii, ierr; FILE *fp; fprintf(stdout, " Wannier90 Geometry file = %s\n", StdI->W90_geom); fp = fopen(StdI->W90_geom, "r"); - /* - Direct lattice vector + /**@brief + Direct lattice vector StdIntList::direct */ for (ii = 0; ii < 3; ii++) ierr = fscanf(fp, "%lf%lf%lf", &StdI->direct[ii][0], &StdI->direct[ii][1], &StdI->direct[ii][2]); if(ierr != 0) printf("%d\n", ierr); - /* - Site position + /**@brief + Intrinsic site position StdIntList::tau and its number StdIntList::NsiteUC */ for (isite = 0; isite < StdI->NsiteUC; isite++) free(StdI->tau[isite]); free(StdI->tau); @@ -65,12 +69,11 @@ static void geometry_W90(struct StdIntList *StdI, int *wan_num) { StdI->tau[isite][0], StdI->tau[isite][1], StdI->tau[isite][2]); }/*static void geometry_W90(struct StdIntList *StdI) */ -/* -* Read Wannier90 hamiltonian file -* -* @author Mitsuaki Kawamura (The University of Tokyo) +/** +@brief Read Wannier90 hamiltonian file (*_hr) +@author Mitsuaki Kawamura (The University of Tokyo) */ -static void read_W90(struct StdIntList *StdI, char *model) +static void read_W90(struct StdIntList *StdI) { FILE *fp; int ierr, nWan, nWSC, iWSC, jWSC, iWan, jWan, iWan0, jWan0, ii; @@ -127,8 +130,8 @@ static void read_W90(struct StdIntList *StdI, char *model) t_tot[iWSC][iWan][jWan] = t0[wan_num[iWan]][wan_num[jWan]]; } } - /* - Inversion symmetry + /**@brief + (1) Apply inversion symmetry and delete duplication */ for (jWSC = 0; jWSC < iWSC; jWSC++) { if ( @@ -150,8 +153,8 @@ static void read_W90(struct StdIntList *StdI, char *model) t_tot[iWSC][iWan][jWan] = 0.0; } } - /* - Max t + /**@brief + (2) Search maximum transfer for appling cutoff later */ for (iWan = 0; iWan < StdI->NsiteUC; iWan++) { for (jWan = 0; jWan < StdI->NsiteUC; jWan++) { @@ -164,10 +167,13 @@ static void read_W90(struct StdIntList *StdI, char *model) fprintf(stdout, " Maximum Hopping = %f\n", tmax); fprintf(stdout, " Threshold for Hopping = %f\n", tmax * StdI->W90_cutoff); - /* - Cut-off of Hopping + /**@brief + (3) Apply cut-off of Hopping + */ + /**@brief + (3-1) Set the number of t with cut-off (StdIntList::W90_nt) + with the inputted cut-off StdIntList::W90_cutoff */ - /* Query */ StdI->W90_nt = 0; for (iWSC = 0; iWSC < nWSC; iWSC++) { for (iWan = 0; iWan < StdI->NsiteUC; iWan++) { @@ -177,12 +183,14 @@ static void read_W90(struct StdIntList *StdI, char *model) } } fprintf(stdout, " Total number of EFFECTIVE Hopping = %d\n", StdI->W90_nt); - + /**@brief + Then malloc and store to the hopping Integeral StdIntList::W90_t and + its site index StdIntList::W90_indx + */ StdI->W90_t = (double complex *)malloc(sizeof(double complex) * StdI->W90_nt); StdI->W90_indx = (int **)malloc(sizeof(int*) * StdI->W90_nt); for (ii = 0; ii < StdI->W90_nt; ii++) StdI->W90_indx[ii] = (int *)malloc(sizeof(int) * 5); - /* Store */ fprintf(stdout, " EFFECTIVE Hoppings:\n"); StdI->W90_nt = 0; for (iWSC = 0; iWSC < nWSC; iWSC++) { @@ -199,9 +207,9 @@ static void read_W90(struct StdIntList *StdI, char *model) creal(StdI->W90_t[StdI->W90_nt]), cimag(StdI->W90_t[StdI->W90_nt])); StdI->W90_nt += 1; } - } - } - } + }/*for (jWan = 0; jWan < StdI->NsiteUC; jWan++)*/ + }/*for (iWan = 0; iWan < StdI->NsiteUC; iWan++)*/ + }/*for (iWSC = 0; iWSC < nWSC; iWSC++)*/ for (iWSC = 0; iWSC < nWSC; iWSC++) { for (iWan = 0; iWan < nWan; iWan++) { @@ -217,14 +225,14 @@ static void read_W90(struct StdIntList *StdI, char *model) free(wan_num); }/*static void read_W90(struct StdIntList *StdI, char *model)*/ /** - * - * Setup a Hamiltonian for the Wannier90 *_hr.dat - * - * @author Mitsuaki Kawamura (The University of Tokyo) - */ -void StdFace_Wannier90(struct StdIntList *StdI, char *model) +@brief Setup a Hamiltonian for the Wannier90 *_hr.dat +@author Mitsuaki Kawamura (The University of Tokyo) +*/ +void StdFace_Wannier90( + struct StdIntList *StdI//!<[inout] +) { - int isite, jsite, ispin; + int isite, jsite, ispin, ntransMax, nintrMax; int iL, iW, iH, kCell, it, ii; double Jtmp[3][3] = { {0.0} }; FILE *fp; @@ -233,8 +241,8 @@ void StdFace_Wannier90(struct StdIntList *StdI, char *model) fprintf(stdout, "\n"); fprintf(stdout, "####### Parameter Summary #######\n"); fprintf(stdout, "\n"); - /* - Initialize Cell + /**@brief + (1) Compute the shape of the super-cell and sites in the super-cell */ fp = fopen("lattice.xsf", "w"); /**/ @@ -244,8 +252,10 @@ void StdFace_Wannier90(struct StdIntList *StdI, char *model) StdFace_PrintVal_d("phase1", &StdI->phase[1], 0.0); StdFace_PrintVal_d("phase2", &StdI->phase[2], 0.0); /**/ - read_W90(StdI, model); - /**/ + read_W90(StdI); + /**@brief + (2) check & store parameters of Hamiltonian + */ fprintf(stdout, "\n @ Hamiltonian \n\n"); StdFace_NotUsed_d("K", StdI->K); StdFace_PrintVal_d("h", &StdI->h, 0.0); @@ -263,8 +273,9 @@ void StdFace_Wannier90(struct StdIntList *StdI, char *model) StdFace_exit(-1); }/*if (model != "spin")*/ fprintf(stdout, "\n @ Numerical conditions\n\n"); - /* - Local Spin + /**@brief + (3) Set local spin flag (StdIntList::locspinflag) and + the number of sites (StdIntList::nsite) */ StdI->nsite = StdI->NsiteUC * StdI->NCell; StdI->locspinflag = (int *)malloc(sizeof(int) * StdI->nsite); @@ -273,25 +284,23 @@ void StdFace_Wannier90(struct StdIntList *StdI, char *model) for (isite = 0; isite < StdI->nsite; isite++) StdI->locspinflag[isite] = StdI->S2; else if(strcmp(StdI->model, "hubbard") == 0 ) for (isite = 0; isite < StdI->nsite; isite++) StdI->locspinflag[isite] = 0; - /* - The number of Transfer & Interaction + /**@brief + (4) Compute the upper limit of the number of Transfer & Interaction and malloc them. */ if (strcmp(StdI->model, "spin") == 0 ) { - StdI->ntrans = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); - StdI->nintr = StdI->NCell * (StdI->NsiteUC/*D*/ + StdI->W90_nt/*J*/) + ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/); + nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + StdI->W90_nt/*J*/) * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1); } else if (strcmp(StdI->model, "hubbard") == 0) { - StdI->ntrans = StdI->NCell * 2/*spin*/ * (2*StdI->NsiteUC/*mu+h+Gamma*/ + StdI->W90_nt * 2/*t*/); - StdI->nintr = StdI->NCell * StdI->NsiteUC/*U*/; + ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + StdI->W90_nt * 2/*t*/); + nintrMax = StdI->NCell * StdI->NsiteUC/*U*/; } /**/ - StdFace_MallocInteractions(StdI); - /* - Set Transfer & Interaction + StdFace_MallocInteractions(StdI, ntransMax, nintrMax); + /**@brief + (5) Set Transfer & Interaction */ - StdI->ntrans = 0; - StdI->nintr = 0; for (kCell = 0; kCell < StdI->NCell; kCell++){ /**/ iW = StdI->Cell[kCell][0]; diff --git a/src/mVMC/calgrn_fsz.c b/src/mVMC/calgrn_fsz.c index 0a361d9b..3eda8725 100644 --- a/src/mVMC/calgrn_fsz.c +++ b/src/mVMC/calgrn_fsz.c @@ -98,7 +98,7 @@ void CalculateGreenFunc_fsz(const double w, const double complex ip, int *eleIdx rk = CisAjsCktAltDCIdx[idx][4]; u = CisAjsCktAltDCIdx[idx][5]; rl = CisAjsCktAltDCIdx[idx][6]; - v = CisAjsCktAltDCIdx[idx][5]; + v = CisAjsCktAltDCIdx[idx][7]; if(s==t && u==v){ tmp = GreenFunc2_fsz(ri,rj,rk,rl,s,u,ip,myEleIdx,eleCfg,myEleNum,eleProjCnt,myEleSpn, diff --git a/src/mVMC/include/stcopt.h b/src/mVMC/include/stcopt.h new file mode 100644 index 00000000..7729aafa --- /dev/null +++ b/src/mVMC/include/stcopt.h @@ -0,0 +1,9 @@ +#ifndef _STCOPT_HEADER +#define _STCOPT_HEADER +#include + +int StochasticOpt(MPI_Comm comm); +void stcOptInit(double *const s, double *const g, const int nSmat, const int *const smatToParaIdx); +int stcOptMain(double *g, const int nSmat, const int *smatToParaIdx, MPI_Comm); + +#endif diff --git a/src/mVMC/include/stcopt_dposv.h b/src/mVMC/include/stcopt_dposv.h index 7ee91b8a..f8b617a9 100644 --- a/src/mVMC/include/stcopt_dposv.h +++ b/src/mVMC/include/stcopt_dposv.h @@ -3,10 +3,7 @@ #define _STCOPT_DPOSV #include -int StochasticOpt(MPI_Comm comm); -void stcOptInit(double *const s, double *const g, const int nSmat, int *const smatToParaIdx); -int stcOptMain(double *const s, double *const g, const int nSmat); - +void stcOptInit(double *const s, double *const g, const int nSmat, const int *const smatToParaIdx); #endif #endif diff --git a/src/mVMC/include/stcopt_pdposv.h b/src/mVMC/include/stcopt_pdposv.h index 67a4a12d..9259ffe4 100644 --- a/src/mVMC/include/stcopt_pdposv.h +++ b/src/mVMC/include/stcopt_pdposv.h @@ -2,8 +2,6 @@ #ifndef _STCOPT_PDPOSV #define _STCOPT_PDPOSV -int StochasticOpt(MPI_Comm comm); -int stcOptMain(double *r, const int nSmat, const int *smatToParaIdx, MPI_Comm comm); int StochasticOptDiag(MPI_Comm comm); int stcOptMainDiag(double *const r, int const nSmat, int *const smatToParaIdx, MPI_Comm comm, int const optNum); diff --git a/src/mVMC/include/vmcmain.h b/src/mVMC/include/vmcmain.h index 7b37753d..8b447045 100644 --- a/src/mVMC/include/vmcmain.h +++ b/src/mVMC/include/vmcmain.h @@ -65,6 +65,8 @@ extern int omp_get_thread_num(void); #include "../vmcclock.c" #include "../workspace.c" +#include "../stcopt.c" + #ifdef _lapack #include "../stcopt_dposv.c" #else diff --git a/src/mVMC/initfile.c b/src/mVMC/initfile.c index 7765b207..815f971f 100644 --- a/src/mVMC/initfile.c +++ b/src/mVMC/initfile.c @@ -106,7 +106,7 @@ void InitFilePhysCal(int i, int rank) { } if(NLanczosMode>0){ - sprintf(fileName, "%s_ls_%03d.dat", CDataFileHead, idx); + sprintf(fileName, "%s_ls_out_%03d.dat", CDataFileHead, idx); FileLS = fopen(fileName, "w"); sprintf(fileName, "%s_ls_qqqq_%03d.dat", CDataFileHead, idx); diff --git a/src/mVMC/matrix.c b/src/mVMC/matrix.c index 7ab4120d..beb6ec92 100644 --- a/src/mVMC/matrix.c +++ b/src/mVMC/matrix.c @@ -472,7 +472,7 @@ int calculateMAll_BF_fcmp_child( /* calculate Pf M */ M_ZSKPFA(&uplo, &mthd, &n, invM, &lda, &pfaff, iwork, work, &lwork, rwork, &info); //TBC if(info!=0) return info; - if(!isfinite(pfaff)) return qpidx+1; + if(!(isfinite(creal(pfaff)) && isfinite(cimag(pfaff)))) return qpidx+1; PfM[qpidx] = pfaff; /* DInv */ diff --git a/src/mVMC/stcopt.c b/src/mVMC/stcopt.c new file mode 100644 index 00000000..d375b2b3 --- /dev/null +++ b/src/mVMC/stcopt.c @@ -0,0 +1,191 @@ +/* +mVMC - A numerical solver package for a wide range of quantum lattice models based on many-variable Variational Monte Carlo method +Copyright (C) 2016 The University of Tokyo, All rights reserved. + +This program is developed based on the mVMC-mini program +(https://github.com/fiber-miniapp/mVMC-mini) +which follows "The BSD 3-Clause License". + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see http://www.gnu.org/licenses/. +*/ +/*------------------------------------------------------------- + * Variational Monte Carlo + * Stochastic Reconfiguration method + *------------------------------------------------------------- + * by Satoshi Morita + *-------------------------------------------------------------*/ + +// #define _DEBUG_STCOPT + +#include "stcopt.h" + +int StochasticOpt(MPI_Comm comm) { + const int nPara=NPara; + const int srOptSize=SROptSize; + const double complex *srOptOO=SROptOO; + //const double *srOptOO= SROptOO_real; + + double *r; /* the parameter change */ + int nSmat; + int smatToParaIdx[2*NPara];//TBC + + int cutNum=0,optNum=0; + double sDiag,sDiagMax,sDiagMin; + double diagCutThreshold; + + int si; /* index for matrix S */ + int pi; /* index for variational parameters */ + + double rmax; + int simax; + int info=0; + +// for real + int int_x,int_y,j,i; + + double complex *para=Para; + + int rank,size; + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&size); + + r = (double*)calloc(2*SROptSize, sizeof(double)); + + StartTimer(50); +//[s] for only real variables TBC + if(AllComplexFlag==0 && iFlgOrbitalGeneral==0){ //real & sz=0 + #pragma omp parallel for default(shared) private(i,int_x,int_y,j) + #pragma loop noalias + for(i=0;i<2*SROptSize*(2*SROptSize+2);i++){ + int_x = i%(2*SROptSize); + int_y = (i-int_x)/(2*SROptSize); + if(int_x%2==0 && int_y%2==0){ + j = int_x/2+(int_y/2)*SROptSize; + SROptOO[i] = SROptOO_real[j];// only real part TBC + }else{ + SROptOO[i] = 0.0+0.0*I; + } + } + } +//[e] + #pragma omp parallel for default(shared) private(pi) + #pragma loop noalias + for(pi=0;pi<2*nPara;pi++) { + //for(pi=0;pisDiagMax) sDiagMax=sDiag; + if(sDiag restricted parameters , pi -> full paramer 0 <-> 2*NPara + si += 1; + } +// e + } + nSmat = si; + for(si=nSmat;si<2*nPara;si++) { + smatToParaIdx[si] = -1; // parameters that will not be optimized + } + +#ifdef _DEBUG_STCOPT + printf("DEBUG in %s (%d): diagCutThreshold = %lg\n", __FILE__, __LINE__, diagCutThreshold); + printf("DEBUG in %s (%d): optNum, cutNum, nSmat, 2*nPara == %d, %d, %d, %d\n", __FILE__, __LINE__, optNum, cutNum, nSmat, 2*nPara); +#endif + + StopTimer(50); + StartTimer(51); + + + //printf("DEBUG: nSmat=%d \n",nSmat); + /* calculate r[i]: global vector [nSmat] */ + info = stcOptMain(r, nSmat, smatToParaIdx, comm); + + StopTimer(51); + StartTimer(52); + + /*** print zqp_SRinfo.dat ***/ + if(rank==0) { + if(info!=0) fprintf(stderr, "StcOpt: DPOSV info=%d\n",info); + rmax = r[0]; simax=0;; + for(si=0;si 2*NPars - - StartTimer(50); - - for(pi=0;pi<2*NPara;pi++) { - /* sDiagElm is temporarily used for diagonal elements of S */ - /* S[i][i] = OO[pi+1][pi+1] - OO[0][pi+1] * OO[0][pi+1]; */ - sDiagElm[pi] = creal(SROptOO[(pi+2)*(2*SROptSize)+(pi+2)]) - creal(SROptOO[pi+2]) * creal(SROptOO[pi+2]); - -#ifdef _DEBUG_STCOPT_DPOSV - fprintf(stderr, "DEBUG in %s (%d): sDiagElm[%d] = %lf\n", __FILE__, __LINE__, pi, sDiagElm[pi]); -#endif - } - - sDiag = sDiagElm[0]; - sDiagMax=sDiag; sDiagMin=sDiag; - for(pi=0;pi<2*NPara;pi++) { - sDiag = sDiagElm[pi]; - if(sDiag>sDiagMax) sDiagMax=sDiag; - if(sDiagsDiagMax) sDiagMax=sDiag; - if(sDiag restricted parameters , pi -> full paramer 0 <-> 2*NPara - si += 1; - } -// e - } - nSmat = si; - for(si=nSmat;si<2*nPara;si++) { - smatToParaIdx[si] = -1; // parameters that will not be optimized - } - -#ifdef _DEBUG_STCOPT_PDPOSV - printf("DEBUG in %s (%d): diagCutThreshold = %lg\n", __FILE__, __LINE__, diagCutThreshold); - printf("DEBUG in %s (%d): optNum, cutNum, nSmat, 2*nPara == %d, %d, %d, %d\n", __FILE__, __LINE__, optNum, cutNum, nSmat, 2*nPara); -#endif - - StopTimer(50); - StartTimer(51); - - - //printf("DEBUG: nSmat=%d \n",nSmat); - /* calculate r[i]: global vector [nSmat] */ - info = stcOptMain(r, nSmat, smatToParaIdx, comm); - - StopTimer(51); - StartTimer(52); - - /*** print zqp_SRinfo.dat ***/ - if(rank==0) { - if(info!=0) fprintf(stderr, "StcOpt: DPOSV info=%d\n",info); - rmax = r[0]; simax=0;; - for(si=0;si