diff --git a/.github/workflows/deploy_docs.yml b/.github/workflows/deploy_docs.yml index 737823f2..b413eb6a 100644 --- a/.github/workflows/deploy_docs.yml +++ b/.github/workflows/deploy_docs.yml @@ -11,29 +11,30 @@ on: jobs: deploy: - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 steps: - name: Inject slug/short variables - uses: rlespinasse/github-slug-action@v3.x + uses: rlespinasse/github-slug-action@v4.x - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v4 with: path: main - name: Checkout gh-pages - uses: actions/checkout@v2 + uses: actions/checkout@v4 with: ref: gh-pages path: gh-pages - name: Set up Python - uses: actions/setup-python@v1 + uses: actions/setup-python@v5 with: python-version: 3.8 - name: Prepare LaTeX env run: | + sudo apt update sudo apt install \ texlive-latex-recommended texlive-latex-extra \ texlive-lang-japanese texlive-fonts-recommended texlive-fonts-extra latexmk @@ -42,7 +43,7 @@ jobs: - name: Install python packages run: | python -m pip install --upgrade pip - pip install sphinx + pip install sphinx sphinxcontrib-spelling - name: Build run: | diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 4783d793..bb5e6bea 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -4,18 +4,24 @@ on: [push] jobs: build: - runs-on: ubuntu-20.04 + runs-on: ${{ matrix.os }} + + strategy: + matrix: + os: ['ubuntu-20.04', 'ubuntu-22.04'] + fail-fast: false steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - name: Setup Python - uses: actions/setup-python@v2 + uses: actions/setup-python@v5 with: python-version: 3.8 - name: apt run: | + sudo apt update sudo apt install \ gcc gfortran cmake wget liblapack-dev \ openmpi-bin libopenmpi-dev @@ -23,6 +29,8 @@ jobs: - name: pip install run: | python -m pip install -U pip + python -m pip install numpy scipy toml + python -m pip install typing_extensions - name: make workspace run: cmake -E make_directory ${{runner.workspace}}/build @@ -30,7 +38,7 @@ jobs: - name: cmake working-directory: ${{runner.workspace}}/build shell: bash - run: cmake -DCMAKE_VERBOSE_MAKEFILE=1 $GITHUB_WORKSPACE + run: cmake -DPYTHON_EXECUTABLE=$(which python3) -DCMAKE_INSTALL_PREFIX=${{runner.workspace}}/usr -DCMAKE_VERBOSE_MAKEFILE=1 $GITHUB_WORKSPACE - name: build working-directory: ${{runner.workspace}}/build @@ -41,3 +49,41 @@ jobs: working-directory: ${{runner.workspace}}/build shell: bash run: ctest -V + + - name: install + working-directory: ${{runner.workspace}}/build + shell: bash + run: make install + + + docs: + runs-on: ubuntu-22.04 + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: 3.8 + + - name: Prepare LaTeX env + run: | + sudo apt update + sudo apt install \ + texlive-latex-recommended texlive-latex-extra \ + texlive-lang-japanese texlive-fonts-recommended texlive-fonts-extra latexmk + kanji-config-updmap-sys ipaex + + - name: Install python packages + run: | + python -m pip install --upgrade pip + pip install sphinx sphinxcontrib-spelling + + - name: Build + run: | + cd ${GITHUB_WORKSPACE} + cmake -E make_directory build + cd build + cmake -DDocument=ON ../ + make doc diff --git a/.gitignore b/.gitignore index 8da692da..9bb07c87 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,9 @@ -### https://raw.github.com/github/gitignore/ec246076319913acee4aaeef8caf86b78e586e7a/C++.gitignore +work +usr +build + +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/C++.gitignore # Prerequisites *.d @@ -32,136 +37,8 @@ *.exe *.out *.app - - -### https://raw.github.com/github/gitignore/ec246076319913acee4aaeef8caf86b78e586e7a/Global/Vim.gitignore - -# Swap -[._]*.s[a-v][a-z] -[._]*.sw[a-p] -[._]s[a-rt-v][a-z] -[._]ss[a-gi-z] -[._]sw[a-p] - -# Session -Session.vim - -# Temporary -.netrwhist -*~ -# Auto-generated tag files -tags -# Persistent undo -[._]*.un~ - - -### https://raw.github.com/github/gitignore/ec246076319913acee4aaeef8caf86b78e586e7a/Global/Emacs.gitignore - -# -*- mode: gitignore; -*- -*~ -\#*\# -/.emacs.desktop -/.emacs.desktop.lock -*.elc -auto-save-list -tramp -.\#* - -# Org-mode -.org-id-locations -*_archive - -# flymake-mode -*_flymake.* - -# eshell files -/eshell/history -/eshell/lastdir - -# elpa packages -/elpa/ - -# reftex files -*.rel - -# AUCTeX auto folder -/auto/ - -# cask packages -.cask/ -dist/ - -# Flycheck -flycheck_*.el - -# server auth directory -/server/ - -# projectiles files -.projectile - -# directory configuration -.dir-locals.el - - -### https://raw.github.com/github/gitignore/ec246076319913acee4aaeef8caf86b78e586e7a/Global/macOS.gitignore - -# General -.DS_Store -.AppleDouble -.LSOverride - -# Icon must end with two \r -Icon - -# Thumbnails -._* - -# Files that might appear in the root of a volume -.DocumentRevisions-V100 -.fseventsd -.Spotlight-V100 -.TemporaryItems -.Trashes -.VolumeIcon.icns -.com.apple.timemachine.donotpresent - -# Directories potentially created on remote AFP share -.AppleDB -.AppleDesktop -Network Trash Folder -Temporary Items -.apdisk - - -### https://raw.github.com/github/gitignore/ec246076319913acee4aaeef8caf86b78e586e7a/Global/Windows.gitignore - -# Windows thumbnail cache files -Thumbs.db -ehthumbs.db -ehthumbs_vista.db - -# Dump file -*.stackdump - -# Folder config file -[Dd]esktop.ini - -# Recycle Bin used on file shares -$RECYCLE.BIN/ - -# Windows Installer files -*.cab -*.msi -*.msix -*.msm -*.msp - -# Windows shortcuts -*.lnk - - -### https://raw.github.com/github/gitignore/ec246076319913acee4aaeef8caf86b78e586e7a/Python.gitignore +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Python.gitignore # Byte-compiled / optimized / DLL files __pycache__/ @@ -185,6 +62,7 @@ parts/ sdist/ var/ wheels/ +share/python-wheels/ *.egg-info/ .installed.cfg *.egg @@ -210,8 +88,10 @@ htmlcov/ nosetests.xml coverage.xml *.cover +*.py,cover .hypothesis/ .pytest_cache/ +cover/ # Translations *.mo @@ -221,6 +101,7 @@ coverage.xml *.log local_settings.py db.sqlite3 +db.sqlite3-journal # Flask stuff: instance/ @@ -233,6 +114,7 @@ instance/ docs/_build/ # PyBuilder +.pybuilder/ target/ # Jupyter Notebook @@ -243,10 +125,38 @@ profile_default/ ipython_config.py # pyenv -.python-version - -# celery beat schedule file +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/#use-with-ide +.pdm.toml + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff celerybeat-schedule +celerybeat.pid # SageMath parsed files *.sage.py @@ -278,4 +188,558 @@ dmypy.json # Pyre type checker .pyre/ +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/CMake.gitignore + +CMakeLists.txt.user +CMakeCache.txt +CMakeFiles +CMakeScripts +Testing +Makefile +cmake_install.cmake +install_manifest.txt +compile_commands.json +CTestTestfile.cmake +_deps +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Global/macOS.gitignore + +# General +.DS_Store +.AppleDouble +.LSOverride + +# Icon must end with two \r +Icon + +# Thumbnails +._* + +# Files that might appear in the root of a volume +.DocumentRevisions-V100 +.fseventsd +.Spotlight-V100 +.TemporaryItems +.Trashes +.VolumeIcon.icns +.com.apple.timemachine.donotpresent + +# Directories potentially created on remote AFP share +.AppleDB +.AppleDesktop +Network Trash Folder +Temporary Items +.apdisk +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Global/Windows.gitignore + +# Windows thumbnail cache files +Thumbs.db +Thumbs.db:encryptable +ehthumbs.db +ehthumbs_vista.db + +# Dump file +*.stackdump + +# Folder config file +[Dd]esktop.ini + +# Recycle Bin used on file shares +$RECYCLE.BIN/ + +# Windows Installer files +*.cab +*.msi +*.msix +*.msm +*.msp + +# Windows shortcuts +*.lnk +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Global/Vim.gitignore + +# Swap +[._]*.s[a-v][a-z] +!*.svg # comment out if you don't need vector files +[._]*.sw[a-p] +[._]s[a-rt-v][a-z] +[._]ss[a-gi-z] +[._]sw[a-p] + +# Session +Session.vim +Sessionx.vim + +# Temporary +.netrwhist +*~ +# Auto-generated tag files +tags +# Persistent undo +[._]*.un~ +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Global/Emacs.gitignore + +# -*- mode: gitignore; -*- +*~ +\#*\# +/.emacs.desktop +/.emacs.desktop.lock +*.elc +auto-save-list +tramp +.\#* + +# Org-mode +.org-id-locations +*_archive + +# flymake-mode +*_flymake.* + +# eshell files +/eshell/history +/eshell/lastdir + +# elpa packages +/elpa/ + +# reftex files +*.rel + +# AUCTeX auto folder +/auto/ + +# cask packages +.cask/ +dist/ + +# Flycheck +flycheck_*.el + +# server auth directory +/server/ + +# projectiles files +.projectile + +# directory configuration +.dir-locals.el + +# network security +/network-security.data + +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Global/JetBrains.gitignore + +# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio, WebStorm and Rider +# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839 + +# User-specific stuff +.idea/**/workspace.xml +.idea/**/tasks.xml +.idea/**/usage.statistics.xml +.idea/**/dictionaries +.idea/**/shelf + +# AWS User-specific +.idea/**/aws.xml + +# Generated files +.idea/**/contentModel.xml + +# Sensitive or high-churn files +.idea/**/dataSources/ +.idea/**/dataSources.ids +.idea/**/dataSources.local.xml +.idea/**/sqlDataSources.xml +.idea/**/dynamic.xml +.idea/**/uiDesigner.xml +.idea/**/dbnavigator.xml + +# Gradle +.idea/**/gradle.xml +.idea/**/libraries + +# Gradle and Maven with auto-import +# When using Gradle or Maven with auto-import, you should exclude module files, +# since they will be recreated, and may cause churn. Uncomment if using +# auto-import. +# .idea/artifacts +# .idea/compiler.xml +# .idea/jarRepositories.xml +# .idea/modules.xml +# .idea/*.iml +# .idea/modules +# *.iml +# *.ipr + +# CMake +cmake-build-*/ + +# Mongo Explorer plugin +.idea/**/mongoSettings.xml + +# File-based project format +*.iws + +# IntelliJ +out/ + +# mpeltonen/sbt-idea plugin +.idea_modules/ + +# JIRA plugin +atlassian-ide-plugin.xml + +# Cursive Clojure plugin +.idea/replstate.xml + +# SonarLint plugin +.idea/sonarlint/ + +# Crashlytics plugin (for Android Studio and IntelliJ) +com_crashlytics_export_strings.xml +crashlytics.properties +crashlytics-build.properties +fabric.properties + +# Editor-based Rest Client +.idea/httpRequests + +# Android studio 3.1+ serialized cache file +.idea/caches/build_file_checksums.ser +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/TeX.gitignore + +## Core latex/pdflatex auxiliary files: +*.aux +*.lof +*.log +*.lot +*.fls +*.out +*.toc +*.fmt +*.fot +*.cb +*.cb2 +.*.lb + +## Intermediate documents: +*.dvi +*.xdv +*-converted-to.* +# these rules might exclude image files for figures etc. +# *.ps +# *.eps +# *.pdf + +## Generated if empty string is given at "Please type another file name for output:" +.pdf + +## Bibliography auxiliary files (bibtex/biblatex/biber): +*.bbl +*.bcf +*.blg +*-blx.aux +*-blx.bib +*.run.xml + +## Build tool auxiliary files: +*.fdb_latexmk +*.synctex +*.synctex(busy) +*.synctex.gz +*.synctex.gz(busy) +*.pdfsync + +## Build tool directories for auxiliary files +# latexrun +latex.out/ + +## Auxiliary and intermediate files from other packages: +# algorithms +*.alg +*.loa + +# achemso +acs-*.bib + +# amsthm +*.thm + +# beamer +*.nav +*.pre +*.snm +*.vrb + +# changes +*.soc + +# comment +*.cut + +# cprotect +*.cpt + +# elsarticle (documentclass of Elsevier journals) +*.spl + +# endnotes +*.ent + +# fixme +*.lox + +# feynmf/feynmp +*.mf +*.mp +*.t[1-9] +*.t[1-9][0-9] +*.tfm + +#(r)(e)ledmac/(r)(e)ledpar +*.end +*.?end +*.[1-9] +*.[1-9][0-9] +*.[1-9][0-9][0-9] +*.[1-9]R +*.[1-9][0-9]R +*.[1-9][0-9][0-9]R +*.eledsec[1-9] +*.eledsec[1-9]R +*.eledsec[1-9][0-9] +*.eledsec[1-9][0-9]R +*.eledsec[1-9][0-9][0-9] +*.eledsec[1-9][0-9][0-9]R + +# glossaries +*.acn +*.acr +*.glg +*.glo +*.gls +*.glsdefs +*.lzo +*.lzs +*.slg +*.slo +*.sls + +# uncomment this for glossaries-extra (will ignore makeindex's style files!) +# *.ist + +# gnuplot +*.gnuplot +*.table + +# gnuplottex +*-gnuplottex-* + +# gregoriotex +*.gaux +*.glog +*.gtex + +# htlatex +*.4ct +*.4tc +*.idv +*.lg +*.trc +*.xref + +# hyperref +*.brf + +# knitr +*-concordance.tex +# TODO Uncomment the next line if you use knitr and want to ignore its generated tikz files +# *.tikz +*-tikzDictionary + +# listings +*.lol + +# luatexja-ruby +*.ltjruby + +# makeidx +*.idx +*.ilg +*.ind + +# minitoc +*.maf +*.mlf +*.mlt +*.mtc[0-9]* +*.slf[0-9]* +*.slt[0-9]* +*.stc[0-9]* + +# minted +_minted* +*.pyg + +# morewrites +*.mw + +# newpax +*.newpax + +# nomencl +*.nlg +*.nlo +*.nls + +# pax +*.pax + +# pdfpcnotes +*.pdfpc + +# sagetex +*.sagetex.sage +*.sagetex.py +*.sagetex.scmd + +# scrwfile +*.wrt + +# svg +svg-inkscape/ + +# sympy +*.sout +*.sympy +sympy-plots-for-*.tex/ + +# pdfcomment +*.upa +*.upb + +# pythontex +*.pytxcode +pythontex-files-*/ + +# tcolorbox +*.listing + +# thmtools +*.loe + +# TikZ & PGF +*.dpth +*.md5 +*.auxlock + +# titletoc +*.ptc + +# todonotes +*.tdo + +# vhistory +*.hst +*.ver + +# easy-todo +*.lod + +# xcolor +*.xcp + +# xmpincl +*.xmpi + +# xindy +*.xdy + +# xypic precompiled matrices and outlines +*.xyc +*.xyd + +# endfloat +*.ttt +*.fff + +# Latexian +TSWLatexianTemp* + +## Editors: +# WinEdt +*.bak +*.sav + +# Texpad +.texpadtmp + +# LyX +*.lyx~ + +# Kile +*.backup + +# gummi +.*.swp + +# KBibTeX +*~[0-9]* + +# TeXnicCenter +*.tps + +# auto folder when using emacs and auctex +./auto/* +*.el + +# expex forward references with \gathertags +*-tags.tex + +# standalone packages +*.sta + +# Makeindex log files +*.lpz + +# xwatermark package +*.xwm + +# REVTeX puts footnotes in the bibliography by default, unless the nofootinbib +# option is specified. Footnotes are the stored in a file with suffix Notes.bib. +# Uncomment the next line to have this generated file ignored. +#*Notes.bib +### Generated by gibo (https://github.com/simonwhitaker/gibo) +### https://raw.github.com/github/gitignore/4488915eec0b3a45b5c63ead28f286819c0917de/Global/VisualStudioCode.gitignore + +.vscode/* +!.vscode/settings.json +!.vscode/tasks.json +!.vscode/launch.json +!.vscode/extensions.json +!.vscode/*.code-snippets + +# Local History for Visual Studio Code +.history/ +# Built Visual Studio Code Extensions +*.vsix diff --git a/CMakeLists.txt b/CMakeLists.txt index bfd54552..017b7859 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ -cmake_minimum_required(VERSION 2.8 FATAL_ERROR) +cmake_minimum_required(VERSION 2.8.12 FATAL_ERROR) project(DSQSS NONE) -set(DSQSS_VERSION 2.0.4) +set(DSQSS_VERSION 2.1-dev) message(STATUS "CMake version: " ${CMAKE_VERSION}) @@ -66,78 +66,48 @@ include_directories(${Boost_INCLUDE_DIRS}) include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src/third-party/plog) -if(BUILD_NEW_GENERATORS OR Testing OR Document) - if(${CMAKE_VERSION} VERSION_LESS 3.12) - find_package(PythonInterp 3.4) - if(NOT PYTHONINTERP_FOUND) - unset(PYTHON_EXECUTABLE CACHE) - find_package(PythonInterp 2.7 REQUIRED) - endif(NOT PYTHONINTERP_FOUND) - set(python_version_mm "${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}") - else(${CMAKE_VERSION} VERSION_LESS 3.12) - find_package(Python3 3.4 COMPONENTS Interpreter) - if(Python3_FOUND) +set(python_version_required 3.8) + +if(NOT PYTHON_EXECUTABLE) + if(BUILD_NEW_GENERATORS OR Testing OR Document) + if(${CMAKE_VERSION} VERSION_LESS 3.12) + find_package(PythonInterp ${python_version_required} REQUIRED) + set(python_version_mm "${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}") + else(${CMAKE_VERSION} VERSION_LESS 3.12) + find_package(Python3 ${python_version_required} COMPONENTS Interpreter REQUIRED) set(PYTHON_EXECUTABLE ${Python3_EXECUTABLE}) set(python_version_mm "${Python3_VERSION_MAJOR}.${Python3_VERSION_MINOR}") - else(Python3_FOUND) - find_package(Python2 2.7 COMPONENTS Interpreter REQUIRED) - set(PYTHON_EXECUTABLE ${Python2_EXECUTABLE}) - set(python_version_mm "${Python2_VERSION_MAJOR}.${Python2_VERSION_MINOR}") - endif(Python3_FOUND) - endif(${CMAKE_VERSION} VERSION_LESS 3.12) - - include(${PROJECT_SOURCE_DIR}/config/FindPythonModule.cmake) - set(pythonpath_prefix "${CMAKE_CURRENT_BINARY_DIR}/pythonmodule") - set(pythonpath_build "${pythonpath_prefix}/lib/python${python_version_mm}/site-packages") - set(ENV{PYTHONPATH} "${pythonpath_build}:$ENV{PYTHONPATH}") - set(ENV{PATH} "${pythonpath_prefix}/bin:$ENV{PATH}") - - set(module_names "") - if(BUILD_NEW_GENERATORS OR Testing) - set(module_names numpy scipy toml ${module_names}) + endif(${CMAKE_VERSION} VERSION_LESS 3.12) endif() - if(Document) - set(module_names ${module_names} sphinx) - set(sphinxcontribs spelling) +else() + # check python version + execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-V" OUTPUT_VARIABLE result ERROR_VARIABLE result) + string(REGEX MATCH "^Python[ ]+[0-9]+\\.[0-9]+\\.[0-9]+" result "${result}") + if(NOT result) + message(FATAL_ERROR "${PYTHON_EXECUTABLE} seems not python executable") endif() - - foreach(module_name ${module_names}) - find_python_module(${module_name}) - string(TOUPPER ${module_name} uppername) - if(NOT ${uppername}_FOUND) - set(missing_modules ${missing_modules} ${module_name}) - endif(NOT ${uppername}_FOUND) - endforeach() - foreach(module_name ${sphinxcontribs}) - find_python_module("sphinxcontrib.${module_name}") - string(TOUPPER ${module_name} uppername) - if(NOT SPHINXCONTRIB.${uppername}_FOUND) - set(missing_modules ${missing_modules} "sphinxcontrib-${module_name}") - endif() - endforeach() - list(LENGTH missing_modules num_missing_modules) - if(num_missing_modules) - find_python_module(pip) - if(NOT PIP_FOUND) - file(DOWNLOAD https://bootstrap.pypa.io/get-pip.py "${CMAKE_CURRENT_BINARY_DIR}/get-pip.py") - execute_process(COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/get-pip.py --prefix=${pythonpath_prefix}) - execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pip install -U pip) - endif(NOT PIP_FOUND) - foreach(module_name ${missing_modules}) - message(STATUS "env PIP_USER=0 PYTHONPATH=${pythonpath_build}:\$PYTHONPATH ${PYTHON_EXECUTABLE} -m pip install --prefix=${pythonpath_prefix} ${module_name}") - execute_process( - COMMAND env PIP_USER=0 PYTHONPATH=${pythonpath_build}:\$PYTHONPATH ${PYTHON_EXECUTABLE} -m pip install --prefix=${pythonpath_prefix} ${module_name} - ) - endforeach(module_name ${missing_modules}) - endif(num_missing_modules) - + message(STATUS "Given Python interpreter: ${PYTHON_EXECUTABLE}") + string(REGEX MATCH "[0-9]+\\.[0-9]+\\.[0-9]+" result "${result}") + message(STATUS "Python version: ${result}") + if("${result}" VERSION_LESS ${python_version_required}) + message(FATAL_ERROR "Python ${result} is not supported (${python_version_required} or later is required)") + endif() + string(REPLACE "." ";" result "${result}") + set(result "${result}") + list(GET result 0 Python3_VERSION_MAJOR) + list(GET result 1 Python3_VERSION_MINOR) + list(GET result 2 Python3_VERSION_MICRO) + set(python_version_mm "${Python3_VERSION_MAJOR}.${Python3_VERSION_MINOR}") endif() add_subdirectory(src/dla) if(MPI_FOUND) add_subdirectory(src/pmwa) endif(MPI_FOUND) -add_subdirectory(tool) + +if(BUILD_NEW_GENERATORS) + add_subdirectory(tool) +endif() if (Testing) enable_testing() diff --git a/README.md b/README.md index 057976eb..54ead129 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,7 @@ -# DSQSS +# DSQSS (Discrete Space Quantum Systems Solver) -[![Build Status](https://travis-ci.org/issp-center-dev/dsqss.svg?branch=master)](https://travis-ci.org/issp-center-dev/dsqss) -[![doc latest_stable en](https://img.shields.io/badge/doc--en-v2.0.1-blue.svg)](https://issp-center-dev.github.io/dsqss/manual/v2.0.1/en/index.html) -[![doc latest_stable jp](https://img.shields.io/badge/doc--jp-v2.0.1-blue.svg)](https://issp-center-dev.github.io/dsqss/manual/v2.0.1/jp/index.html) +[![doc latest_stable en](https://img.shields.io/badge/doc--en-v2.0.4-blue.svg)](https://issp-center-dev.github.io/dsqss/manual/v2.0.4/en/index.html) +[![doc latest_stable jp](https://img.shields.io/badge/doc--jp-v2.0.4-blue.svg)](https://issp-center-dev.github.io/dsqss/manual/v2.0.4/jp/index.html) [![doc master](https://img.shields.io/badge/doc-dev-blue.svg)](https://issp-center-dev.github.io/dsqss/manual/master/en/index.html) [![doc master jp](https://img.shields.io/badge/doc--jp-dev-blue.svg)](https://issp-center-dev.github.io/dsqss/manual/master/jp/index.html) @@ -10,6 +9,7 @@ DSQSS (Discrete Space Quantum Systems Solver) is a software package for calculat DSQSS implements the path-integral Monte Carlo method with the directed loop algorithm. ## Structure of directories + - src - source codes of the main programs, DSQSS/DLA and DSQSS/PMWA - tool @@ -28,8 +28,11 @@ DSQSS implements the path-integral Monte Carlo method with the directed loop alg ### prerequisite - C++ Compiler -- CMake 2.8=< -- Python 2.7 or 3.4=< +- CMake >=2.8.12 +- Python >=3.8 + - numpy + - scipy + - toml ### Simple build @@ -42,10 +45,8 @@ make ``` #### Notes -DSQSS requires some python packages, and if not found DSQSS installs them automatically into local while CMake process. -If `pip` is too old and does not know `--prefix` option of `install`, please update `pip` . -To change compiler, add `-DCMAKE_CXX_COMPILER` like +To change compiler, add `-DCMAKE_CXX_COMPILER` like ``` bash cmake ../ -DCMAKE_CXX_COMPILER=`which icpc` @@ -96,9 +97,9 @@ cmake -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} ../ make install ``` -After this process, executable files such as `dla` are installed to `${INSTALL_DIR}/bin`, +After this process, executable files such as `dla` are installed to `${INSTALL_DIR}/bin`, a configuration file `dsqssvars-${DSQSS_VERSION}.sh` to `${INSTALL_DIR}/share/dsqss/`, -sample files to `${INSTALL_DIR}/share/dsqss/dsqss-${DSQSS_VERSION}/sample`, +sample files to `${INSTALL_DIR}/share/dsqss/dsqss-${DSQSS_VERSION}/sample`, and documents (if built) to `${INSTALL_DIR}/share/dsqss/dsqss-${DSQSS_VERSION}/doc`. ### Example @@ -113,14 +114,25 @@ cat sample.log | grep ene For details, [see manual page](https://issp-center-dev.github.io/dsqss/manual/master/en/dla/tutorial/spindimer.html). +## Paper + +We would appreciate if you cite the following article in your research with DSQSS, +[Y. Motoyama, K. Yoshimi, A. Masaki-Kato, T. Kato, and N. Kawashima, Comput. Phys. Commun. 264, 107944 (2021)](https://www.sciencedirect.com/science/article/pii/S0010465521000692). + +The preprint is [arXiv:2007.11329](https://arxiv.org/abs/2007.11329). + ## License + ### License of DSQSS + DSQSS is distributed under the GNU GPL v3. ### License of the bundled libs + - [Boost C++ library](https://www.boost.org/) is redistributed under the Boost software license. - [Plog](https://github.com/SergiusTheBest/plog) is redistributed under the Mozilla Public License 2.0 - [`FindPythonModule.cmake`](https://github.com/openturns/openturns/tree/master/cmake.FindPythonModule.cmake) is redistributed under the OSI approved BSD license. ## Acknowledgement + DSQSS v1.2 and v2.0 are developed under the support of "Project for advancement of software usability in materials science" in fiscal year 2018 by The Institute for Solid State Physics, The University of Tokyo. diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 57af432a..1b0a9092 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -9,6 +9,12 @@ endif() find_package(LATEX) +# For Apple silicon, pyenchant (used by sphinxcontrib.spelling) may fail to find the enchant library. +# https://github.com/pyenchant/pyenchant/issues/265 +find_library(ENCHANT_LIBRARIES + NAMES enchant-2 + ) + add_subdirectory(en) add_subdirectory(jp) diff --git a/doc/en/CMakeLists.txt b/doc/en/CMakeLists.txt index 12ba5c36..27b896ea 100644 --- a/doc/en/CMakeLists.txt +++ b/doc/en/CMakeLists.txt @@ -6,25 +6,40 @@ set(DSQSS_DOC_DIR share/dsqss/${DSQSS_VERSION}/doc/en) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/conf.py ${CMAKE_CURRENT_BINARY_DIR}) +set(cmd_list + env PYTHONPATH="${pythonpath_build}:$ENV{PYTHONPATH}" + ${SPHINX_EXECUTABLE} + -b html + -d ${SPHINX_CACHE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR} + ${SPHINX_HTML_DIR} +) + +if(ENCHANT_LIBRARIES) + list(INSERT cmd_list 0 env PYENCHANT_LIBRARY_PATH=${ENCHANT_LIBRARIES}) +endif() + add_custom_target(doc-en-html ALL COMMAND - env PYTHONPATH="${pythonpath_build}:$ENV{PYTHONPATH}" - ${SPHINX_EXECUTABLE} - -b html - -d ${SPHINX_CACHE_DIR} - ${CMAKE_CURRENT_SOURCE_DIR} - ${SPHINX_HTML_DIR} + ${cmd_list} ) if(LATEX_FOUND) - add_custom_target(doc-en-pdf ALL - COMMAND + set(cmd_list env PYTHONPATH="${pythonpath_build}:$ENV{PYTHONPATH}" + ${ENC_ENV} ${SPHINX_EXECUTABLE} -b latex -d ${SPHINX_CACHE_DIR} ${CMAKE_CURRENT_SOURCE_DIR} ${SPHINX_PDF_DIR} + ) + if(ENCHANT_LIBRARIES) + list(INSERT cmd_list 0 env PYENCHANT_LIBRARY_PATH=${ENCHANT_LIBRARIES}) + endif() + add_custom_target(doc-en-pdf ALL + COMMAND + ${cmd_list} COMMAND cd pdf && make ) diff --git a/doc/en/dla/tutorial/spinchain.rst b/doc/en/dla/tutorial/spinchain.rst index 4b34dc80..0135003c 100644 --- a/doc/en/dla/tutorial/spinchain.rst +++ b/doc/en/dla/tutorial/spinchain.rst @@ -55,7 +55,7 @@ Before executing this script, ``source`` a configuring file ``dsqssvars-VERSION. The result of :math:`S=1/2,1` will be written to ``xmzu_1.dat`` and ``xmzu_2.dat``, respectively (:numref:`fig_spinchain`). The :math:`S=1/2` chain is gapless and so the magnetic susceptibility remains finite at absolute zero (note that in the simulation the finite size effect opens energy gap). -On the other hand the magnetic susceptibility of the :math:`S=1` chain drops to zero at finite temperature due to the spin gap. +On the other hand, the magnetic susceptibility of the :math:`S=1` chain drops to zero at finite temperature due to the spin gap. .. figure:: ../../../image/dla/tutorial/spinchain.* :name: fig_spinchain diff --git a/doc/en/dla/tutorial/standard_mode b/doc/en/dla/tutorial/standard_mode new file mode 100644 index 00000000..f5cdd9d9 --- /dev/null +++ b/doc/en/dla/tutorial/standard_mode @@ -0,0 +1,226 @@ +Definition of models and lattices of antiferromagnetic Heisenberg dimers by standard mode of DSQSS/DLA +=============================================================== + +In this tutorial, we will consider a system in which an external magnetic field is applied to a :math:`S=1/2` antiferromagnetic Heisenberg dimer :math:`\mathcal{H}= \sum_{\langle i,j \rangle}-J_z S^z_i S^z_j -\frac{J_{xy}}{2}(S^+_i S^-_j +S^-_i S^+_j)-h\sum_i S^z_i ,~J_z<0, J_xy<0, h>0` +by calculating the approximate form of the lattice, +you learn how to use the standard mode of DSQSS/DLA. + +The calculations in the standard mode of DSQSS/DLA are as follows + +1. prepare an input file for simple mode, +2. create an input file for standard mode, +3. run the calculation. + + + +Preparing the input file for simple mode +******************** + +To run DSQSS/DLA, you will need the simple mode file std.toml in TOML format with the following information + +#. Hamiltonian Information, +#. Information on lattices, +#. Information on parameters such as inverse temperature and number of Monte Carlo steps. + +Thus, the first step is to prepare this input file. +The file of antiferromagnetic Heisenberg dimer is available in the following location +(sample/dla/01_spindimer/std.toml) +:: + + [hamiltonian] + model = "spin" + M = 1 # S=1/2 + Jz = -1.0 # coupling constant, negative for AF + Jxy = -1.0 # coupling constant, negative for AF + h = 0.0 # magnetic field + [lattice] + lattice = "hypercubic" # hypercubic, periodic + dim = 1 # dimension + L = 2 # number of sites along each direction + bc = false # open boundary + [parameter] + beta = 100 # inverse temperature + nset = 5 # set of Monte Carlo sweeps + npre = 10 # MCSteps to estimate hyperparameter + ntherm = 10 # MCSweeps for thermalization + nmcs = 100 # MCSweeps for measurement + seed = 31415 # seed of RNG + +To this file, we add the effect of applying an external magnetic field in the z direction, +and set the lattice dimension and the number of sites to 2 and 4, respectively +:: + + [hamiltonian] + model = "spin" + M = 1 # S=1/2 + Jz = -1.0 # coupling constant, negative for AF + Jxy = -1.0 # coupling constant, negative for AF + h = 0.5 # magnetic field + [lattice] + lattice = "hypercubic" # hypercubic, periodic + dim = 2 # dimension + L = 4 # number of sites along each direction + bc = false # open boundary + [parameter] + beta = 100 # inverse temperature + nset = 5 # set of Monte Carlo sweeps + npre = 10 # MCSteps to estimate hyperparameter + ntherm = 10 # MCSweeps for thermalization + nmcs = 100 # MCSweeps for measurement + seed = 31415 # seed of RNG + + +Writing an input file for standard mode +********** + +This file is given to dla_hamgen +:: + + $ dla_hamgen std.toml + +As a result, a text file in TOML format, hamiltonian.toml, +is generated to describe the information of the site Hamiltonian and the many-body interaction. + + +Next, you give std.toml to dla_latgen +:: + + $ dla_latgen -t lattice.toml std.toml + + +As a result, a text file in TOML format, lattice.toml, +is generated to describe the information of the space using unit cells and basic translation vectors. + + +Running the calculation +**************** + +With lattice.toml as input file, you can create a gnuplot file lattice.plt:: +:: + + $ dla_latgen -g lattice.plt lattice.toml + +Giving lattce.plt to gnuplot will produce a two-dimensional schematic of the lattice:: + + +.. figure:: ../../../image/dla/tutorial/2Dlattice.png + :name: fig_2Dlattice + :alt: Two-dimensional lattice + + Schematic of the lattice of an antiferromagnetic Heisenberg spin chain in an external magnetic field in the z direction. + + + +Writing lattice.dat and kpoints.dat +******************** + +The lattice.toml file defines the lattice information using unit cells and fundamental translation vectors, +while the lattice.dat file defines the lattice information using the number of sites and interactions between sites. +The lattice.dat file also can be generated using the simple mode file std.toml in TOML format. +You give std.toml to dla_latgen, which contains the effect of applying an external magnetic field in the z direction as used earlier +:: + + $ dla_latgen -o lattice.dat std.toml + +As a result, you will get the lattice.dat file as follows +:: + + name + hypercubic + + lattice + 2 # dim + 4 4 # size + 0 0 # 0:open boundary, 1:periodic boundary + 0 1.0 0.0 # latvec_0 + 1 0.0 1.0 # latvec_1 + + directions + 2 # ndirections + # id, coords... + 0 1.0 0.0 + 1 0.0 1.0 + + sites + 16 # nsites + # id, type, coord... + 0 0 0.0 0.0 + 1 0 1.0 0.0 + 2 0 2.0 0.0 + 3 0 3.0 0.0 + 4 0 0.0 1.0 + 5 0 1.0 1.0 + 6 0 2.0 1.0 + 7 0 3.0 1.0 + 8 0 0.0 2.0 + 9 0 1.0 2.0 + 10 0 2.0 2.0 + 11 0 3.0 2.0 + 12 0 0.0 3.0 + 13 0 1.0 3.0 + 14 0 2.0 3.0 + 15 0 3.0 3.0 + + interactions + 24 # nints + # id, type, nbody, sites..., edge_flag, direction + 0 0 2 0 1 0 0 + 1 0 2 0 4 0 1 + 2 1 2 1 2 0 0 + 3 2 2 1 5 0 1 + 4 3 2 2 3 0 0 + 5 2 2 2 6 0 1 + 6 0 2 3 7 0 1 + 7 2 2 4 5 0 0 + 8 1 2 4 8 0 1 + 9 4 2 5 6 0 0 + 10 4 2 5 9 0 1 + 11 5 2 6 7 0 0 + 12 4 2 6 10 0 1 + 13 1 2 7 11 0 1 + 14 2 2 8 9 0 0 + 15 3 2 8 12 0 1 + 16 4 2 9 10 0 0 + 17 5 2 9 13 0 1 + 18 5 2 10 11 0 0 + 19 5 2 10 14 0 1 + 20 3 2 11 15 0 1 + 21 0 2 12 13 0 0 + 22 1 2 13 14 0 0 + 23 3 2 14 15 0 0 + + +Next, you can create kpoints.dat, +which is a text file that specifies the coefficients of the reciprocal grid vector expansion of the wavenumber vector. +You can specify the size by using -s. Please give the std.toml file to dla_wvgen -s :: + + $ dla_wvgen -s "8 8" std.toml + +As a result, you will get the following kpoints.dat file +:: + + dim + 2 + + kpoints + 4 + 0 0 0 + 1 0 4 + 2 4 0 + 3 4 4 + + + + + + + + + + + + + + + + diff --git a/doc/en/dla/users-manual/generator.rst b/doc/en/dla/users-manual/generator.rst index c7144f49..f8611e06 100644 --- a/doc/en/dla/users-manual/generator.rst +++ b/doc/en/dla/users-manual/generator.rst @@ -21,7 +21,7 @@ Simple mode tool ``dla_pre`` The meanings of the parameters are following. ``paramfile`` - The name of the parameter file to be generated (default: ``qmc.inp``.) + The name of the parameter file to be generated (default: ``param.in``.) ``inputfile`` The name of the input file. diff --git a/doc/en/dla/users-manual/input_expert.rst b/doc/en/dla/users-manual/input_expert.rst index ae84b0e2..11ede5fa 100644 --- a/doc/en/dla/users-manual/input_expert.rst +++ b/doc/en/dla/users-manual/input_expert.rst @@ -12,7 +12,7 @@ The list of input files :header-rows: 0 :widths: 1,4 - qmc.inp, "Parameter list for the simulation, e.g., the number of Monte Carlo sets." + param.in, "Parameter list for the simulation, e.g., the number of Monte Carlo sets." lattice.xml, "Definition of the lattice." algorithm.xml, "Definition of the algorithm (e.g., scattering rate of a worm)." wavevector.xml, "Indication of wave vectors for structure factors. (optional)" diff --git a/doc/en/dla/users-manual/output.rst b/doc/en/dla/users-manual/output.rst index 85190a98..8d6ea4c1 100644 --- a/doc/en/dla/users-manual/output.rst +++ b/doc/en/dla/users-manual/output.rst @@ -118,6 +118,13 @@ Notations :math:`\beta` The inverse temperature. +:math:`E_0` + The imaginary time average of the expectation value of the unperturbed Hamiltonian + :math:`\displaystyle \frac{1}{\beta}\int d\tau \langle \phi(\tau)|\mathcal{H}_0|\phi(\tau)\rangle`. + +:math:`N_v` + The number of vertices, i.e., the order of the perturbation. + :math:`h` The conjugate field to the operator :math:`M^z` . The longitudinal magnetic field for spin systems and the chemical potential for the Bose-Hubbard models. @@ -129,13 +136,17 @@ Main results ***************** Main results are written in a file with the name specified by ``outfile`` keyword in the input parameter file. +NOTICE: In general, Monte Carlo simulations have systematic errors of :math:`O(1/N)` with respect to the number of samples :math:`N` (``nmcs``) for the expectation values including nonlinear functions of the sample average like the specific heat and the susceptibility. +For example, in the region where the specific heat becomes very small, e.g., in the low-temperature region below the energy gap, the calculated value may be negative. +For precise analysis, we need to take into account not only the statistical errors but also the systematic errors. + ``sign`` The sign of the weights. :math:`\sum_i W_i / \sum_i |W_i` ``anv`` - The mean number of the vertices. + The number of the vertices per site. :math:`\displaystyle \frac{\langle N_v \rangle}{N_s}` ``ene`` @@ -145,7 +156,7 @@ Main results are written in a file with the name specified by ``outfile`` keywor ``spe`` The specific heat - :math:`\displaystyle C_V \equiv \frac{\partial \epsilon}{\partial T}` + :math:`\displaystyle C_V \equiv \frac{\partial \epsilon}{\partial T} = \frac{1}{N_s T^2} \left[\left\langle\left(E_0 - TN_v\right)^2\right\rangle - \left\langle\left(E_0 - TN_v\right)\right\rangle^2 - T^2\left\langle N_v \right\rangle\right]` ``som`` The ratio of the specific heat and the temperature. @@ -210,6 +221,10 @@ Main results are written in a file with the name specified by ``outfile`` keywor :math:`\displaystyle \chi^{zz}(\vec{k}, \omega=0) = \beta N_s \left[\left\langle (\tilde{m}_K^z)^2 \right\rangle - \left\langle \tilde{m}_K^z \right\rangle^2 \right]` +``ds1`` + The derivative of the "magnetization" (``amzu``) with respect to the temperature. + + :math:`\displaystyle T\frac{\partial \left\langle \tilde{m}^z \right\rangle}{\partial T} = -\beta\frac{\partial \left\langle \tilde{m}^z \right\rangle}{\partial\beta}` ``wi2`` The winding number. @@ -270,6 +285,8 @@ Displacement :math:`\vec{r}_{ij}` and imaginary time :math:`\tau` are specified where ```` is an index of the displacement specified by ``kind`` (the first element of each ``R`` tag) in the relative coordinate XML file, and ```` is an index of the discretized imaginary time. +NOTE: The current version, this works only for :math:`S=1/2` model. + Momentum space temperature Green's function output ************************************************** The momentum space temperature Green's function is written into a file with the name specified by ``ckoutfile`` keyword in the input file. @@ -281,3 +298,5 @@ The momentum space temperature Green's function is defined as the following: Wave vector :math:`\vec{r}_{ij}` and imaginary time :math:`\tau` are specified by the name ``Ct`` as the same way of structure factor, where ```` is an index of the displacement specified by ``kind`` (the last element of each ``RK`` tag) in the wavevector XML file, and ```` is an index of the discretized imaginary time. + +NOTE: The current version, this works only for :math:`S=1/2` model. diff --git a/doc/en/dsqss/about_dsqss.rst b/doc/en/dsqss/about_dsqss.rst index 5f87a5be..454608ff 100644 --- a/doc/en/dsqss/about_dsqss.rst +++ b/doc/en/dsqss/about_dsqss.rst @@ -6,7 +6,7 @@ About DSQSS Overview **************** -DSQSS is a program package for solving quantum many-body problems defined on lattices. It is based on the quantum Monte Carlo method in Feynman's path integral representation. It covers a broad range of problems written by flexible input files that define arbitrary unit cells in arbitrary dimensions, and arbitrary matrix elements of the interactions among arbitrary number of degrees of freedom. +DSQSS is a program package for solving quantum many-body problems defined on lattices. It is based on the quantum Monte Carlo method in Feynman's path integral representation. It covers a broad range of problems written by flexible input files that define arbitrary unit cells in arbitrary dimensions and arbitrary matrix elements of the interactions among arbitrary number of degrees of freedom. For example, you can perform finite temperature calculation of XXZ spin model by specifying parameters such as dimension, size of lattice, anisotropic coupling constants, length of spin, strength of magnetic field, and temperature. You can calculate Bose-Hubbard model as well as quantum spin model. PMWA (Parallel Multi Worm Algorithm) suits for large-scale non-trivial parallel calculation by domain parallelization. @@ -50,11 +50,11 @@ Collaborators License **************** - GNU General Public License (GPL) -The users are kindly requested to acknowledge the usage of this software in their publication, if any, based on the software, and let the developers know its reference information. +- The users are kindly requested to cite the following paper in their publication, if any, based on the software. -.. topic:: Acknowledgment Sample +.. topic:: Reference - Numerical results in the present paper were obtained by the quantum Monte Carlo program DSQSS(https://github.com/qmc/dsqss/wiki). This package is distributed under GNU General Public License version 3 (GPL v3) or later. + Y. Motoyama, K. Yoshimi, A. Masaki-Kato, T. Kato, and N. Kawashima, Comput. Phys. Commun. 264, 107944 (2021). Acknowledgment **************** diff --git a/doc/en/dsqss/install.rst b/doc/en/dsqss/install.rst index d743cb82..d4441825 100644 --- a/doc/en/dsqss/install.rst +++ b/doc/en/dsqss/install.rst @@ -8,7 +8,7 @@ Requirements ******************** - (Optional) MPI (essential for PMWA) -- python 2.7 or 3.4+ +- python 3.8+ - numpy - scipy @@ -84,6 +84,15 @@ It is noted that the default install directory is set as ``/usr/local/bin`` . For details, see https://github.com/issp-center-dev/HPhi/wiki/FAQ . + +.. note:: + + CMake automatically search for an interpreter of Python. + If you want to use a specific one, + you should tell it to CMake by using ``-DPYTHON_EXECUTABLE`` option:: + + $ cmake ../ -DPYTHON_EXECUTABLE=`which python3` + Each binary files for dsqss will be made in ``src`` and ``tool`` directories. To check whether the binary files are correctly made or not, please type the following command: diff --git a/doc/en/pmwa/tutorial/intro.rst b/doc/en/pmwa/tutorial/intro.rst index 49b8a89f..8454ba07 100644 --- a/doc/en/pmwa/tutorial/intro.rst +++ b/doc/en/pmwa/tutorial/intro.rst @@ -3,4 +3,4 @@ About DSQSS/PMWA --------------------------- -``DSQSS / PMWA`` is a numerical package of the world-line quantum Monte Carlo method with non-local updates algorithm for large-scale parallel computing. PMWA can treat :math:`S=1/2` quantum many-body system such as XXZ spin model, Heisenberg model, transverse Ising model model and hard-core boson system with large lattice sizes at finite temperatures. Basically, any physical quantities can be calculated within the statistical error. By default, the energy under a longitudinal worm source field, longitudinal magnetization, transverse magnetization, winding number, and imaginary time :math:`S^z S^z` spin correlation function are calculated. +``DSQSS / PMWA`` is a numerical package of the world-line quantum Monte Carlo method with non-local updates algorithm for large-scale parallel computing. PMWA can treat :math:`S=1/2` quantum many-body system such as XXZ spin model, Heisenberg model, transverse Ising model and hard-core boson system with large lattice sizes at finite temperatures. Basically, any physical quantities can be calculated within the statistical error. By default, the energy under a longitudinal worm source field, longitudinal magnetization, transverse magnetization, winding number, and imaginary time :math:`S^z S^z` spin correlation function are calculated. diff --git a/doc/en/pmwa/tutorial/lattgene.rst b/doc/en/pmwa/tutorial/lattgene.rst index f618f9c6..f205d3b7 100644 --- a/doc/en/pmwa/tutorial/lattgene.rst +++ b/doc/en/pmwa/tutorial/lattgene.rst @@ -25,11 +25,11 @@ The parameters to be specified are as follows. Examples of use are described below. -1. One-dimensional 8-sites chain, inverse temperature 10.0, definition of lattice file when division number is 1: +1. One-dimensional 8-sites chain, inverse temperature 10.0, the definition of lattice file when division number is 1: ``$ lattgene_P 1 8 10.0 1 1 0`` -2. Two dimensional 4*4 sites lattice, reverse temperature 10.0, definition of lattice file when division number is 1: +2. Two-dimensional 4*4 sites lattice, reverse temperature 10.0, the definition of lattice file when division number is 1: ``$ lattgene_P 2 4 4 10.0 1 1 0`` diff --git a/doc/image/dla/tutorial/2Dlattice.png b/doc/image/dla/tutorial/2Dlattice.png new file mode 100644 index 00000000..c79e68f5 Binary files /dev/null and b/doc/image/dla/tutorial/2Dlattice.png differ diff --git a/doc/jp/CMakeLists.txt b/doc/jp/CMakeLists.txt index 4798d326..a4571947 100644 --- a/doc/jp/CMakeLists.txt +++ b/doc/jp/CMakeLists.txt @@ -6,25 +6,42 @@ set(DSQSS_DOC_DIR_JP share/dsqss/${DSQSS_VERSION}/doc/jp) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/conf.py ${CMAKE_CURRENT_BINARY_DIR}) +set(cmd_list + env PYTHONPATH="${pythonpath_build}:$ENV{PYTHONPATH}" + ${SPHINX_EXECUTABLE} + -b html + -d ${SPHINX_CACHE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR} + ${SPHINX_HTML_DIR} +) + +if(ENCHANT_LIBRARIES) + list(INSERT cmd_list 0 env PYENCHANT_LIBRARY_PATH=${ENCHANT_LIBRARIES}) +endif() + add_custom_target(doc-jp-html ALL COMMAND - env PYTHONPATH="${pythonpath_build}:$ENV{PYTHONPATH}" - ${SPHINX_EXECUTABLE} - -b html - -d ${SPHINX_CACHE_DIR} - ${CMAKE_CURRENT_SOURCE_DIR} - ${SPHINX_HTML_DIR} + ${cmd_list} ) if(LATEX_FOUND) - add_custom_target(doc-jp-pdf ALL - COMMAND + + set(cmd_list env PYTHONPATH="${pythonpath_build}:$ENV{PYTHONPATH}" + ${ENC_ENV} ${SPHINX_EXECUTABLE} -b latex -d ${SPHINX_CACHE_DIR} ${CMAKE_CURRENT_SOURCE_DIR} ${SPHINX_PDF_DIR} + ) + if(ENCHANT_LIBRARIES) + list(INSERT cmd_list 0 env PYENCHANT_LIBRARY_PATH=${ENCHANT_LIBRARIES}) + endif() + + add_custom_target(doc-jp-pdf ALL + COMMAND + ${cmd_list} COMMAND cd pdf && make ) diff --git a/doc/jp/dla/tutorial/spindimer.rst b/doc/jp/dla/tutorial/spindimer.rst index b3427188..0271a3e2 100644 --- a/doc/jp/dla/tutorial/spindimer.rst +++ b/doc/jp/dla/tutorial/spindimer.rst @@ -70,7 +70,7 @@ DSQSS/DLA を実行するには, $ mpiexec -np 4 dla param.in -乱数並列計算では並列数 (この例では4) だけ独立に計算が行われ, その分モンテカルロサンプル数を増えるために計算精度が向上します. [#fn_ompi_macos]_ +乱数並列計算では並列数 (この例では4) だけ独立に計算が行われ, その分モンテカルロサンプル数が増えるために計算精度が向上します. [#fn_ompi_macos]_ 計算結果の解釈 diff --git a/doc/jp/dla/tutorial/standard_mode.rst b/doc/jp/dla/tutorial/standard_mode.rst new file mode 100644 index 00000000..cc216bbf --- /dev/null +++ b/doc/jp/dla/tutorial/standard_mode.rst @@ -0,0 +1,220 @@ +DSQSS/DLA のスタンダードモードによる反強磁性ハイゼンベルグダイマーの模型や格子の定義 +=============================================================== + +このチュートリアルでは, :math:`S=1/2` 反強磁性ハイゼンベルグダイマーに外部磁場が印加されている系 :math:`\mathcal{H}= \sum_{\langle i,j \rangle}-J_z S^z_i S^z_j -\frac{J_{xy}}{2}(S^+_i S^-_j +S^-_i S^+_j)-h\sum_i S^z_i ,~J_z<0, J_xy<0, h>0` の格子の概形を計算をすることで, +DSQSS/DLA のスタンダードモード使い方を学びます. + +DSQSS/DLA のスタンダードモードによる計算は, + +1. シンプルモードの入力ファイルの準備 +2. スタンダードモードの入力ファイルの作成 +3. 計算の実行 + +の3段階に分かれます. + + +シンプルモードの入力ファイルの準備 +******************** + +DSQSS/DLA を実行するには, + +#. ハミルトニアンの情報 +#. 格子の情報 +#. 逆温度やモンテカルロステップ数などのパラメータの情報 + +を指定したTOML形式のシンプルモードファイル std.tomlが必要です. +そのため, まずはこの入力ファイルを準備します. +反強磁性ハイゼンベルグダイマーのファイルは以下の場所に用意されています +(sample/dla/01_spindimer/std.toml) +:: + + [hamiltonian] + model = "spin" + M = 1 # S=1/2 + Jz = -1.0 # coupling constant, negative for AF + Jxy = -1.0 # coupling constant, negative for AF + h = 0.0 # magnetic field + [lattice] + lattice = "hypercubic" # hypercubic, periodic + dim = 1 # dimension + L = 2 # number of sites along each direction + bc = false # open boundary + [parameter] + beta = 100 # inverse temperature + nset = 5 # set of Monte Carlo sweeps + npre = 10 # MCSteps to estimate hyperparameter + ntherm = 10 # MCSweeps for thermalization + nmcs = 100 # MCSweeps for measurement + seed = 31415 # seed of RNG + +このファイルにz方向の外部磁場の印加の効果を付け加え, 格子の次元とサイトの数をそれぞれ2と4にします +:: + + [hamiltonian] + model = "spin" + M = 1 # S=1/2 + Jz = -1.0 # coupling constant, negative for AF + Jxy = -1.0 # coupling constant, negative for AF + h = 0.5 # magnetic field + [lattice] + lattice = "hypercubic" # hypercubic, periodic + dim = 2 # dimension + L = 4 # number of sites along each direction + bc = false # open boundary + [parameter] + beta = 100 # inverse temperature + nset = 5 # set of Monte Carlo sweeps + npre = 10 # MCSteps to estimate hyperparameter + ntherm = 10 # MCSweeps for thermalization + nmcs = 100 # MCSweeps for measurement + seed = 31415 # seed of RNG + + +スタンダードモードの入力ファイルの作成 +********** + +このファイルを dla_hamgen に与えます. +:: + + $ dla_hamgen std.toml + +この結果, サイトハミルトニアンや多体相互作用の情報を記述するTOML形式のテキストファイル hamiltonian.toml が生成されます. + + +次に, std.toml を dla_latgen に与えます. +:: + + $ dla_latgen -t lattice.toml std.toml + + +この結果, ユニットセルと基本並進ベクトルを用いて空間の情報を記述するTOML形式のテキストファイルlattice.toml が生成されます. + + +計算の実行 +**************** + +lattice.toml を入力ファイルとして, gnuplotファイル lattice.pltを作成できます:: +:: + + $ dla_latgen -g lattice.plt lattice.toml + +lattce.plt を gnuplot に与えると2次元の格子の概形が得られます:: + + +.. figure:: ../../../image/dla/tutorial/2Dlattice.png + :name: fig_2Dlattice + :alt: 2次元格子 + + z方向の外部磁場中の反強磁性ハイゼンベルグスピン鎖の格子の概形. + + + +lattice.datとkpoints.datの作成 +******************** + +lattice.tomlファイルでは格子の情報をユニットセルと基本並進ベクトルを用いて定義しました。一方, lattice.datファイルは格子の情報をサイトの数やサイト間の相互作用を用いて定義します。lattice.datもTOML形式のシンプルモードファイル std.tomlを用いて作成できます。 +上で用いたz方向の外部磁場の印加の効果が入っているstd.tomlを dla_latgen に与えます. +:: + + $ dla_latgen -o lattice.dat std.toml + +この結果, 以下のようなlattice.datファイルが得られます +:: + + name + hypercubic + + lattice + 2 # dim + 4 4 # size + 0 0 # 0:open boundary, 1:periodic boundary + 0 1.0 0.0 # latvec_0 + 1 0.0 1.0 # latvec_1 + + directions + 2 # ndirections + # id, coords... + 0 1.0 0.0 + 1 0.0 1.0 + + sites + 16 # nsites + # id, type, coord... + 0 0 0.0 0.0 + 1 0 1.0 0.0 + 2 0 2.0 0.0 + 3 0 3.0 0.0 + 4 0 0.0 1.0 + 5 0 1.0 1.0 + 6 0 2.0 1.0 + 7 0 3.0 1.0 + 8 0 0.0 2.0 + 9 0 1.0 2.0 + 10 0 2.0 2.0 + 11 0 3.0 2.0 + 12 0 0.0 3.0 + 13 0 1.0 3.0 + 14 0 2.0 3.0 + 15 0 3.0 3.0 + + interactions + 24 # nints + # id, type, nbody, sites..., edge_flag, direction + 0 0 2 0 1 0 0 + 1 0 2 0 4 0 1 + 2 1 2 1 2 0 0 + 3 2 2 1 5 0 1 + 4 3 2 2 3 0 0 + 5 2 2 2 6 0 1 + 6 0 2 3 7 0 1 + 7 2 2 4 5 0 0 + 8 1 2 4 8 0 1 + 9 4 2 5 6 0 0 + 10 4 2 5 9 0 1 + 11 5 2 6 7 0 0 + 12 4 2 6 10 0 1 + 13 1 2 7 11 0 1 + 14 2 2 8 9 0 0 + 15 3 2 8 12 0 1 + 16 4 2 9 10 0 0 + 17 5 2 9 13 0 1 + 18 5 2 10 11 0 0 + 19 5 2 10 14 0 1 + 20 3 2 11 15 0 1 + 21 0 2 12 13 0 0 + 22 1 2 13 14 0 0 + 23 3 2 14 15 0 0 + + +次に, kpoints.datを作成します。kpoint.datは波数ベクトルを逆格子ベクトルで展開した際の係数を指定するテキストファイルです. -sを使うことでサイズを指定できます. dla_wvgen -s に std.tomlファイルを与えます:: + + $ dla_wvgen -s "8 8" std.toml + +この結果, 以下の kpoints.datファイルが得られます +:: + + dim + 2 + + kpoints + 4 + 0 0 0 + 1 0 4 + 2 4 0 + 3 4 4 + + + + + + + + + + + + + + + + diff --git a/doc/jp/dla/users-manual/generator.rst b/doc/jp/dla/users-manual/generator.rst index aae6a31d..debc5c09 100644 --- a/doc/jp/dla/users-manual/generator.rst +++ b/doc/jp/dla/users-manual/generator.rst @@ -21,7 +21,7 @@ DLA は入力ファイルとして格子定義ファイル,アルゴリズム定 パラメータは以下の通り. ``paramfile`` - 出力されるパラメータファイル名.デフォルトは ``qmc.inp`` です. + 出力されるパラメータファイル名.デフォルトは ``param.in`` です. ``inputfile`` 入力ファイル名. ファイル形式の詳細は :ref:`simple_mode_file` を参照してください. diff --git a/doc/jp/dla/users-manual/input_expert.rst b/doc/jp/dla/users-manual/input_expert.rst index f9cf1238..c710235a 100644 --- a/doc/jp/dla/users-manual/input_expert.rst +++ b/doc/jp/dla/users-manual/input_expert.rst @@ -16,7 +16,7 @@ DSQSS/DLA のエキスパートモード入力ファイルは, DSQSS/DLA の実 :header-rows: 0 :widths: 1,4 - qmc.inp, "モンテカルロの繰り返し回数など,計算制御のためのパラメータファイル." + param.in, "モンテカルロの繰り返し回数など,計算制御のためのパラメータファイル." lattice.xml, "格子の定義ファイル." algorithm.xml, "アルゴリズム（ワームの散乱確率など）の定義ファイル." wavevetor.xml, "波数定義ファイル." @@ -422,7 +422,7 @@ Algorithm/Vertex/InitialConfiguration/Channel 波数ベクトルXMLファイルは, スタッガード秩序変数 .. math:: - M^{z}(\vec{k}) \equiv \frac{1}{N} \sum_i e^{-i\vec{k}\vec\vec{r}_i} \left\langle M^{z}_i \right\rangle + M^{z}(\vec{k}) \equiv \frac{1}{N} \sum_i e^{-i\vec{k}\cdot\vec{r}_i} \left\langle M^{z}_i \right\rangle や動的構造因子 diff --git a/doc/jp/dla/users-manual/output.rst b/doc/jp/dla/users-manual/output.rst index 077d08f7..11e807de 100644 --- a/doc/jp/dla/users-manual/output.rst +++ b/doc/jp/dla/users-manual/output.rst @@ -116,6 +116,13 @@ DLA は計算結果を行区切りのプレーンテキストファイルで出 :math:`\beta` 逆温度. +:math:`E_0` + 非摂動ハミルトニアンの期待値の虚時間平均 + :math:`\displaystyle \frac{1}{\beta}\int d\tau \langle \phi(\tau)|\mathcal{H}_0|\phi(\tau)\rangle`. + +:math:`N_v` + バーテックスの数, すなわち摂動の次数. + :math:`h` :math:`M^z` に共役な外場. スピン系では縦磁場, ボース粒子系では化学ポテンシャル. @@ -127,23 +134,27 @@ DLA は計算結果を行区切りのプレーンテキストファイルで出 ***************** メイン出力ファイルは, 入力パラメータファイルの ``outfile`` キーワードで指定した名前で出力されます. +NOTICE: 一般にモンテカルロ法では, 比熱や感受率など, サンプル平均の非線形関数を含むような量の期待値には, モンテカルロサンプル数 :math:`N` (``nmcs``)に対して :math:`O(1/N)` の系統誤差も含まれます. +たとえば, エネルギーギャップ以下の極低温領域など, 比熱の値が非常に小さくなるような場合には計算結果が負になることがあります. +精密な解析には統計誤差だけではなく, 系統誤差にも注意する必要があります. + ``sign`` 重みの符号. :math:`\frac{\sum_i W_i }{ \sum_i |W_i| }`, ここで :math:`i` はモンテカルロサンプルの番号. ``anv`` - 平均バーテックス数. + サイトあたりの平均バーテックス数. :math:`\displaystyle \frac{\langle N_v \rangle}{N_s}` ``ene`` エネルギー密度. - :math:`\displaystyle \epsilon \equiv \frac{1}{N_s}\left(E_0 - T\langle N_v\rangle\right)` + :math:`\displaystyle \epsilon \equiv \frac{1}{N_s}\left(\langle E_0 \rangle - T\langle N_v\rangle\right)` ``spe`` 比熱. - :math:`\displaystyle C_V \equiv \frac{\partial \epsilon}{\partial T}` + :math:`\displaystyle C_V \equiv \frac{\partial \epsilon}{\partial T} = \frac{1}{N_s T^2} \left[\left\langle\left(E_0 - TN_v\right)^2\right\rangle - \left\langle\left(E_0 - TN_v\right)\right\rangle^2 - T^2\left\langle N_v \right\rangle\right]` ``som`` 比熱と温度の比. @@ -206,6 +217,11 @@ DLA は計算結果を行区切りのプレーンテキストファイルで出 :math:`\displaystyle \chi^{zz}(\vec{k}, \omega=0) = \beta N_s \left[\left\langle (\tilde{m}_K^z)^2 \right\rangle - \left\langle \tilde{m}_K^z \right\rangle^2 \right]` +``ds1`` + 「磁化 ``amzu`` 」の温度微分. + + :math:`\displaystyle T\frac{\partial \left\langle \tilde{m}^z \right\rangle}{\partial T} = -\beta\frac{\partial \left\langle \tilde{m}^z \right\rangle}{\partial\beta}` + ``wi2`` ワインディングナンバー. @@ -268,6 +284,9 @@ DLA は計算結果を行区切りのプレーンテキストファイルで出 ここで ```` は変位XML ファイルの ``kind`` (``R`` タグの第一要素) で指定される変位のインデックスで, ```` は離散化した虚時間のインデックス. +NOTICE: +現在のバージョンでは :math:`S=1/2` のモデル以外はうまく計算できません. + 波数表示温度グリーン関数出力ファイル **************************************** 波数表示温度グリーン関数出力ファイルは, 入力パラメータファイルの ``ckoutfile`` キーワードで指定した名前で出力されます. @@ -281,3 +300,6 @@ DLA は計算結果を行区切りのプレーンテキストファイルで出 ``Ct`` という形で物理量名によって区別されます. ここで ```` は波数ベクトルXMLファイルの ``kindex`` (``RK`` タグの最終要素) で指定される波数のインデックスで, ```` は離散化した虚時間のインデックス. + +NOTICE: +現在のバージョンでは :math:`S=1/2` のモデル以外はうまく計算できません. diff --git a/doc/jp/dsqss/about_dsqss.rst b/doc/jp/dsqss/about_dsqss.rst index c85d8e40..4c0e5cf0 100644 --- a/doc/jp/dsqss/about_dsqss.rst +++ b/doc/jp/dsqss/about_dsqss.rst @@ -94,12 +94,11 @@ DSQSSは, 以下のサブシステム群から構成されています. ライセンス **************** - GNU General Public License (GPL)に基づきます. -- 利用のための必須条件ではありませんが, 利用実態を把握したいので, 科学計算などに使用した場合, 関連公表論文の書誌情報などをアプリケーション管理者までお知らせ下さることを希望します. また, 論文などによる成果公開に際して謝辞に記載していただければ幸いです. +- 利用のための必須条件ではありませんが, DSQSS を利用した論文を執筆した際には, 以下の文献を引いていただけると幸いです. -.. topic:: Acknowledgment Sample +.. topic:: Reference - Numerical results in the present paper were obtained by the quantum Monte Carlo program DSQSS(https://github.com/issp-center-dev/dsqss). - This package is distributed under GNU General Public License version 3 (GPL v3) or later. + Y. Motoyama, K. Yoshimi, A. Masaki-Kato, T. Kato, and N. Kawashima, Comput. Phys. Commun. 264, 107944 (2021). 謝辞 diff --git a/doc/jp/dsqss/install.rst b/doc/jp/dsqss/install.rst index 686a64d3..3bd995b4 100644 --- a/doc/jp/dsqss/install.rst +++ b/doc/jp/dsqss/install.rst @@ -9,7 +9,7 @@ DSQSSの使用には以下のプログラム・ライブラリが必要です. - (Optional) MPI (PMWAを使用する場合には必須) -- python 2.7 or 3.4+ +- python 3.8+ - numpy - scipy @@ -66,7 +66,7 @@ DSQSSのダウンロード後にzipファイルを解凍すると, ファイル :: $ mkdir dsqss.build && cd dsqss.build - $ cmake ../ -DCMAKE_INSTALL_PREFIX=/path/to/install/to + $ cmake ../ -DCMAKE_INSTALL_PREFIX=/path/to/install/to $ make ``/path/to/install/to`` をインストールしたい先のパスに設定してください（例： ``$HOME/opt/dsqss`` ）. @@ -87,6 +87,12 @@ DSQSSのダウンロード後にzipファイルを解凍すると, ファイル 詳細については https://github.com/issp-center-dev/HPhi/wiki/FAQ をご覧ください. +.. note:: + + CMake は自動的に Python のインタープリタを探します。 + 明示的に指定したいときは ``-DPYTHON_EXECUTABLE`` オプションを用いてください:: + + $ cmake ../ -DPYTHON_EXECUTABLE=`which python3` これにより各実行ファイルが ``dsqss.build/src`` ディレクトリ以下に, 入力ファイルの生成ツールが ``dsqss.build/tool`` ディレクトリ以下に作成されます. diff --git a/src/common/version.h b/src/common/version.h index 470c4bde..192ffc87 100644 --- a/src/common/version.h +++ b/src/common/version.h @@ -17,6 +17,6 @@ #ifndef SRC_COMMON_VERSION_H_ #define SRC_COMMON_VERSION_H_ -#define DSQSS_VERSION "v2.0.4" +#define DSQSS_VERSION "2.1-dev" #endif // SRC_COMMON_VERSION_H_ diff --git a/src/dla/CMakeLists.txt b/src/dla/CMakeLists.txt index 3cfec297..b19bb003 100644 --- a/src/dla/CMakeLists.txt +++ b/src/dla/CMakeLists.txt @@ -1,5 +1,4 @@ # include guard -cmake_minimum_required(VERSION 2.8) if(${CMAKE_PROJECT_NAME} STREQUAL "Project") message(FATAL_ERROR "cmake should be executed not for 'src' subdirectory, but for the top directory of DLA.") endif(${CMAKE_PROJECT_NAME} STREQUAL "Project") @@ -11,13 +10,6 @@ set(dla_SOURCES random.cc dla.cc) add_executable(dla ${dla_SOURCES}) install(TARGETS dla RUNTIME DESTINATION bin) -# add_executable(dla_B ${dla_SOURCES}) -# add_executable(dla_H ${dla_SOURCES}) -# install(TARGETS dla_B RUNTIME DESTINATION bin) -# install(TARGETS dla_H RUNTIME DESTINATION bin) -# target_compile_definitions(dla_B PRIVATE MEASURE_BOSON=1) -# target_compile_definitions(dla_H PRIVATE MEASURE_SPIN=1) - if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel") target_link_libraries(dla rt) endif(CMAKE_CXX_COMPILER_ID STREQUAL "Intel") @@ -25,22 +17,12 @@ endif(CMAKE_CXX_COMPILER_ID STREQUAL "Intel") if(MPI_FOUND) add_definitions(-DMULTI) target_include_directories(dla PUBLIC ${MPI_CXX_INCLUDE_DIRS}) - # target_include_directories(dla_H PUBLIC ${MPI_CXX_INCLUDE_DIRS}) - # target_include_directories(dla_B PUBLIC ${MPI_CXX_INCLUDE_DIRS}) target_link_libraries(dla ${MPI_CXX_LIBRARIES}) - # target_link_libraries(dla_H ${MPI_CXX_LIBRARIES}) - # target_link_libraries(dla_B ${MPI_CXX_LIBRARIES}) endif(MPI_FOUND) -# set(dla_deps dla_H dla_B) set(dla_deps dla) if(BUILD_OLD_GENERATORS) add_subdirectory(generators) set(dla_deps ${dla_deps} dla_alg_old lattgene_C lattgene_T hamgen_H hamgen_B cfgene sfgene) endif(BUILD_OLD_GENERATORS) - -# add_custom_target(dla -# DEPENDS -# ${dla_deps} -# ) diff --git a/src/dla/SPECIFIC/Boson/measure_specific.cc b/src/dla/SPECIFIC/Boson/measure_specific.cc deleted file mode 100644 index 525e9d9f..00000000 --- a/src/dla/SPECIFIC/Boson/measure_specific.cc +++ /dev/null @@ -1,191 +0,0 @@ -#include "../../util.hpp" - -void Measurement::measure(double sgn) { - using namespace Specific; - int NV = LAT.countVertices(); - - double MZUA = 0.0; // uniform, tau=0 - double MZUB = 0.0; // uniform, integrated - double MZSA = 0.0; // staggered, tau=0 - double MZSB = 0.0; // staggered, integrated - - const double T = 1.0 / LAT.BETA; - const double invV = 1.0 / LAT.NSITE; - - std::vector phase(LAT.NSMTYPE); - for (int m = 0; m < LAT.NSMTYPE; ++m) { - phase[m] = std::cos(2 * M_PI * m / LAT.NSMTYPE); - } - - for (int s = 0; s < LAT.NSITE; s++) { - Site& SITE = LAT.S(s); - int mt = SITE.getMTYPE(); - double ph = phase[mt]; - Segment& S0 = SITE.first(); - double mz0 = dvals[S0.X()]; - Site::iterator p(SITE); - double mza0 = 0.0; - - while (!(++p).atOrigin()) { - Segment& S = *p; - double mz = dvals[S.X()]; - mza0 += mz * S.length(); - } - - MZUA += mz0; - MZUB += mza0; - MZSA += ph * mz0; - MZSB += ph * mza0; - } - MZUA *= invV; - MZSA *= invV; - MZUB *= invV; - MZSB *= invV; - MZUB *= T; - MZSB *= T; - - double EBSAMP = -(double)NV; - - for (int b = 0; b < LAT.NINT; b++) { - Interaction& I = LAT.I(b); - InteractionProperty& IP = I.property(); - // VertexProperty& VP = IP.getVertexProperty(); - int NBODY = IP.NBODY; - std::vector tau(NBODY); - std::vector x(NBODY); - std::vector p(NBODY); - - for (int i = 0; i < NBODY; i++) { - Site& S = I.site(i); - p[i].init(S); - ++p[i]; - tau[i] = p[i]->topTime(); - x[i] = p[i]->X(); - } - - double t = 0.0; - - while (t < LAT.BETA) { - int it = util::min_index(tau); - EBSAMP -= (tau[it] - t) * IP.VertexDensity(x); - - if (p[it]->top().isTerminal()) break; - t = tau[it]; - ++p[it]; - tau[it] = p[it]->topTime(); - x[it] = p[it]->X(); - } - } - - ACC[SGN].accumulate(sgn); - - ACC[NV1].accumulate(sgn * NV); - - ACC[MZUA1].accumulate(sgn * MZUA); - ACC[MZUA2].accumulate(sgn * MZUA * MZUA); - ACC[MZUB1].accumulate(sgn * MZUB); - ACC[MZUB2].accumulate(sgn * MZUB * MZUB); - ACC[MZSA1].accumulate(sgn * MZSA); - ACC[MZSA2].accumulate(sgn * MZSA * MZSA); - ACC[MZSB1].accumulate(sgn * MZSB); - ACC[MZSB2].accumulate(sgn * MZSB * MZSB); - - ACC[EB1].accumulate(sgn * EBSAMP); - ACC[EB2].accumulate(sgn * EBSAMP * EBSAMP); - ACC[NH1].accumulate(sgn * MZUA * EBSAMP); - - // WINDING NUMBER - const int Dim = LAT.BD; - if (Dim > 0) { - std::vector wind(Dim); - - for (int xi = 0; xi < LAT.NEDGE; xi++) { - int Asite = LAT.EDGE(xi).A; - int Bsite = LAT.EDGE(xi).B; - int dim = LAT.EDGE(xi).bd; - - Site& SITE = LAT.S(Asite); - Segment& S0 = SITE.first(); - Site::iterator p(SITE); - int xlast = S0.X(); - - // count winding - - while (!(++p).atOrigin()) { - Segment& S = *p; - if (xlast != S.X()) { - if ((*S.bottom().S(0).getONSITE()).id() - 1 == Bsite) { - wind[dim] += static_cast(S.bottom().S(1).X() - S.bottom().S(0).X()); - } else if ((*S.bottom().S(2).getONSITE()).id() - 1 == Bsite) { - wind[dim] += static_cast(S.bottom().S(3).X() - S.bottom().S(2).X()); - } - xlast = S.X(); - } - } - } // end of for(xi) - - double wind2 = 0.0; - for (int di = 0; di < LAT.D; di++) { - double W = 0; - for (int db = 0; db < Dim; db++) { - W += wind[db] * LAT.vec[di][db] * LAT.L[di]; - } - wind2 += W * W; - } - ACC[Wxy2].accumulate(sgn * wind2); - } // end of if(Dim) -} - -void Measurement::setsummary() { - using namespace Specific; - double WDIAG = ALG.getBlock( - "WDIAG", (double)1.0); // ALG.X["General"]["WDIAG" ].getDouble(); // 0.25 - - std::vector X(NACC); - - for (int i = 0; i < NACC; i++) { - ACC[i].average(); - X[i] = ACC[i].mean(); - } - - double B = LAT.BETA; - double T = 1.0 / B; - double V = LAT.NSITE; - double invV = 1.0 / V; - double D = LAT.D; - - double invsign = 1.0 / X[SGN]; - - for (int i = 0; i < NACC; ++i) { - X[SGN] *= invsign; - } - - Q[SIGN] = X[SGN]; - - Q[ANV] = X[NV1] * invV; - Q[ENE] = (EBASE + X[EB1] * T) * invV; - - Q[SPE] = (X[EB2] - X[EB1] * X[EB1] - X[NV1]) * invV; - - Q[LEN] = X[LE1]; - Q[XMX] = WDIAG * X[LE1] / B; - - Q[AMZU] = X[MZUA1]; - Q[BMZU] = X[MZUB1]; - Q[SMZU] = (X[MZUA2] - X[MZUA1] * X[MZUA1]) * V; - Q[XMZU] = (X[MZUB2] - X[MZUB1] * X[MZUB1]) * B * V; - - Q[AMZS] = X[MZSA1]; - Q[BMZS] = X[MZSB1]; - Q[SMZS] = (X[MZSA2] - X[MZSA1] * X[MZSA1]) * V; - Q[XMZS] = (X[MZSB2] - X[MZSB1] * X[MZSB1]) * B * V; - - Q[DS1] = B * (X[NH1] - X[MZUA1] * X[EB1]) / V; - Q[W2] = X[Wxy2]; - Q[RHOS] = X[Wxy2] * 0.5 / D / V / B; - Q[COMP] = Q[XMZU] / (X[MZUB1] * X[MZUB1]); - - for (int i = 0; i < NPHY; i++) { - PHY[i].accumulate(Q[i]); - } -} diff --git a/src/dla/SPECIFIC/Boson/measure_specific.h b/src/dla/SPECIFIC/Boson/measure_specific.h deleted file mode 100644 index d95ca46d..00000000 --- a/src/dla/SPECIFIC/Boson/measure_specific.h +++ /dev/null @@ -1,82 +0,0 @@ -#ifndef SRC_DLA_SPECIFIC_BOSON_MEASURE_SPECIFIC_IMPL_H_ -#define SRC_DLA_SPECIFIC_BOSON_MEASURE_SPECIFIC_IMPL_H_ - -#include - -namespace Specific { - -std::vector diagonal_operators(int NXMAX) { - std::vector ret(NXMAX); - for (int i = 0; i < NXMAX; ++i) { - ret[i] = i; - } - return ret; -} - -// accumulator specifier -enum { - SGN, // Sign of weight - NV1, // Number of Vertices - EB1, // Number of kinks - EB2, // Number of kinks^2 - LE1, // LEngth of worm pair - - // observables about Magnetization Z per site - // A: tau=0 B: tau-integrated - // U: Uniform S: Staggered - // 1: linear 2: squared - MZUA1, - MZUA2, - MZUB1, - MZUB2, - MZSA1, - MZSA2, - MZSB1, - MZSB2, - - NH1, - Wxy2, - NACC, // number of quantities measured at each MC step -}; -static const char ANAME[][NACC] = { - "sgn", "nv1", "eb1", "eb2", "le1", "mzua1", "mzua2", "mzub1", - "mzub2", "mzsa1", "mzsa2", "mzsb1", "mzsb2", "nh1", "wxy2"}; - -// observable specifier -enum { - SIGN, // Sign of weight - ANV, // (site) Average of the Number of Vertices - ENE, // ENErgy per site - SPE, // SPEcific heat - LEN, // LENgth of worm pair - XMX, - - // observables about Magnetization Z per site - // prefix: - // A: tau=0 B: tau-integrated - // S: Variance X: Susceptibility - // postfix: - // U: Uniform, S: Staggered - AMZU, - BMZU, - SMZU, - XMZU, - AMZS, - BMZS, - SMZS, - XMZS, - - DS1, - W2, // Winding number^2 - RHOS, // Superfluid density - COMP, // COMPressibility - - NPHY, // number of quantities computed at each set -}; -static const char PNAME[][NPHY] = { - "sign", "anv", "ene", "spe", "len", "xmx", "amzu", "bmzu", "smzu", - "xmzu", "amzs", "bmzs", "smzs", "xmzs", "ds1", "wi2", "rhos", "comp"}; - -} // namespace Specific - -#endif // SRC_DLA_SPECIFIC_BOSON_MEASURE_SPECIFIC_IMPL_H_ diff --git a/src/dla/SPECIFIC/Heisenberg/measure_specific.cc b/src/dla/SPECIFIC/Heisenberg/measure_specific.cc deleted file mode 100644 index 6e09a323..00000000 --- a/src/dla/SPECIFIC/Heisenberg/measure_specific.cc +++ /dev/null @@ -1,143 +0,0 @@ -#include "../../util.hpp" - -void Measurement::measure(double sgn) { - using namespace Specific; - int NV = LAT.countVertices(); - - double MZUA = 0.0; // uniform, tau=0 - double MZUB = 0.0; // uniform, integrated - double MZSA = 0.0; // staggered, tau=0 - double MZSB = 0.0; // staggered, integrated - - const double T = 1.0 / LAT.BETA; - const double invV = 1.0 / LAT.NSITE; - - std::vector phase(LAT.NSMTYPE); - for (int m = 0; m < LAT.NSMTYPE; ++m) { - phase[m] = std::cos(2.0 * M_PI * m / LAT.NSMTYPE); - } - - for (int s = 0; s < LAT.NSITE; s++) { - Site& SITE = LAT.S(s); - int mt = SITE.getMTYPE(); - double ph = phase[mt]; - Segment& S0 = SITE.first(); - double mz0 = dvals[S0.X()]; - Site::iterator p(SITE); - double mza0 = 0.0; - - while (!(++p).atOrigin()) { - Segment& S = *p; - double mz = dvals[S.X()]; - mza0 += mz * S.length(); - } - - MZUA += mz0; - MZUB += mza0; - MZSA += ph * mz0; - MZSB += ph * mza0; - } - MZUA *= invV; - MZSA *= invV; - MZUB *= invV; - MZSB *= invV; - MZUB *= T; - MZSB *= T; - - double EBSAMP = -1.0*NV; - - for (int b = 0; b < LAT.NINT; b++) { - Interaction& I = LAT.I(b); - InteractionProperty& IP = I.property(); - // VertexProperty& VP = IP.getVertexProperty(); - int NBODY = IP.NBODY; - std::vector tau(NBODY); - std::vector x(NBODY); - std::vector x2(2 * NBODY); - std::vector p(NBODY); - - for (int i = 0; i < NBODY; i++) { - Site& S = I.site(i); - p[i].init(S); - ++p[i]; - tau[i] = p[i]->topTime(); - x[i] = p[i]->X(); - } - - double t = 0.0; - int it; - - while (t < LAT.BETA) { - it = util::min_index(tau); - EBSAMP -= (tau[it] - t) * IP.VertexDensity(x); - - if (p[it]->top().isTerminal()) break; - t = tau[it]; - ++p[it]; - tau[it] = p[it]->topTime(); - x[it] = p[it]->X(); - } - } - - ACC[SGN].accumulate(sgn); - - ACC[NV1].accumulate(sgn * NV); - - ACC[MZUA1].accumulate(sgn * MZUA); - ACC[MZUA2].accumulate(sgn * MZUA * MZUA); - ACC[MZUB1].accumulate(sgn * MZUB); - ACC[MZUB2].accumulate(sgn * MZUB * MZUB); - ACC[MZSA1].accumulate(sgn * MZSA); - ACC[MZSA2].accumulate(sgn * MZSA * MZSA); - ACC[MZSB1].accumulate(sgn * MZSB); - ACC[MZSB2].accumulate(sgn * MZSB * MZSB); - - ACC[EB1].accumulate(sgn * EBSAMP); - ACC[EB2].accumulate(sgn * EBSAMP * EBSAMP); -} - -void Measurement::setsummary() { - using namespace Specific; - const double WDIAG = ALG.getBlock( - "WDIAG", (double)1.0); // ALG.X["General"]["WDIAG" ].getDouble(); // 0.25 - - std::vector X(NACC); - - for (int i = 0; i < NACC; i++) { - ACC[i].average(); - X[i] = ACC[i].mean(); - } - - const double B = LAT.BETA; - const double T = 1.0 / B; - const double V = LAT.NSITE; - const double invV = 1.0 / V; - const double D = LAT.D; - - double invsign = 1.0 / X[SGN]; - - for (int i = 0; i < NACC; ++i) { - X[i] *= invsign; - } - - Q[SIGN] = X[SGN]; - Q[ANV] = X[NV1] * invV; - Q[ENE] = (EBASE + X[EB1] / B) * invV; - - Q[SPE] = (X[EB2] - X[EB1] * X[EB1] - X[NV1]) * invV; - - Q[LEN] = X[LE1]; - Q[XMX] = WDIAG * X[LE1] * T; - - Q[AMZU] = X[MZUA1]; - Q[BMZU] = X[MZUB1]; - Q[SMZU] = (X[MZUA2] - X[MZUA1] * X[MZUA1]) * V; - Q[XMZU] = (X[MZUB2] - X[MZUB1] * X[MZUB1]) * B * V; - - Q[AMZS] = X[MZSA1]; - Q[BMZS] = X[MZSB1]; - Q[SMZS] = (X[MZSA2] - X[MZSA1] * X[MZSA1]) * V; - Q[XMZS] = (X[MZSB2] - X[MZSB1] * X[MZSB1]) * B * V; - - for (int i = 0; i < NPHY; i++) PHY[i].accumulate(Q[i]); -} diff --git a/src/dla/SPECIFIC/Heisenberg/measure_specific.h b/src/dla/SPECIFIC/Heisenberg/measure_specific.h deleted file mode 100644 index f6cdd167..00000000 --- a/src/dla/SPECIFIC/Heisenberg/measure_specific.h +++ /dev/null @@ -1,77 +0,0 @@ -#ifndef SRC_DLA_SPECIFIC_HEISENBERG_MEASURE_SPECIFIC_H_ -#define SRC_DLA_SPECIFIC_HEISENBERG_MEASURE_SPECIFIC_H_ - -#include - -namespace Specific { - -std::vector diagonal_operators(int NXMAX) { - std::vector ret(NXMAX); - for (int i = 0; i < NXMAX; ++i) { - ret[i] = i - 0.5 * (NXMAX - 1); - } - return ret; -} - -// accumulator specifier -enum { - SGN, // Sign of the weight - NV1, // Number of Vertices - EB1, // Number of kinks - EB2, // Number of kinks^2 - LE1, // LEngth of worm pair - - // observables about Magnetization Z per site - // A: tau=0 B: tau-integrated - // U: Uniform S: Staggered - // 1: linear 2: squared - MZUA1, - MZUA2, - MZUB1, - MZUB2, - MZSA1, - MZSA2, - MZSB1, - MZSB2, - - NACC, // number of quantities measured at each MC step -}; -static const char ANAME[][NACC] = {"sgn", "nv1", "eb1", "eb2", "le1", - "mzua1", "mzua2", "mzub1", "mzub2", "mzsa1", - "mzsa2", "mzsb1", "mzsb2"}; - -// observable specifier -enum { - SIGN, // Sign of weight - - ANV, // (site) Average of the Number of Vertices - ENE, // ENErgy per site - SPE, // SPEcific heat - LEN, // LENgth of worm pair - XMX, - - // observables about Magnetization Z per site - // prefix: - // A: tau=0 B: tau-integrated - // S: Variance X: Susceptibility - // postfix: - // U: Uniform, S: Staggered - AMZU, - BMZU, - SMZU, - XMZU, - AMZS, - BMZS, - SMZS, - XMZS, - - NPHY, // number of quantities computed at each set -}; -static const char PNAME[][NPHY] = { - "sign", "anv", "ene", "spe", "len", "xmx", "amzu", - "bmzu", "smzu", "xmzu", "amzs", "bmzs", "smzs", "xmzs", -}; - -} // namespace Specific - -#endif // SRC_DLA_SPECIFIC_HEISENBERG_MEASURE_SPECIFIC_H_ diff --git a/src/dla/algorithm.hpp b/src/dla/algorithm.hpp index 60d69403..c6747a1e 100644 --- a/src/dla/algorithm.hpp +++ b/src/dla/algorithm.hpp @@ -14,23 +14,23 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -//###################################################################### -//#### -//#### World-Line Monte Carlo simulation -//#### by the Directed-Loop Algorithm -//#### -//#### Mar.03 / 2005, Naoki Kawashima -//#### -//###################################################################### - -//###################################################################### -//#### -//#### World-Line Monte Carlo simulation -//#### by the non-Vertex Directed-Loop Algorithm -//#### -//#### Nov.11 / 2007, Yasuyuki Kato -//#### -//###################################################################### +// ###################################################################### +// #### +// #### World-Line Monte Carlo simulation +// #### by the Directed-Loop Algorithm +// #### +// #### Mar.03 / 2005, Naoki Kawashima +// #### +// ###################################################################### + +// ###################################################################### +// #### +// #### World-Line Monte Carlo simulation +// #### by the non-Vertex Directed-Loop Algorithm +// #### +// #### Nov.11 / 2007, Yasuyuki Kato +// #### +// ###################################################################### #ifndef SRC_DLA_ALGORITHM_HPP_ #define SRC_DLA_ALGORITHM_HPP_ @@ -55,9 +55,9 @@ class VertexProperty; class VertexInitialConfiguration; class ScatteringChannel; -//###################################################################### -//#### Object Property Declarations -//###################################################################### +// ###################################################################### +// #### Object Property Declarations +// ###################################################################### class SiteInitialConfiguration { public: @@ -224,13 +224,13 @@ class Algorithm { int NVTYPE; // Number of bond types int NXMAX; // Maximum number of segment states - // double WDIAG; // Artifitial weight attached to the diagonal state + // double WDIAG; // Artificial weight attached to the diagonal state // Used in relating the worm-head mean-path to the susceptibility void read(); void initialize(); - int MXNIC1; // Maximum number of Vertex InitialConfigureations + int MXNIC1; // Maximum number of Vertex InitialConfigurations void set_i(int n) { ix = n; } void set_alg(int n); @@ -241,9 +241,9 @@ class Algorithm { void dump() { X.dump(); } }; -//###################################################################### -//########### Member Functions ####################################### -//###################################################################### +// ###################################################################### +// ########### Member Functions ####################################### +// ###################################################################### inline Algorithm::Algorithm(std::string const& FNAME) { AutoDebugDump("Algorithm::Algorithm"); @@ -304,9 +304,9 @@ void Algorithm::read() { if (name == "Site") { int id = B["STYPE"].getInteger(); // id=0 - SPROP(id).initialize( - B); // SitePropertyClass setting is finished at this point - // SPROP(0) returns SiteProperty& val[0] + // SitePropertyClass setting is finished at this point + // note: SPROP(0) returns SiteProperty& val[0] + SPROP(id).initialize(B); } if (name == "Interaction") { @@ -333,11 +333,11 @@ void Algorithm::initialize() { // scattering prob for (int i = 0; i < NSTYPE; i++) { // 1 SiteProperty& SP = SPROP(i); - for (int j = 0; j < SP.NIC; j++) { // NIC=NX=Numberofstates 2 + for (int j = 0; j < SP.NIC; j++) { // NIC=NX=Number of states 2 SiteInitialConfiguration& IC = SP.IC[j]; // operator return val[i] double p = 0.0; - for (int k = 0; k < IC.NCH; k++) { // numberofchannnels + for (int k = 0; k < IC.NCH; k++) { // number of channnels ScatteringChannel& SC = IC.CH[k]; p += SC.PROB; SC.PROB = p; @@ -434,7 +434,7 @@ void Algorithm::initialize() { } double p = 0.0; - bool isVertex = false; // false: kink or term , true: means interaction + bool isVertex = false; // false: kink or term, true: means interaction int i_type = 0; for (i_type = 0; i_type < NITYPE; i_type++) { if (IPROP[i_type].VTYPE == i) { @@ -445,7 +445,6 @@ void Algorithm::initialize() { // InteractionProperty with density of vertex. if (isKink || (!isVertex)) { - // printf("i=%d j=%d iskink \n",i,j); IC.dRHO = 1.0; for (int k = 0; k < IC.NCH; k++) { ScatteringChannel& SC = IC.CH[k]; @@ -457,7 +456,6 @@ void Algorithm::initialize() { double Vdensity = IPROP[i_type].VertexDensity(x); int count_ch = 0; double dR = Vdensity; - // printf("i=%d j=%d dR=%f \n",i,j,dR); for (int k = 0; k < IC.NCH; k++) { ScatteringChannel& SC = IC.CH[k]; diff --git a/src/dla/array.h b/src/dla/array.h index f66365d7..68013197 100644 --- a/src/dla/array.h +++ b/src/dla/array.h @@ -23,28 +23,28 @@ class IndexSystem { public: void init(const int d, const int* l, const std::string& LBL0 = ""); - IndexSystem() { INI = false; }; + IndexSystem() { INI = false; } IndexSystem(const int d, const int* l, const std::string& LBL0 = "") { INI = false; init(d, l, LBL0); - }; + } IndexSystem operator=(IndexSystem& I) { printf("IndexSystem::=()> Error.\n"); printf(" ... The copy operation is prohibited.\n"); exit(0); - }; + } IndexSystem(IndexSystem& I) { printf("IndexSystem::IndexSystem> Error.\n"); printf(" ... The copy operation is prohibited.\n"); exit(0); - }; + } ~IndexSystem() { if (initialized()) delete[] L; - }; + } bool initialized() const { return INI; }; int coord(const int ist, const int d); @@ -55,21 +55,21 @@ class IndexSystem { exit(0); } return D; - }; + } int size() { if (!initialized()) { printf("IndexSystem::size()> Error. Not yet initialized.\n"); exit(0); } return N; - }; + } int size(int i) { if (!initialized()) { printf("IndexSystem::size(int)> Error. Not yet initialized.\n"); exit(0); } return L[i]; - }; + } int operator()(const int* x) const; int operator()(std::vector const& x) const; int operator()(const int M, va_list& ap) const; @@ -87,7 +87,7 @@ class IndexSystem { printf(" %d", L[i]); } printf(") = %d\n", N); - }; + } }; void IndexSystem::init(const int d, const int* l, const std::string& LBL0) { @@ -109,7 +109,7 @@ void IndexSystem::init(const int d, const int* l, const std::string& LBL0) { printf("IndexSystem::init> Error. N = 0.\n"); for (int i = 0; i < d; i++) printf(" L[%d] = %d\n", i, L[i]); exit(0); - }; + } } int IndexSystem::coord(const int ist, const int d) { @@ -209,7 +209,6 @@ int IndexSystem::operator()(const int M, va_list& ap) const { exit(0); } int x[D]; - // int* x = new int[D]; //edit sakakura x[0] = M; for (int i = 1; i < D; i++) x[i] = va_arg(ap, int); return (*this)(x); @@ -245,16 +244,15 @@ class Array { Array() { LBL = ""; val = 0; - }; + } Array(const std::string& LBL0) { LBL = LBL0; val = 0; - }; - // Array(int d, ...); + } ~Array(); void set_all(C v); void set_value(double* v); - void setLabel(const char* label) { LBL = label; }; + void setLabel(const char* label) { LBL = label; } int index(int i, int d); C& operator()(const int M, ...); C& operator()(int* x); diff --git a/src/dla/chainjob.hpp b/src/dla/chainjob.hpp index 5f60872e..6b80c2d5 100644 --- a/src/dla/chainjob.hpp +++ b/src/dla/chainjob.hpp @@ -87,8 +87,8 @@ void Simulation::save() { end_job(); } - int NVER_ = - P.NVERMAX - TheVertexPool.count() - LAT.NSITE - 1; // -1 means warm's tail + // -1 means warm's tail + int NVER_ = P.NVERMAX - TheVertexPool.count() - LAT.NSITE - 1; save(cjobout, NVER_); int NVER_count = 0; for (int interaction_ = 0; interaction_ < LAT.NINT; interaction_++) { @@ -113,7 +113,6 @@ void Simulation::save() { } ++NVER_count; } - // } } if ((NVER_) != NVER_count) { std::cout << "ERROR: We can't keep all vertices." << std::endl; @@ -162,12 +161,12 @@ void Simulation::save() { void Simulation::load() { using Serialize::load; - std::map > - newID2V; // first:new segmentID, second: old vertexIDs - std::map - oldID_S; // first:old segmentID, second:new segment address - std::map - oldID_V; // first:old vertexID, second:new vertex address + // new segmentID -> old vertex IDs + std::map > newID2V; + // old segmentID -> new segment address + std::map oldID_S; + // old vertexID -> new vertex address + std::map oldID_V; load(cjobin, isEnd); if (isEnd) { @@ -205,7 +204,6 @@ void Simulation::load() { oldID_S.insert(std::map::value_type(ID_, (&S))); int VID_0 = load(cjobin); - int VID_1 = load(cjobin); std::pair VID2(VID_0, VID_1); diff --git a/src/dla/displacement.hpp b/src/dla/displacement.hpp index b4f1a355..da2b429a 100644 --- a/src/dla/displacement.hpp +++ b/src/dla/displacement.hpp @@ -38,9 +38,6 @@ struct Displacement { int nkinds; std::vector NR; // the number of pairs with the same dR - // std::vector > IR; // IR[isite][jsite] == disp_index - - // IR[std::make_pair(isite,jsite)] == disp_index typedef boost::unordered_map, int, boost::hash > > IR_type; diff --git a/src/dla/dla.cc b/src/dla/dla.cc index 1f051108..bb0ae227 100644 --- a/src/dla/dla.cc +++ b/src/dla/dla.cc @@ -67,6 +67,10 @@ int main(int argc, char* argv[]) { printf("\n\n>>> The program is being run on DEBUG mode.\n\n\n"); #endif + if (I_PROC == 0) { + printf("This is DSQSS ver.%s\n\n", DSQSS_VERSION); + } + Parameter P(argc, argv); TheSegmentPool.init(P.NSEGMAX); diff --git a/src/dla/generators/CMakeLists.txt b/src/dla/generators/CMakeLists.txt index 2d5c6df0..dbece9b0 100644 --- a/src/dla/generators/CMakeLists.txt +++ b/src/dla/generators/CMakeLists.txt @@ -1,5 +1,4 @@ # include guard -cmake_minimum_required(VERSION 2.8) if(${CMAKE_PROJECT_NAME} STREQUAL "Project") message(FATAL_ERROR "cmake should be executed not for 'src' subdirectory, but for the top directory of DLA.") endif(${CMAKE_PROJECT_NAME} STREQUAL "Project") diff --git a/src/dla/lattice.hpp b/src/dla/lattice.hpp index ffbe6c6e..fe59dcd8 100644 --- a/src/dla/lattice.hpp +++ b/src/dla/lattice.hpp @@ -75,9 +75,9 @@ class Lattice { void initialize(); - Site& S(int i) { return site[i]; }; + Site& S(int i) { return site[i]; } - Interaction& I(int i) { return interaction[i]; }; + Interaction& I(int i) { return interaction[i]; } Edge& EDGE(int i) { return edge[i]; } diff --git a/src/dla/measure_specific.cc b/src/dla/measure_specific.cc deleted file mode 100644 index 36d1719c..00000000 --- a/src/dla/measure_specific.cc +++ /dev/null @@ -1,223 +0,0 @@ -// DSQSS (Discrete Space Quantum Systems Solver) -// Copyright (C) 2018- The University of Tokyo -// -// 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 . - -#if 0 -#if defined(MEASURE_BOSON) && !defined(MEASURE_SPIN) && \ - !defined(MEASURE_USER_SPECIFIC) -#include "SPECIFIC/Boson/measure_specific.cc" -#elif !defined(MEASURE_BOSON) && defined(MEASURE_SPIN) && \ - !defined(MEASURE_USER_SPECIFIC) -#include "SPECIFIC/Heisenberg/measure_specific.cc" -#elif !defined(MEASURE_BOSON) && !defined(MEASURE_SPIN) && \ - defined(MEASURE_USER_SPECIFIC) -#include "SPECIFIC/User/measure_specific.cc" -#else -#error You should define ONE of the following tokens: MEASURE_BOSON, MEASURE_SPIN, or MEASURE_USER_SPECIFIC. -#endif -#endif - -#include "measure_specific.h" -#include "util.hpp" - -void Measurement::measure(double sgn) { - using namespace Specific; - int NV = LAT.countVertices(); - - double MZUA = 0.0; // uniform, tau=0 - double MZUB = 0.0; // uniform, integrated - double MZSA = 0.0; // staggered, tau=0 - double MZSB = 0.0; // staggered, integrated - - const double T = 1.0 / LAT.BETA; - const double invV = 1.0 / LAT.NSITE; - - std::vector phase(LAT.NSMTYPE); - for (int m = 0; m < LAT.NSMTYPE; ++m) { - phase[m] = std::cos(2 * M_PI * m / LAT.NSMTYPE); - } - - for (int s = 0; s < LAT.NSITE; s++) { - Site& SITE = LAT.S(s); - int mt = SITE.getMTYPE(); - double ph = phase[mt]; - Segment& S0 = SITE.first(); - double mz0 = dvals[S0.X()]; - Site::iterator p(SITE); - double mza0 = 0.0; - - while (!(++p).atOrigin()) { - Segment& S = *p; - double mz = dvals[S.X()]; - mza0 += mz * S.length(); - } - - MZUA += mz0; - MZUB += mza0; - MZSA += ph * mz0; - MZSB += ph * mza0; - } - MZUA *= invV; - MZSA *= invV; - MZUB *= invV; - MZSB *= invV; - MZUB *= T; - MZSB *= T; - - double EBSAMP = -(double)NV; - - for (int b = 0; b < LAT.NINT; b++) { - Interaction& I = LAT.I(b); - InteractionProperty& IP = I.property(); - // VertexProperty& VP = IP.getVertexProperty(); - int NBODY = IP.NBODY; - std::vector tau(NBODY); - std::vector x(NBODY); - std::vector p(NBODY); - - for (int i = 0; i < NBODY; i++) { - Site& S = I.site(i); - p[i].init(S); - ++p[i]; - tau[i] = p[i]->topTime(); - x[i] = p[i]->X(); - } - - double t = 0.0; - - while (t < LAT.BETA) { - int it = util::min_index(tau); - EBSAMP -= (tau[it] - t) * IP.VertexDensity(x); - - if (p[it]->top().isTerminal()) break; - t = tau[it]; - ++p[it]; - tau[it] = p[it]->topTime(); - x[it] = p[it]->X(); - } - } - - ACC[SGN].accumulate(sgn); - - ACC[NV1].accumulate(sgn * NV); - - ACC[MZUA1].accumulate(sgn * MZUA); - ACC[MZUA2].accumulate(sgn * MZUA * MZUA); - ACC[MZUB1].accumulate(sgn * MZUB); - ACC[MZUB2].accumulate(sgn * MZUB * MZUB); - ACC[MZSA1].accumulate(sgn * MZSA); - ACC[MZSA2].accumulate(sgn * MZSA * MZSA); - ACC[MZSB1].accumulate(sgn * MZSB); - ACC[MZSB2].accumulate(sgn * MZSB * MZSB); - - ACC[EB1].accumulate(sgn * EBSAMP); - ACC[EB2].accumulate(sgn * EBSAMP * EBSAMP); - ACC[NH1].accumulate(sgn * MZUA * EBSAMP); - - // WINDING NUMBER - const int Dim = LAT.BD; - if (Dim > 0) { - std::vector wind(Dim); - - for (int xi = 0; xi < LAT.NEDGE; xi++) { - int Asite = LAT.EDGE(xi).A; - int Bsite = LAT.EDGE(xi).B; - int dim = LAT.EDGE(xi).bd; - - Site& SITE = LAT.S(Asite); - Segment& S0 = SITE.first(); - Site::iterator p(SITE); - int xlast = S0.X(); - - // count winding - - while (!(++p).atOrigin()) { - Segment& S = *p; - if (xlast != S.X()) { - if ((*S.bottom().S(0).getONSITE()).id() - 1 == Bsite) { - wind[dim] += (double)(S.bottom().S(1).X() - S.bottom().S(0).X()); - } else if ((*S.bottom().S(2).getONSITE()).id() - 1 == Bsite) { - wind[dim] += (double)(S.bottom().S(3).X() - S.bottom().S(2).X()); - } - xlast = S.X(); - } - } - } // end of for(xi) - - double wind2 = 0.0; - for (int di = 0; di < LAT.D; di++) { - double W = 0; - for (int db = 0; db < Dim; db++) { - W += wind[db] * LAT.vec[di][db] * LAT.L[di]; - } - wind2 += W * W; - } - ACC[Wxy2].accumulate(sgn * wind2); - } // end of if(Dim) -} - -void Measurement::setsummary() { - using namespace Specific; - double WDIAG = ALG.getBlock( - "WDIAG", (double)1.0); // ALG.X["General"]["WDIAG" ].getDouble(); // 0.25 - - std::vector X(NACC); - - for (int i = 0; i < NACC; i++) { - ACC[i].average(); - X[i] = ACC[i].mean(); - } - - double B = LAT.BETA; - double T = 1.0 / B; - double V = LAT.NSITE; - double invV = 1.0 / V; - double D = LAT.D; - - double invsign = 1.0 / X[SGN]; - - for (int i = 0; i < NACC; ++i) { - X[SGN] *= invsign; - } - - Q[SIGN] = X[SGN]; - - Q[ANV] = X[NV1] * invV; - Q[ENE] = (EBASE + X[EB1] * T) * invV; - - Q[SPE] = (X[EB2] - X[EB1] * X[EB1] - X[NV1]) * invV; - - Q[LEN] = X[LE1]; - Q[XMX] = WDIAG * X[LE1] / B; - - Q[AMZU] = X[MZUA1]; - Q[BMZU] = X[MZUB1]; - Q[SMZU] = (X[MZUA2] - X[MZUA1] * X[MZUA1]) * V; - Q[XMZU] = (X[MZUB2] - X[MZUB1] * X[MZUB1]) * B * V; - - Q[AMZS] = X[MZSA1]; - Q[BMZS] = X[MZSB1]; - Q[SMZS] = (X[MZSA2] - X[MZSA1] * X[MZSA1]) * V; - Q[XMZS] = (X[MZSB2] - X[MZSB1] * X[MZSB1]) * B * V; - - Q[DS1] = B * (X[NH1] - X[MZUA1] * X[EB1]) / V; - Q[W2] = X[Wxy2]; - Q[RHOS] = X[Wxy2] * 0.5 / D / V / B; - Q[COMP] = Q[XMZU] / (X[MZUB1] * X[MZUB1]); - - for (int i = 0; i < NPHY; i++) { - PHY[i].accumulate(Q[i]); - } -} diff --git a/src/dla/measure_specific.h b/src/dla/measure_specific.h index 1166792f..338f54b5 100644 --- a/src/dla/measure_specific.h +++ b/src/dla/measure_specific.h @@ -14,22 +14,11 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#if 0 -#if defined(MEASURE_BOSON) && !defined(MEASURE_SPIN) && \ - !defined(MEASURE_USER_SPECIFIC) -#include "SPECIFIC/Boson/measure_specific.h" -#elif !defined(MEASURE_BOSON) && defined(MEASURE_SPIN) && \ - !defined(MEASURE_USER_SPECIFIC) -#include "SPECIFIC/Heisenberg/measure_specific.h" -#elif !defined(MEASURE_BOSON) && !defined(MEASURE_SPIN) && \ - defined(MEASURE_USER_SPECIFIC) -#include "SPECIFIC/User/measure_specific.h" -#else -#error You should define ONE of the following tokens: MEASURE_BOSON, MEASURE_SPIN, or MEASURE_USER_SPECIFIC. -#endif -#endif +#ifndef SRC_DLA_MEASURE_SPECIFIC_H +#define SRC_DLA_MEASURE_SPECIFIC_H #include +#include namespace Specific { @@ -92,3 +81,5 @@ static std::string PNAME[NPHY] = {"sign", "anv", "ene", "spe", "som", "len", "wi2", "rhos", "rhof", "comp"}; } // namespace Specific + +#endif // SRC_DLA_MEASURE_SPECIFIC_H diff --git a/src/dla/objects.hpp b/src/dla/objects.hpp index 594aca0f..3096700f 100644 --- a/src/dla/objects.hpp +++ b/src/dla/objects.hpp @@ -115,14 +115,14 @@ class BareVertex { // Used when PlaceWorm void init(double t, VertexProperty& VP) { TIME = t; - _s.init(1, VP.NLEG); // koko !!! + _s.init(1, VP.NLEG); _s.set_all(0); _VP = &VP; } void init(Interaction* O_I, double t, VertexProperty& VP) { TIME = t; - _s.init(1, VP.NLEG); // koko !!! + _s.init(1, VP.NLEG); _s.set_all(0); _VP = &VP; ONINTERACTION = O_I; @@ -342,7 +342,6 @@ class Site : public Ring { } void idclear() { lastID = 0; } - // void idclear() { cout << "dddd " << endl; }; Vertex& getVterm() { return (*_vterm); } }; @@ -457,7 +456,6 @@ inline double BareSegment::length() { return topTime() - bottomTime(); } //====================================================================== inline void BareSegment::dump() { - // printf("\nBareSegment::dump> Start.\n"); fprintf(FERR, " Segment [%3d] : ", id()); double bt = 0.0; int bid = -1; @@ -678,15 +676,7 @@ inline void Vertex::reconnect() { inline void Vertex::erase() { reconnect(); - // if (ALERT) { - // printf("\nVertex::erase> ### Erasing\n"); - // dump(); - //} Linked::remove(); - // if (ALERT) { - // printf("\nVertex::erase> ### Done\n"); - // dump(); - //} TheVertexPool.push(*this); } @@ -760,11 +750,6 @@ inline Site::Site() { //====================================================================== inline Site::~Site() { - /* - for(int iCI = 0; iCI") + sys.exit(1) + +for spdir in site.getsitepackages([sys.argv[1]]): + p = pathlib.Path(spdir) + if p.exists(): + print(p.resolve()) + break diff --git a/tool/cmake/install.sh.in b/tool/cmake/install.sh.in index 70dfd01d..ed7905fa 100644 --- a/tool/cmake/install.sh.in +++ b/tool/cmake/install.sh.in @@ -1,23 +1,18 @@ -export PYTHONPATH=@DSQSS_PYTHONPATH@ +PREFIX=@CMAKE_INSTALL_PREFIX@ cd @CMAKE_CURRENT_BINARY_DIR@ -if ! @PYTHON_EXECUTABLE@ -m pip --version >/dev/null 2>/dev/null ;then - set -e - wget https://bootstrap.pypa.io/get-pip.py -O get-pip.py - PIP_USER=0 @PYTHON_EXECUTABLE@ get-pip.py --prefix=@CMAKE_INSTALL_PREFIX@ -fi +echo "-- Installing dsqss tools to $PREFIX" +@PYTHON_EXECUTABLE@ -m pip install --prefix=$PREFIX --no-deps @CMAKE_CURRENT_SOURCE_DIR@ -set -e -PIP_USER=0 @PYTHON_EXECUTABLE@ -m pip install -U pip --prefix=@CMAKE_INSTALL_PREFIX@ -PIP_USER=0 @PYTHON_EXECUTABLE@ -m pip install wheel --prefix=@CMAKE_INSTALL_PREFIX@ --no-warn-script-location -DSQSS_VERSION=$(@PYTHON_EXECUTABLE@ setup.py --version) -@PYTHON_EXECUTABLE@ setup.py bdist_wheel --universal -mkdir -p temp -cd temp +echo "export PATH=${PREFIX}/bin:\$PATH" > dsqssvars.sh -# To reinstall dsqss forcely, remove old dsqss package -# `pip uninstall dsqss` removes ALL libraries and binaries -rm -rf @DSQSS_PYTHONPATH@/dsqss* +# Some systems like Ubuntu 22.04 install python packages under ${PREFIX}/local instead of ${PREFIX} +if [ -d ${PREFIX}/local ]; then + echo "export PATH=${PREFIX}/local/bin:\$PATH" >> dsqssvars.sh +fi +DSQSS_PYTHONPATH=$(@PYTHON_EXECUTABLE@ @CMAKE_CURRENT_SOURCE_DIR@/cmake/get_site_packages.py $PREFIX) -PIP_USER=0 @PYTHON_EXECUTABLE@ -m pip install ../dist/dsqss-${DSQSS_VERSION}-py2.py3-none-any.whl --prefix=@CMAKE_INSTALL_PREFIX@ --no-warn-script-location +echo "export PYTHONPATH=${DSQSS_PYTHONPATH}:\$PYTHONPATH" >> dsqssvars.sh +mkdir -p ${PREFIX}/share/dsqss +cp dsqssvars.sh ${PREFIX}/share/dsqss/dsqssvars-@DSQSS_VERSION@.sh diff --git a/tool/cmake/pmwa_pre.in b/tool/cmake/pmwa_pre.in index daaf070f..006ea9ed 100755 --- a/tool/cmake/pmwa_pre.in +++ b/tool/cmake/pmwa_pre.in @@ -18,7 +18,6 @@ import sys sys.path.append('@CMAKE_CURRENT_BINARY_DIR@') -sys.path.append('@pythonpath_build@') from dsqss.pmwa_pre import main if __name__ == '__main__' : main() diff --git a/tool/dsqss/__init__.py b/tool/dsqss/__init__.py index e69de29b..b7ca06cb 100644 --- a/tool/dsqss/__init__.py +++ b/tool/dsqss/__init__.py @@ -0,0 +1 @@ +__version__ = "2.1-dev" diff --git a/tool/dsqss/algorithm.py b/tool/dsqss/algorithm.py index db8a734f..8ab34312 100644 --- a/tool/dsqss/algorithm.py +++ b/tool/dsqss/algorithm.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,27 +14,41 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import codecs +from typing import List, Iterable, TextIO, Sequence, Optional, NamedTuple, Tuple + import math -from collections import namedtuple from copy import deepcopy from itertools import chain, product import numpy as np -from .hamiltonian import append_matelem, keystate -from .prob_kernel import suwa_todo -from .util import tagged +import dsqss +from dsqss import prob_kernel +import dsqss.hamiltonian +import dsqss.util +from dsqss.util import tagged + + +class Channel(NamedTuple): + outleg: int + state: int + prob: float -Channel = namedtuple("Channel", ["outleg", "state", "prob"]) -SiteInitialConfiguration = namedtuple("SiteInitialConfiguration", ["state", "channels"]) -InitialConfiguration = namedtuple( - "InitialConfiguration", ["istate", "fstate", "direction", "newstate", "channels"] -) +class SiteInitialConfiguration(NamedTuple): + state: int + channels: List[Channel] -def num_channels(sources, nx): +class InitialConfiguration(NamedTuple): + istate: Tuple[int, ...] + fstate: Tuple[int, ...] + direction: int + newstate: int + channels: List[Channel] + + +def num_channels(sources: dsqss.hamiltonian.MatElems, nx: int) -> np.ndarray: nchs = np.zeros(nx, dtype=int) for source in sources: st = source[0][0] @@ -43,14 +57,20 @@ def num_channels(sources, nx): class AlgSite: - def __init__(self, hamsite, cutoff=1e-10): + id: int + N: int + values: np.ndarray + vtype: int + initialconfigurations: List[SiteInitialConfiguration] + + def __init__(self, hamsite: dsqss.hamiltonian.Site, cutoff=1e-10) -> None: self.id = hamsite.id self.N = hamsite.N - self.values = hamsite.values + self.values = np.array(hamsite.values) self.vtype = self.id self.initialconfigurations = [] for st in range(self.N): - channels = [] + channels: List[Channel] = [] bouncew = 1.0 for source in hamsite.sources: if source[0][0] != st: @@ -58,14 +78,14 @@ def __init__(self, hamsite, cutoff=1e-10): weight = hamsite.sources[source] ** 2 t = source[1][0] bouncew -= weight - channels.append(Channel(0, t, 0.5 * weight)) - channels.append(Channel(1, t, 0.5 * weight)) + channels.append(Channel(outleg=0, state=t, prob=0.5 * weight)) + channels.append(Channel(outleg=1, state=t, prob=0.5 * weight)) if bouncew > cutoff: - channels.append(Channel(-1, -1, bouncew)) + channels.append(Channel(outleg=-1, state=-1, prob=bouncew)) self.initialconfigurations.append(SiteInitialConfiguration(st, channels)) self.vertex = WormVertex(hamsite) - def write_xml(self, f, level): + def write_xml(self, f: TextIO, level: int) -> None: indent = " " f.write(indent * level + "\n") level += 1 @@ -93,13 +113,20 @@ def write_xml(self, f, level): class Vertex(object): - def __init__(self, vtype, category, nbody, ics): + id: int + category: int + nbody: int + initialconfigurations: List[InitialConfiguration] + + def __init__( + self, vtype: int, category: int, nbody: int, ics: Iterable[InitialConfiguration] + ): self.id = vtype self.category = category self.nbody = nbody - self.initialconfigurations = ics + self.initialconfigurations = list(ics) - def write_xml(self, f, level): + def write_xml(self, f: TextIO, level: int) -> None: indent = " " f.write(indent * level + "\n") level += 1 @@ -134,7 +161,12 @@ def write_xml(self, f, level): class WormVertex(Vertex): - def __init__(self, hamsite): + id: int + category: int + nbody: int + initialconfigurations: List[InitialConfiguration] + + def __init__(self, hamsite: dsqss.hamiltonian.Site) -> None: self.id = hamsite.id self.category = 1 self.nbody = 1 @@ -142,30 +174,62 @@ def __init__(self, hamsite): for source in hamsite.sources: self.initialconfigurations.append( InitialConfiguration( - source[0], source[1], 0, source[1][0], [Channel(-1, -1, 1.0)] + source.initial, + source.final, + 0, + source.final[0], + [Channel(outleg=-1, state=-1, prob=1.0)], ) ) self.initialconfigurations.append( InitialConfiguration( - source[0], source[1], 1, source[0][0], [Channel(-1, -1, 1.0)] + source.initial, + source.final, + 1, + source.initial[0], + [Channel(outleg=-1, state=-1, prob=1.0)], ) ) class IntVertex(Vertex): - def __init__(self, vtype, nbody, ics): + id: int + category: int + nbody: int + initialconfigurations: List[InitialConfiguration] + + def __init__( + self, vtype: int, nbody: int, ics: Iterable[InitialConfiguration] + ) -> None: self.id = vtype self.category = 2 self.nbody = nbody - self.initialconfigurations = ics + self.initialconfigurations = list(ics) class AlgInteraction: - def __init__(self, hamint, hamsites, prob_kernel=suwa_todo, ebase_extra=0.0): + itype: int + vtype: int + nbody: int + kernel: prob_kernel.KernelCallBack + hamint: dsqss.hamiltonian.IndeedInteraction + sites: List[dsqss.hamiltonian.Site] + sitesources: List[dsqss.hamiltonian.MatElems] + intelements: dsqss.hamiltonian.MatElems + ebase_negsign: float + signs: dsqss.hamiltonian.MatElems + + def __init__( + self, + hamint: dsqss.hamiltonian.IndeedInteraction, + hamsites: Sequence[dsqss.hamiltonian.Site], + kernel: prob_kernel.KernelCallBack = prob_kernel.suwa_todo, + ebase_extra: float = 0.0, + ): self.itype = hamint.itype self.vtype = hamint.itype + len(hamsites) self.nbody = hamint.nbody - self.prob_kernel = prob_kernel + self.kernel = kernel self.hamint = hamint self.sites = [hamsites[stype] for stype in hamint.stypes] @@ -185,7 +249,7 @@ def __init__(self, hamint, hamsites, prob_kernel=suwa_todo, ebase_extra=0.0): self.ebase_negsign = min(self.ebase_negsign, v) maxd = max(maxd, v) else: - append_matelem( + dsqss.hamiltonian.append_matelem( self.signs, istate=elem[0], fstate=elem[1], value=np.sign(v) ) v = abs(v) @@ -198,12 +262,14 @@ def __init__(self, hamint, hamsites, prob_kernel=suwa_todo, ebase_extra=0.0): ndiag = 0 for st in product(*map(lambda site: range(site.N), self.sites)): ndiag += 1 - k = (st, st) + k = dsqss.hamiltonian.keystate(st, st) if k in self.intelements: self.intelements[k] += self.ebase_negsign sumw += self.intelements[k] else: - append_matelem(self.intelements, state=st, value=self.ebase_negsign) + dsqss.hamiltonian.append_matelem( + self.intelements, state=st, value=self.ebase_negsign + ) sumw += self.intelements[k] self.ebase_nobounce = max((2 * maxo - sumw) / ndiag, 0.0) @@ -213,7 +279,7 @@ def __init__(self, hamint, hamsites, prob_kernel=suwa_todo, ebase_extra=0.0): self.ebase_extra = 0.5 * maxo for st in product(*map(lambda site: range(site.N), self.sites)): - k = (st, st) + k = dsqss.hamiltonian.keystate(st, st) self.intelements[k] += self.ebase_nobounce + self.ebase_extra self.intelements = {k: v for k, v in self.intelements.items() if v > 0.0} @@ -224,13 +290,18 @@ def __init__(self, hamint, hamsites, prob_kernel=suwa_todo, ebase_extra=0.0): insite = inleg // 2 for instate in range(self.sites[insite].N): ic = self.make_initialconfiguration(st, inleg, instate) - if not ic: + if ic is None: continue else: initialconfigurations.append(ic) self.vertex = IntVertex(self.vtype, self.nbody, initialconfigurations) - def make_initialconfiguration(self, states, inleg, instate): + def make_initialconfiguration( + self, + states: dsqss.hamiltonian.States, + inleg: int, + instate: int, + ) -> Optional[InitialConfiguration]: nbody = self.nbody insite = inleg // 2 intau = inleg % 2 @@ -238,22 +309,24 @@ def make_initialconfiguration(self, states, inleg, instate): if intau == 0: p = self.sitesources[insite].get( - keystate(istate=instate, fstate=states[0][insite]), 0.0 + dsqss.hamiltonian.keystate(istate=instate, fstate=states[0][insite]), + 0.0, ) else: p = self.sitesources[insite].get( - keystate(istate=states[1][insite], fstate=instate), 0.0 + dsqss.hamiltonian.keystate(istate=states[1][insite], fstate=instate), + 0.0, ) if p == 0.0: - return False + return None midstates = [list(states[0]), list(states[1])] midstates[intau][insite] = instate + incomeindex = 100000 # very large integer probs = [] outlegs = [] outstates = [] - # import pdb; pdb.set_trace() for outleg in range(2 * nbody): outsite = outleg // 2 outtau = outleg % 2 @@ -261,14 +334,18 @@ def make_initialconfiguration(self, states, inleg, instate): st = deepcopy(midstates) oldoutstate = st[outtau][outsite] st[outtau][outsite] = outstate - p = self.intelements.get(keystate(istate=st[0], fstate=st[1]), 0.0) + p = self.intelements.get( + dsqss.hamiltonian.keystate(istate=st[0], fstate=st[1]), 0.0 + ) if outtau == 0: p *= self.sitesources[outsite].get( - keystate(istate=oldoutstate, fstate=outstate), 0.0 + dsqss.hamiltonian.keystate(istate=oldoutstate, fstate=outstate), + 0.0, ) else: p *= self.sitesources[outsite].get( - keystate(istate=outstate, fstate=oldoutstate), 0.0 + dsqss.hamiltonian.keystate(istate=outstate, fstate=oldoutstate), + 0.0, ) if p == 0.0: continue @@ -278,14 +355,16 @@ def make_initialconfiguration(self, states, inleg, instate): if inleg == outleg and outstate == oldinstate: incomeindex = len(probs) - 1 if len(probs) == 0: - return False - W = self.prob_kernel(probs)[incomeindex, :] + return None + W = self.kernel(probs)[incomeindex, :] channels = [ - Channel(outleg, outstate, p) + Channel(outleg=outleg, state=outstate, prob=p) for outleg, outstate, p in zip(outlegs, outstates, W) if p > 0.0 ] - return InitialConfiguration(states[0], states[1], inleg, instate, channels) + return InitialConfiguration( + states.initial, states.final, inleg, instate, channels + ) def write_xml(self, f, level): indent = " " @@ -317,7 +396,17 @@ def write_xml(self, f, level): class Algorithm: - def __init__(self, ham, prob_kernel=suwa_todo, ebase_extra=0.0): + name: str + nstypes: int + nitypes: int + nxmax: int + + def __init__( + self, + ham: dsqss.hamiltonian.GraphedHamiltonian, + prob_kernel: prob_kernel.KernelCallBack = prob_kernel.suwa_todo, + ebase_extra: float = 0.0, + ): self.name = ham.name self.nstypes = ham.nstypes self.nitypes = ham.nitypes @@ -335,13 +424,13 @@ def __init__(self, ham, prob_kernel=suwa_todo, ebase_extra=0.0): self.sites = [AlgSite(site) for site in hamsites] self.interactions = [ AlgInteraction( - hamint, hamsites, prob_kernel=prob_kernel, ebase_extra=ebase_extra + hamint, hamsites, kernel=prob_kernel, ebase_extra=ebase_extra ) for hamint in ham.indeed_interactions ] - def write_xml(self, filename): - with codecs.open(filename, "w", "utf_8") as f: + def write_xml(self, filename: dsqss.util.Filename) -> None: + with open(filename, "w", encoding="utf_8") as f: indent = " " f.write('\n') f.write("\n") diff --git a/tool/dsqss/bosehubbard.py b/tool/dsqss/bosehubbard.py index 06195223..60e490af 100644 --- a/tool/dsqss/bosehubbard.py +++ b/tool/dsqss/bosehubbard.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,22 +14,24 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from math import sqrt +from typing import MutableMapping -from .hamiltonian import Hamiltonian, Interaction, Site, append_matelem -from .util import extend_list, get_as_list +import numpy as np +from dsqss import util +from dsqss import hamiltonian -def creator_boson(n): - return sqrt(n + 1.0) +def creator_boson(n: float) -> float: + return np.sqrt(n + 1.0) -def annihilator_boson(n): - return sqrt(n) +def annihilator_boson(n: float) -> float: + return np.sqrt(n) -class BoseSite(Site): - def __init__(self, id, M, U, mu): + +class BoseSite(hamiltonian.Site): + def __init__(self, id: int, M: int, U: float, mu: float) -> None: """ M: max n U: n_i (n_i - 1)/2 @@ -39,27 +41,26 @@ def __init__(self, id, M, U, mu): NX = M + 1 values = [] - sources = {} - elements = {} + sources: hamiltonian.MatElems = {} + elements: hamiltonian.MatElems = {} for n in range(NX): value = -mu * n + 0.5 * U * n * (n - 1) values.append(n) - append_matelem(elements, state=n, value=value) + hamiltonian.append_matelem(elements, state=n, value=value) if n > 0: # annihilator - append_matelem( + hamiltonian.append_matelem( sources, istate=n, fstate=n - 1, value=annihilator_boson(n) ) if n < M: # creator value = creator_boson(n) - append_matelem(sources, istate=n, fstate=n + 1, value=value) - super(BoseSite, self).__init__(id=id, N=NX, values=values, - elements=elements, sources=sources) + hamiltonian.append_matelem(sources, istate=n, fstate=n + 1, value=value) + super().__init__(id=id, N=NX, values=values, elements=elements, sources=sources) -class BoseBond(Interaction): - def __init__(self, id, M, t, V): +class BoseBond(hamiltonian.Interaction): + def __init__(self, id: int, M: int, t: float, V: float) -> None: """ M: max N t: c_1 a_2 + c_2 a_1 @@ -74,44 +75,42 @@ def __init__(self, id, M, t, V): c = [creator_boson(n) for n in N] c[-1] = 0.0 a = [annihilator_boson(n) for n in N] - elements = {} + elements: hamiltonian.MatElems = {} for i in range(nx): for j in range(nx): # diagonal w = V * i * j if w != 0.0: - append_matelem(elements, state=[i, j], value=w) + hamiltonian.append_matelem(elements, state=[i, j], value=w) # offdiagonal w = -t * c[i] * a[j] if w != 0.0: - append_matelem( + hamiltonian.append_matelem( elements, istate=[i, j], fstate=[i + 1, j - 1], value=w ) w = -t * a[i] * c[j] if w != 0.0: - append_matelem( + hamiltonian.append_matelem( elements, istate=[i, j], fstate=[i - 1, j + 1], value=w ) - super(BoseBond, self).__init__(id=id, nbody=nbody, - Ns=Ns, elements=elements) + super().__init__(id=id, nbody=nbody, Ns=Ns, elements=elements) -class BoseHubbard_hamiltonian(Hamiltonian): - def __init__(self, param): +class BoseHubbard_hamiltonian(hamiltonian.Hamiltonian): + def __init__(self, param: MutableMapping) -> None: M = param["M"] - Us = get_as_list(param, "U", 0.0) - mus = get_as_list(param, "mu", 0.0) - nstypes = max(len(Us), len(mus)) - extend_list(Us, nstypes) - extend_list(mus, nstypes) - self.sites = [BoseSite(i, M, U, mu) - for i, (U, mu) in enumerate(zip(Us, mus))] - - ts = get_as_list(param, "t", 0.0) - Vs = get_as_list(param, "V", 0.0) - nitypes = max(len(ts), len(Vs)) - extend_list(ts, nitypes) - extend_list(Vs, nitypes) + Us = util.get_as_list(param, "U", 0.0) + mus = util.get_as_list(param, "mu", 0.0) + self.nstypes = max(len(Us), len(mus)) + util.extend_list(Us, self.nstypes) + util.extend_list(mus, self.nstypes) + self.sites = [BoseSite(i, M, U, mu) for i, (U, mu) in enumerate(zip(Us, mus))] + + ts = util.get_as_list(param, "t", 0.0) + Vs = util.get_as_list(param, "V", 0.0) + self.nitypes = max(len(ts), len(Vs)) + util.extend_list(ts, self.nitypes) + util.extend_list(Vs, self.nitypes) self.interactions = [ BoseBond(i, M, t, V) for i, (t, V) in enumerate(zip(ts, Vs)) ] diff --git a/tool/dsqss/displacement.py b/tool/dsqss/displacement.py index 49485074..e4cc0cf9 100644 --- a/tool/dsqss/displacement.py +++ b/tool/dsqss/displacement.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,21 +14,27 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import codecs -from itertools import product +import itertools import numpy as np -from .util import tagged +import dsqss +import dsqss.util +import dsqss.lattice class Displacement: - def __init__(self, lat, distance_only=False): - self.lat = lat - self.nr = 0 + latname: str + nsites: int + nkinds: int + displacements: np.ndarray + + def __init__(self, lat: dsqss.lattice.Lattice, distance_only: bool = False): + self.latname = lat.name + self.nsites = lat.nsites + self.nkinds = 0 sources = range(lat.nsites) - self.displacements = -np.ones((lat.nsites, lat.nsites), dtype=int) - self.nelems = lat.nsites * lat.nsites + self.displacements = -np.ones((lat.nsites, lat.nsites), dtype=np.int64) rdict = {} for s in sources: sr = np.array(lat.sites[s].coord) @@ -44,24 +50,23 @@ def __init__(self, lat, distance_only=False): else: r = tuple(dr) if r not in rdict: - rdict[r] = self.nr - self.nr += 1 + rdict[r] = self.nkinds + self.nkinds += 1 self.displacements[s, t] = rdict[r] - def write_xml(self, filename): - lat = self.lat - with codecs.open(filename, "w", "utf-8") as f: + def write_xml(self, filename: dsqss.util.Filename): + tagged = dsqss.util.tagged + with open(filename, "w", encoding="utf-8") as f: f.write('\n') f.write("\n") - f.write(tagged("Comment", lat.name)) - f.write(tagged("NumberOfKinds", self.nr)) - f.write(tagged("NumberOfSites", lat.nsites)) + f.write(tagged("Comment", self.latname)) + f.write(tagged("NumberOfKinds", self.nkinds)) + f.write(tagged("NumberOfSites", self.nsites)) f.write("\n") f.write("\n") f.write("\n") - for s, t in product(range(lat.nsites), range(lat.nsites)): - f.write(tagged("R", - (self.displacements[s, t], s, t))) + for s, t in itertools.product(range(self.nsites), range(self.nsites)): + f.write(tagged("R", (self.displacements[s, t], s, t))) f.write("\n") diff --git a/tool/dsqss/dla_alg.py b/tool/dsqss/dla_alg.py index 9e8184e0..e86b0f5e 100644 --- a/tool/dsqss/dla_alg.py +++ b/tool/dsqss/dla_alg.py @@ -6,7 +6,7 @@ # 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. +# (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 @@ -16,45 +16,46 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from __future__ import print_function +import sys -import argparse - -from .algorithm import Algorithm -from .displacement import Displacement -from .hamiltonian import GraphedHamiltonian -from .lattice import Lattice -from .prob_kernel import (heat_bath, metropolis, reversible_suwa_todo, - suwa_todo) -from .util import ERROR -from .wavevector import Wavevector +import dsqss +import dsqss.util +import dsqss.algorithm +import dsqss.hamiltonian +import dsqss.displacement +import dsqss.lattice +import dsqss.wavevector +import dsqss.prob_kernel def dla_alg( lat, - ham=None, - wv=None, - latxml=None, - algxml=None, - wvxml=None, - dispxml=None, - prob_kernel=suwa_todo, - ebase_extra=0.0, - distance_only=False, + ham: dsqss.hamiltonian.GraphedHamiltonian = None, + wv: dsqss.wavevector.Wavevector = None, + latxml: str = None, + algxml: str = None, + wvxml: str = None, + dispxml: str = None, + prob_kernel: dsqss.prob_kernel.KernelCallBack = dsqss.prob_kernel.suwa_todo, + ebase_extra: float = 0.0, + distance_only: bool = False, ): if latxml is not None: lat.write_xml(latxml) if ham is not None and algxml is not None: - alg = Algorithm(ham, prob_kernel=prob_kernel, ebase_extra=ebase_extra) + alg = dsqss.algorithm.Algorithm( + ham, prob_kernel=prob_kernel, ebase_extra=ebase_extra + ) alg.write_xml(algxml) if wv is not None and wvxml is not None: wv.write_xml(wvxml, lat) if dispxml is not None: - disp = Displacement(lat, distance_only=distance_only) + disp = dsqss.displacement.Displacement(lat, distance_only=distance_only) disp.write_xml(dispxml) def main(): + import argparse parser = argparse.ArgumentParser( description="Generate algorithm and lattice XML files for DSQSS/DLA", @@ -134,27 +135,36 @@ def main(): default=0.0, help="extra energy shift", ) + parser.add_argument("--version", action="version", version=dsqss.__version__) args = parser.parse_args() if args.kernel == "suwa todo": - kernel = suwa_todo + kernel = dsqss.prob_kernel.suwa_todo elif args.kernel == "reversible suwa todo": - kernel = reversible_suwa_todo + kernel = dsqss.prob_kernel.reversible_suwa_todo elif args.kernel == "heat bath": - kernel = heat_bath + kernel = dsqss.prob_kernel.heat_bath elif args.kernel == "metropolis": - kernel = metropolis + kernel = dsqss.prob_kernel.metropolis elif args.kernel == "metropolice": - ERROR('kenel = "metropolice" is now an invalid option because this is a typographic error in old DSQSS. Use "metropolis".') + dsqss.util.ERROR( + 'kernel = "metropolice" is now an invalid option' + ' because this is a typographic error in old DSQSS.' + ' Use "metropolis".' + ) + sys.exit(1) else: - ERROR("unknown kernel: {0}".format(args.kernel)) + dsqss.util.ERROR(f"unknown kernel: {args.kernel}") + sys.exit(1) - lat = Lattice(args.lat) - ham = GraphedHamiltonian(args.ham, lat) if not args.woalg else None + lat = dsqss.lattice.Lattice(args.lat) + ham = ( + dsqss.hamiltonian.GraphedHamiltonian(args.ham, lat) if not args.woalg else None + ) wv = None if args.kpoint is not None: - wv = Wavevector() + wv = dsqss.wavevector.Wavevector() wv.load(args.kpoint) dla_alg( @@ -170,21 +180,6 @@ def main(): distance_only=args.distance_only, ) - # lat = Lattice(args.lat) - # if not args.wolat: - # lat.write_xml(args.latxml) - # if not args.woalg: - # ham = GraphedHamiltonian(args.ham, lat) - # alg = Algorithm(ham, prob_kernel=kernel, ebase_extra=args.extra_shift) - # alg.write_xml(args.algxml) - # if args.kpoint is not None: - # wv = Wavevector() - # wv.load(args.kpoint) - # wv.write_xml(args.wv, lat) - # if args.disp is not None: - # disp = Displacement(lat, distance_only=args.distance_only) - # disp.write_xml(args.disp) - if __name__ == "__main__": main() diff --git a/tool/dsqss/dla_pre.py b/tool/dsqss/dla_pre.py index b1bb1482..e51d5c40 100644 --- a/tool/dsqss/dla_pre.py +++ b/tool/dsqss/dla_pre.py @@ -6,7 +6,7 @@ # 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. +# (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 @@ -16,61 +16,66 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from __future__ import print_function +import sys -import codecs +from dsqss import ( + algorithm, + displacement, + hamiltonian, + parameter, + prob_kernel, + std_lattice, + std_model, + util, + wavevector, + __version__, +) -from .algorithm import Algorithm -from .displacement import Displacement -from .hamiltonian import GraphedHamiltonian -from .lattice import Lattice -from .parameter import Parameter -from .prob_kernel import (heat_bath, metropolis, reversible_suwa_todo, - suwa_todo) -from .std_lattice import std_lattice -from .std_model import std_model -from .util import ERROR, INFO -from .wavevector import Wavevector +def dla_pre(param, pfile: str) -> None: -def dla_pre(param, pfile): + p = parameter.Parameter(param) - p = Parameter(param) - - lat = std_lattice(param) + lat, _ = std_lattice.std_lattice(param) lat.write_xml(p["latfile"]) - h = std_model(param["hamiltonian"]) - ham = GraphedHamiltonian(h, lat) + h = std_model.std_model(param["hamiltonian"]) + ham = hamiltonian.GraphedHamiltonian(h, lat) palg = param.get("algorithm", {}) extra_shift = palg.get("extra_shift", 0.0) kernel = palg.get("kernel", "suwa todo") if kernel == "suwa todo": - kernel = suwa_todo + kernel = prob_kernel.suwa_todo elif kernel == "reversible suwa todo": - kernel = reversible_suwa_todo + kernel = prob_kernel.reversible_suwa_todo elif kernel == "heat bath": - kernel = heat_bath + kernel = prob_kernel.heat_bath elif kernel == "metropolis": - kernel = metropolis - elif args.kernel == "metropolice": - ERROR('kenel = "metropolice" is now an invalid option because this is a typographic error in old DSQSS. Use "metropolis".') + kernel = prob_kernel.metropolis + elif kernel == "metropolice": + util.ERROR( + 'kernel = "metropolice" is now an invalid option' + " because this is a typographic error in old DSQSS." + ' Use "metropolis".' + ) + sys.exit(1) else: - ERROR("unknown kernel: {0}".format(kernel)) + util.ERROR("unknown kernel: {0}".format(kernel)) + sys.exit(1) - alg = Algorithm(ham, prob_kernel=kernel, ebase_extra=extra_shift) + alg = algorithm.Algorithm(ham, prob_kernel=kernel, ebase_extra=extra_shift) alg.write_xml(p["algfile"]) if p["wvfile"] != "": - wv = Wavevector() + wv = wavevector.Wavevector() wv.generate(param.get("kpoints", {}), size=lat.size) wv.write_xml(p["wvfile"], lat) if p["dispfile"] != "": pdisp = param.get("displacement", {}) - disp = Displacement(lat, pdisp.get("distance_only", False)) + disp = displacement.Displacement(lat, pdisp.get("distance_only", False)) disp.write_xml(p["dispfile"]) p.write_param(pfile) @@ -95,10 +100,11 @@ def main(): parser.add_argument( "-p", "--paramfile", dest="pfile", default="param.in", help="Parameter file" ) + parser.add_argument("--version", action="version", version=__version__) args = parser.parse_args() if args.input is sys.stdin: - INFO("waiting for standard input...") + util.INFO("waiting for standard input...") dla_pre(toml.load(args.input), args.pfile) diff --git a/tool/dsqss/hamiltonian.py b/tool/dsqss/hamiltonian.py index ad6b9de2..6c532120 100644 --- a/tool/dsqss/hamiltonian.py +++ b/tool/dsqss/hamiltonian.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,50 +14,80 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import codecs +import typing +from typing import ( + Union, + Iterable, + List, + Tuple, + Optional, + Sequence, + MutableMapping, + Dict, + NamedTuple, +) + import itertools from copy import deepcopy import toml -from .lattice import Lattice -from .util import tagged +import dsqss +import dsqss.lattice +import dsqss.util -def keystate(istate, fstate): - if type(istate) is not tuple: - if type(istate) is not list: - istate = [istate] - ti = tuple(istate) - else: - ti = istate +State = Union[int, Iterable[int]] - if type(fstate) is not tuple: - if type(fstate) is not list: - fstate = [fstate] - tf = tuple(fstate) - else: - tf = fstate - return (ti, tf) + +class States(NamedTuple): + initial: Tuple[int, ...] + final: Tuple[int, ...] + + +MatElems = Dict[States, float] + + +def keystate(istate: State, fstate: State) -> States: + ret: List[Tuple[int, ...]] = [] + for st in (istate, fstate): + if isinstance(st, int): + ret.append((st,)) + else: + ret.append(tuple(st)) + return States(initial=ret[0], final=ret[1]) def append_matelem( - matelems, state=None, istate=None, fstate=None, value=None, param=None -): + matelems: MatElems, + state: Optional[State] = None, + istate: Optional[State] = None, + fstate: Optional[State] = None, + value: Optional[float] = None, + param: Optional[MutableMapping] = None, +) -> MatElems: if param is None: - if istate is not None: - matelems[keystate(istate, fstate)] = value + if value is None: + raise RuntimeError("value is required") + if istate is not None and fstate is not None: + matelems[keystate(istate, fstate)] = float(value) + elif state is not None: + matelems[keystate(state, state)] = float(value) else: - matelems[keystate(state, state)] = value + raise RuntimeError("(istate and fstate) or state is required") else: - if "istate" in param: - matelems[keystate(param["istate"], param["fstate"])] = param["value"] + if "value" not in param: + raise RuntimeError("value is required") + if "istate" in param and "fstate" in param: + matelems[keystate(param["istate"], param["fstate"])] = float(param["value"]) + elif "state" in param: + matelems[keystate(param["state"], param["state"])] = float(param["value"]) else: - matelems[keystate(param["state"], param["state"])] = param["value"] + raise RuntimeError("(istate and fstate) or state is required") return matelems -def matelems_todict(matelems, keepzero=False): +def matelems_todict(matelems: MatElems, keepzero: bool = False) -> List[Dict]: return [ {"istate": list(key[0]), "fstate": list(key[1]), "value": value} for key, value in matelems.items() @@ -66,9 +96,21 @@ def matelems_todict(matelems, keepzero=False): class Site(object): - def __init__(self, param=None, - id=None, N=None, values=None, - elements=None, sources=None): + id: int + N: int + values: List + elements: MatElems + sources: MatElems + + def __init__( + self, + param: Optional[MutableMapping] = None, + id=None, + N=None, + values=None, + elements=None, + sources=None, + ): if param is not None: self.id = param["type"] self.N = param["N"] @@ -86,7 +128,7 @@ def __init__(self, param=None, self.elements = elements self.sources = sources - def to_dict(self): + def to_dict(self) -> Dict: return { "type": self.id, "N": self.N, @@ -97,7 +139,19 @@ def to_dict(self): class Interaction(object): - def __init__(self, param=None, id=None, nbody=None, Ns=None, elements=None): + id: int + nbody: int + Ns: List[int] + elements: MatElems + + def __init__( + self, + param: Optional[MutableMapping] = None, + id=None, + nbody=None, + Ns=None, + elements=None, + ): if param is not None: self.id = param["type"] self.nbody = param["nbody"] @@ -111,7 +165,7 @@ def __init__(self, param=None, id=None, nbody=None, Ns=None, elements=None): self.Ns = Ns self.elements = elements - def to_dict(self): + def to_dict(self) -> Dict: return { "type": self.id, "nbody": self.nbody, @@ -121,7 +175,17 @@ def to_dict(self): class IndeedInteraction(object): - def __init__(self, sites, ints, v): + itype: int + stypes: List[int] + nbody: int + elements: MatElems + + def __init__( + self, + sites: Sequence[Site], + ints: Sequence[Interaction], + v: dsqss.lattice.Vertex, + ): inter = ints[v.int_type] zs = v.zs @@ -143,59 +207,76 @@ def __init__(self, sites, ints, v): self.elements[key] = val -class Hamiltonian(object): - def __init__(self, ham_dict): - if type(ham_dict) == str: +class Hamiltonian: + name: str + nstypes: int + sites: List[Site] + nitypes: int + interactions: List[Interaction] + + def __init__(self, ham_dict: Union[str, MutableMapping]) -> None: + if isinstance(ham_dict, str): ham_dict = toml.load(ham_dict) self.name = ham_dict.get("name", "") self.nstypes = len(ham_dict["sites"]) - self.sites = [None for i in range(self.nstypes)] - self.nitypes = len(ham_dict["interactions"]) - self.interactions = [None for i in range(self.nitypes)] - + sites: List[Optional[Site]] = [None for i in range(self.nstypes)] for site in ham_dict["sites"]: S = Site(site) - self.sites[S.id] = S + sites[S.id] = S + for site in sites: + if site is None: + raise RuntimeError("") + self.sites = typing.cast(List[Site], sites) + + self.nitypes = len(ham_dict["interactions"]) + interactions: List[Optional[Interaction]] = [None for i in range(self.nitypes)] for inter in ham_dict["interactions"]: Int = Interaction(inter) - self.interactions[Int.id] = Int + interactions[Int.id] = Int + + for inter in interactions: + if inter is None: + raise RuntimeError("") + self.interactions = typing.cast(List[Interaction], interactions) self.nxmax = 0 for site in self.sites: self.nxmax = max(site.N, self.nxmax) - def to_dict(self): + def to_dict(self) -> Dict: return { "name": self.name, "sites": list(map(lambda x: x.to_dict(), self.sites)), "interactions": list(map(lambda x: x.to_dict(), self.interactions)), } - def write_toml(self, filename): - with codecs.open(filename, "w", "utf_8") as f: + def write_toml(self, filename: dsqss.util.Filename) -> None: + with open(filename, "w", encoding="utf_8") as f: toml.dump(self.to_dict(), f) class GraphedHamiltonian(Hamiltonian): - def __init__(self, ham, lat): - if not isinstance(ham, Hamiltonian): - super(GraphedHamiltonian, self).__init__(ham) + indeed_interactions: List[IndeedInteraction] + + def __init__(self, ham, lat) -> None: + if isinstance(ham, Hamiltonian): + super().__init__(ham.to_dict()) else: - super(GraphedHamiltonian, self).__init__(ham.to_dict()) + super().__init__(ham) self.load_lattice(lat) - def load_lattice(self, lat): + def load_lattice(self, lat) -> None: if type(lat) == str: - lat = Lattice(lat) + lat = dsqss.lattice.Lattice(lat) self.indeed_interactions = [ IndeedInteraction(self.sites, self.interactions, v) for v in lat.vertices ] self.nitypes = len(self.indeed_interactions) - def write_xml(self, filename): - with codecs.open(filename, "w", "utf_8") as f: - + def write_xml(self, filename: dsqss.util.Filename) -> None: + tagged = dsqss.util.tagged + with open(filename, "w", encoding="utf_8") as f: f.write('\n') f.write("\n") indent = " " diff --git a/tool/dsqss/lattice.py b/tool/dsqss/lattice.py index c4b35e3d..e5e6f111 100644 --- a/tool/dsqss/lattice.py +++ b/tool/dsqss/lattice.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,9 +14,9 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from __future__ import print_function +import typing +from typing import List, MutableMapping, Tuple, Dict, Iterable -import codecs import sys from itertools import chain @@ -24,94 +24,170 @@ import numpy as np -from .util import ERROR, WARN, coord2index, get_as_list, index2coord, tagged +from dsqss import util +from dsqss.util import tagged class Site: - def __init__(self, site_id, site_type, coordinate): + """Site class + + Attributes: + id: index of site + stype: site type + coord: coordination + z: coordination number (the number of the nearest neighbor sites) + """ + + id: int + stype: int + coord: np.ndarray + z: int + + def __init__(self, site_id: int, site_type: int, coordinate: np.ndarray) -> None: self.id = site_id self.stype = site_type self.coord = coordinate self.z = 0 - def __str__(self): - s = "{0} {1}".format(self.id, self.stype) + def __str__(self) -> str: + s = f"{self.id} {self.stype}" for c in self.coord: - s += " {0}".format(c) + s += f" {c}" return s class Interaction: - def __init__(self, int_id, int_type, nbody, site_indices, edge_flag, direction): + """Interaction class + + Attributes: + id: index of interaction + itype: interaction type + itype_org: + nbody: the number of involved sites + sites: list of index of involved sites + edge: on boundary (1) or not (0) + direction: index of bond direction + """ + + id: int + itype: int + itype_org: int + nbody: int + sites: List[int] + edge: int + direction: int + + def __init__( + self, + int_id: int, + int_type: int, + nbody: int, + site_indices: Iterable[int], + edge: int, + direction: int, + ) -> None: self.id = int_id self.itype = int_type self.itype_org = self.itype self.nbody = nbody - self.sites = site_indices - self.edge = edge_flag + self.sites = list(site_indices) + self.edge = edge self.dir = direction self.vtype = -1 - def __str__(self): - s = "{0} {1} {2}".format(self.id, self.itype_org, self.nbody) + def __str__(self) -> str: + s = f"{self.id} {self.itype_org} {self.nbody}" edgeflag = 1 if self.edge >= 0 else 0 for site in self.sites: - s += " {0}".format(site) - s += " {0} {1}".format(edgeflag, self.dir) + s += f" {site}" + s += f" {edgeflag} {self.dir}" return s class Vertex: - def __init__(self, v_id, int_type, site_types, zs): + """Vertex (sites and interaction) + + Attributes: + v_id: index of vertex + int_type: type of interaction + site_types: type of sites + zs: coodination numbers + + """ + + v_id: int + int_type: int + site_types: List[int] + zs: List[int] + + def __init__( + self, v_id: int, int_type: int, site_types: Iterable[int], zs: Iterable[int] + ) -> None: self.v_id = v_id self.int_type = int_type - self.site_types = site_types - self.zs = zs + self.site_types = list(site_types) + self.zs = list(zs) class Lattice: - def __init__(self, inp=None): + name: str + dim: int + size: List[int] + bc: List[bool] # periodic if True, open if False + latvec: np.ndarray + ndir: int + dirs: List[np.ndarray] + sites: List[Site] + ints: List[Interaction] + nstypes: int + nedges: int + nitypes: int + nvtypes: int + vertices: List[Vertex] + + def __init__(self, inp=None) -> None: if inp is not None: self.load(inp) - def load(self, inp): - if type(inp) is dict: + def load(self, inp) -> None: + if isinstance(inp, dict): self.load_dict(inp) - elif inp.endswith(".toml"): + elif isinstance(inp, str) and inp.endswith(".toml"): self.load_dict(toml.load(inp)) else: self.load_dat(inp) - def load_dict(self, param): + def load_dict(self, param: MutableMapping) -> None: if "lattice" in param and not isinstance(param["lattice"], str): - return self.load_dict(param["lattice"]) + self.load_dict(param["lattice"]) + return parameter = param["parameter"] self.name = parameter["name"] self.dim = parameter["dim"] - self.size = get_as_list(parameter, "L", extendto=self.dim) - self.bc = get_as_list(parameter, "bc", default=True, extendto=self.dim) - ncells = np.product(self.size) + self.size = util.get_as_list(parameter, "L", extendto=self.dim) + self.bc = util.get_as_list(parameter, "bc", default=True, extendto=self.dim) self.latvec = np.array(parameter["basis"], dtype=float).transpose() unitcell = param["unitcell"] nsites_cell = len(unitcell["sites"]) - nbonds_cell = len(unitcell["bonds"]) - localsitecoords = [] + localsitecoords: List[np.ndarray] = [] for site in unitcell["sites"]: localsitecoords.append(np.array(site["coord"], dtype=float)) - directionmap = {} - directions = [] + directionmap: Dict[Tuple, int] = {} + directions: List[int] = [] self.dirs = [] for ib, bond in enumerate(unitcell["bonds"]): if "offset" in bond["source"]: if not all(bond["source"]["offset"] == 0): - ERROR("bonds['source'] has nonzero offset") + util.ERROR("bonds['source'] has nonzero offset") offset = np.array( - get_as_list(bond["target"], "offset", default=0.0, extendto=self.dim), + util.get_as_list( + bond["target"], "offset", default=0.0, extendto=self.dim + ), dtype=float, ) offset += localsitecoords[bond["target"]["siteid"]] @@ -123,45 +199,12 @@ def load_dict(self, param): directions.append(directionmap[key]) self.ndir = len(self.dirs) - def bondsites(center_cellcoord, bond): - ret = [] - offset = center_cellcoord + np.array( - get_as_list(bond["target"], "offset", default=0, extendto=self.dim) - ) - edge = 0 - for dim in range(self.dim): - if offset[dim] < 0: - if self.bc[dim]: - edge = 1 - else: - return False, 0 - elif offset[dim] >= self.size[dim]: - if self.bc[dim]: - edge = 1 - else: - return False, 0 - - for site in (bond["source"], bond["target"]): - cell_coord = center_cellcoord + np.array( - get_as_list(site, "offset", default=0, extendto=self.dim) - ) - for dim in range(self.dim): - if cell_coord[dim] < 0: - cell_coord[dim] += self.size[dim] - elif cell_coord[dim] >= self.size[dim]: - cell_coord[dim] -= self.size[dim] - icell = coord2index(cell_coord, self.size) - sid = site["siteid"] + icell * nsites_cell - ret.append(sid) - if ret[0] == ret[1]: - WARN("selfloop {0}=>{0} appears.".format(ret[0])) - return ret, edge - self.sites = [] self.ints = [] + ncells = np.product(self.size) bid = 0 for icell in range(ncells): - cell_coord = np.array(index2coord(icell, self.size)) + cell_coord = np.array(util.index2coord(icell, self.size)) for lid, site in enumerate(unitcell["sites"]): sid = site["siteid"] + icell * nsites_cell coord = cell_coord + np.array(site["coord"], dtype=float) @@ -169,70 +212,104 @@ def bondsites(center_cellcoord, bond): self.sites.append(S) for ib, bond in enumerate(unitcell["bonds"]): nbody = 2 - sites, edge = bondsites(cell_coord, bond) - if not sites: + sites, edge = self._bondsites(cell_coord, bond, nsites_cell) + if len(sites) == 0: continue INT = Interaction(bid, bond["type"], nbody, sites, edge, directions[ib]) self.ints.append(INT) bid += 1 self.nsites = len(self.sites) self.nints = len(self.ints) - self.update() - - def save_dat(self, out): - if type(out) is str: - with codecs.open(out, "w", "utf-8") as f: - self.save_dat(f) - return + self.check_all() + self._update() + + def _bondsites( + self, + center_cellcoord: np.ndarray, + bond: Dict, + nsites_cell: int, + ) -> Tuple[List[int], int]: + offset = center_cellcoord + np.array( + util.get_as_list(bond["target"], "offset", default=0, extendto=self.dim) + ) + edge = 0 + for dim in range(self.dim): + if offset[dim] < 0: + if self.bc[dim]: + edge = 1 + else: + return [], 0 + elif offset[dim] >= self.size[dim]: + if self.bc[dim]: + edge = 1 + else: + return [], 0 + ret = [] + for site in (bond["source"], bond["target"]): + cell_coord = center_cellcoord + np.array( + util.get_as_list(site, "offset", default=0, extendto=self.dim) + ) + for dim in range(self.dim): + if cell_coord[dim] < 0: + cell_coord[dim] += self.size[dim] + elif cell_coord[dim] >= self.size[dim]: + cell_coord[dim] -= self.size[dim] + icell = util.coord2index(cell_coord, self.size) + sid = site["siteid"] + icell * nsites_cell + ret.append(sid) + if ret[0] == ret[1]: + util.WARN("selfloop {0}=>{0} appears.".format(ret[0])) + return ret, edge + + def save_dat(self, filename: util.Filename) -> None: + out = open(filename, "w", encoding="utf-8") out.write("name\n") - out.write("{0}\n".format(self.name)) + out.write(f"{self.name}\n") out.write("\n") out.write("lattice\n") - out.write("{0} # dim\n".format(self.dim)) + out.write(f"{self.dim} # dim\n") for x in self.size: - out.write("{0} ".format(x)) + out.write(f"{x} ") out.write("# size\n") for x in self.bc: - b = 1 if x else 0 - out.write("{0} ".format(b)) + out.write(f"{int(x)} ") out.write("# 0:open boundary, 1:periodic boundary\n") for d in range(self.dim): - out.write("{0} ".format(d)) + out.write(f"{d} ") for x in self.latvec[:, d]: - out.write("{0} ".format(x)) - out.write("# latvec_{0}\n".format(d)) + out.write(f"{x} ") + out.write(f"# latvec_{d}\n") out.write("\n") out.write("directions\n") - out.write("{0} # ndirections\n".format(self.ndir)) + out.write(f"{self.ndir} # ndirections\n") out.write("# id, coords...\n") for i, dir in enumerate(self.dirs): - out.write("{0} ".format(i)) + out.write(f"{i} ") for x in dir: - out.write("{0} ".format(x)) + out.write(f"{x} ") out.write("\n") out.write("\n") out.write("sites\n") - out.write("{0} # nsites\n".format(self.nsites)) + out.write(f"{self.nsites} # nsites\n") out.write("# id, type, coord...\n") - for i, site in enumerate(self.sites): - out.write("{0}\n".format(site)) + for site in self.sites: + out.write(f"{site}\n") out.write("\n") out.write("interactions\n") - out.write("{0} # nints\n".format(self.nints)) + out.write(f"{self.nints} # nints\n") out.write("# id, type, nbody, sites..., edge_flag, direction\n") - for i, inter in enumerate(self.ints): - out.write("{0}\n".format(inter)) + for inter in self.ints: + out.write(f"{inter}\n") + out.close() + # end of save_dat - def load_dat(self, inp): - if type(inp) is str: - with open(inp) as f: - self.load_dat(f) - return + def load_dat(self, filename: util.Filename) -> None: + inp = open(filename) state = "waiting" count = 0 @@ -247,35 +324,44 @@ def load_dat(self, inp): continue if state == "waiting" and body.lower() == "lattice": state = "dim" + continue elif state == "dim": self.dim = int(body) - self.size = None - self.bc = [None] * self.dim self.latvec = np.eye(self.dim) state = "size" + continue elif state == "size": - self.size = list(map(int, body.split())) - if len(self.size) != self.dim: - if len(self.size) < self.dim: - ERROR("too few elements ({0})".format(body)) + elem = body.split() + if len(elem) != self.dim: + if len(elem) < self.dim: + util.ERROR(f"too few elements ({body})") else: - ERROR("too many elements ({0})".format(body)) + util.ERROR(f"too many elements ({body})") + self.size = list(map(int, elem)) state = "bc" + continue elif state == "bc": elem = body.split() if len(elem) != self.dim: if len(elem) < self.dim: - ERROR("too few elements ({0})".format(body)) + util.ERROR(f"too few elements ({body})") else: - ERROR("too many elements ({0})".format(body)) - self.bc[:] = list(map(int, elem)) + util.ERROR(f"too many elements ({body})") + self.bc = list(map(bool, elem)) state = "latvec" + continue elif state == "latvec": - self.load_latvec(body, count) + self._load_latvec(body, count) count += 1 if count == self.dim: count = 0 break + # end of for line in inp: + + # working vars + dirs = [] + sites = [] + ints = [] # search for others inp.seek(0) @@ -295,141 +381,151 @@ def load_dat(self, inp): state = "nsites" elif body.lower() == "interactions": state = "nints" + continue elif state == "name": self.name = body state = "waiting" + continue elif state == "ndir": self.ndir = int(body) - self.dirs = [None for i in range(self.ndir)] + dirs = [None for i in range(self.ndir)] state = "dirs" count = 0 + continue elif state == "dirs": - self.load_direction(body, count) + self._load_direction(dirs, body) count += 1 if count == self.ndir: + self.dirs = typing.cast(List[np.ndarray], dirs) count = 0 state = "waiting" + continue elif state == "nsites": self.nsites = int(body) state = "sites" - self.sites = [None for i in range(self.nsites)] + sites = [None for i in range(self.nsites)] + continue elif state == "sites": - self.load_site(body, count) + self._load_site(sites, body) count += 1 if count == self.nsites: + self.sites = typing.cast(List[Site], sites) state = "waiting" count = 0 + continue elif state == "nints": self.nints = int(body) state = "ints" - self.ints = [None for i in range(self.nints)] + ints = [None for i in range(self.nints)] + continue elif state == "ints": - self.load_int(body, count) + self._load_int(ints, body) count += 1 if count == self.nints: + self.ints = typing.cast(List[Interaction], ints) state = "waiting" count = 0 - # end of for line in inp: - + inp.close() self.check_all() - self.update() + self._update() - def load_latvec(self, body, count): + def _load_latvec(self, body, count): elem = body.split() if len(elem) != self.dim + 1: if len(elem) < self.dim + 1: - ERROR("too few elements ({0})".format(body)) + util.ERROR("too few elements ({0})".format(body)) else: - ERROR("too many elements ({0})".format(body)) + util.ERROR("too many elements ({0})".format(body)) self.latvec[:, int(elem[0])] = list(map(float, elem[1:])) - def load_direction(self, body, count): + def _load_direction(self, dirs:List, body: str) -> None: elem = body.split() if len(elem) != self.dim + 1: if len(elem) < self.dim + 1: - ERROR("too few elements ({0})".format(body)) + util.ERROR(f"too few elements ({body})") else: - ERROR("too many elements ({0})".format(body)) - self.dirs[int(elem[0])] = list(map(float, elem[1:])) + util.ERROR(f"too many elements ({body})") + dirs[int(elem[0])] = np.array(list(map(float, elem[1:]))) - def load_site(self, body, count): + def _load_site(self, sites:List, body: str) -> None: elem = body.split() if len(elem) != self.dim + 2: if len(elem) < self.dim + 2: - ERROR("too few elements ({0})".format(body)) + util.ERROR(f"too few elements ({body})") else: - ERROR("too many elements ({0})".format(body)) - self.sites[int(elem[0])] = Site( + util.ERROR(f"too many elements ({body})") + sites[int(elem[0])] = Site( site_id=int(elem[0]), site_type=int(elem[1]), - coordinate=list(map(float, elem[2:])), + coordinate=np.array(list(map(float, elem[2:]))), ) - def load_int(self, body, count): + def _load_int(self, ints:List, body: str) -> None: elem = body.split() nbody = int(elem[2]) if len(elem) != nbody + 5: if len(elem) < nbody + 5: - ERROR("too few elements ({0})".format(body)) + util.ERROR(f"too few elements ({body})") else: - ERROR("too many elements ({0})".format(body)) - self.ints[int(elem[0])] = Interaction( + util.ERROR(f"too many elements ({body})") + ints[int(elem[0])] = Interaction( int_id=int(elem[0]), int_type=int(elem[1]), nbody=nbody, site_indices=list(map(int, elem[3 : (3 + nbody)])), - edge_flag=int(elem[-2]), + edge=int(elem[-2]), direction=int(elem[-1]), ) def check_all(self): - self.check_latvec() - self.check_dirs() - self.check_sites() - self.check_ints() + num_errors = 0 + num_errors += self._check_latvec() + num_errors += self._check_dirs() + num_errors += self._check_sites() + num_errors += self._check_ints() + if num_errors > 0: + sys.exit(1) - def check_latvec(self): - invalid = False + def _check_latvec(self) -> int: + num_errors = 0 for i in range(self.dim): if self.latvec[i] is None: - ERROR("lattice vector {0} is not defined".format(i), True) - invalid = True - if invalid: - sys.exit(1) + util.ERROR(f"lattice vector {i} is not defined", to_be_continued=True) + num_errors += 1 + return num_errors - def check_dirs(self): - invalid = False + def _check_dirs(self) -> int: + num_errors = 0 for i in range(self.ndir): if self.dirs[i] is None: - ERROR("bond direction {0} is not defined".format(i), True) - invalid = True - if invalid: - sys.exit(1) + util.ERROR(f"bond direction {i} is not defined", to_be_continued=True) + num_errors += 1 + return num_errors - def check_sites(self): - invalid = False + def _check_sites(self) -> int: + num_errors = 0 for i in range(self.nsites): if self.sites[i] is None: - ERROR("site {0} is not defined".format(i), True) - invalid = True - if invalid: - sys.exit(1) + util.ERROR(f"site {i} is not defined", to_be_continued=True) + num_errors += 1 + return num_errors - def check_ints(self): - invalid = False + def _check_ints(self) -> int: + num_errors = 0 for i in range(self.nints): if self.ints[i] is None: - ERROR("interaction {0} is not defined".format(i), True) - invalid = True - if invalid: - sys.exit(1) + util.ERROR(f"interaction {i} is not defined", to_be_continued=True) + num_errors += 1 + return num_errors + + def _update(self) -> None: + """ generate vertices""" - def update(self): self.nstypes = 0 - for s in self.sites: - s.z = 0 - self.nstypes = max(self.nstypes, s.stype) + for site in self.sites: + site.z = 0 + self.nstypes = max(self.nstypes, site.stype) self.nstypes += 1 self.nedges = 0 @@ -437,17 +533,15 @@ def update(self): if inter.edge == 1: inter.edge = self.nedges self.nedges += 1 - pass else: inter.edge = -1 - pass for s in inter.sites: self.sites[s].z += 1 self.nitypes = 0 self.nvtypes = 0 self.vertices = [] - vmap = {} + vmap: Dict = {} for inter in self.ints: itype = inter.itype self.nitypes = max(self.nitypes, itype) @@ -457,12 +551,14 @@ def update(self): inter.itype = vmap[(itype, stypes, zs)] else: inter.itype = vmap[(itype, stypes, zs)] = self.nvtypes - self.vertices.append(Vertex(self.nvtypes, itype, stypes, zs)) + self.vertices.append( + Vertex(self.nvtypes, itype, list(stypes), list(zs)) + ) self.nvtypes += 1 self.nitypes += 1 - def write_xml(self, filename): - with codecs.open(filename, "w", "utf_8") as f: + def write_xml(self, filename: util.Filename) -> None: + with open(filename, "w", encoding="utf_8") as f: f.write('\n') f.write("\n") @@ -517,41 +613,45 @@ def write_xml(self, filename): f.write(tagged("Direction", chain([i], bond))) f.write("\n") + # end of write_xml - def write_gnuplot(self, filename): + def write_gnuplot(self, filename: util.Filename) -> None: with open(filename, "w") as f: for st in range(self.nstypes): - f.write("$SITES_{0} << EOD\n".format(st)) + f.write(f"$SITES_{st} << EOD\n") for site in self.sites: if site.stype != st: continue v = np.dot(self.latvec, site.coord) for x in v: - f.write("{0} ".format(x)) - f.write("{0}\n".format(site.stype)) + f.write(f"{x} ") + f.write(f"{site.stype}\n") f.write("EOD\n") for bt in range(self.nitypes): - f.write("$BONDS_{0} << EOD\n".format(bt)) + f.write(f"$BONDS_{bt} << EOD\n") for bond in self.ints: if bond.nbody != 2: continue if bond.itype_org != bt: continue - v = np.dot(self.latvec, np.array(self.sites[bond.sites[0]].coord, dtype=float)) + v = np.dot( + self.latvec, + np.array(self.sites[bond.sites[0]].coord, dtype=float), + ) for x in v: - f.write("{0} ".format(x)) + f.write(f"{x} ") f.write("\n") v += np.dot(self.latvec, np.array(self.dirs[bond.dir], dtype=float)) for x in v: - f.write("{0} ".format(x)) + f.write(f"{x} ") f.write("\n\n") f.write("EOD\n") f.write("plot ") for st in range(self.nstypes): - f.write('$SITES_{0} w p pt {1} ps 2 t "" , \\\n'.format(st, st + 4)) + f.write(f'$SITES_{st} w p pt {st+4} ps 2 t "" , \\\n') for bt in range(self.nitypes): - f.write('$BONDS_{0} w l lw 2 lt {1} t "" , \\\n'.format(bt, bt + 1)) - f.write('\n') - f.write('pause -1\n') + f.write(f'$BONDS_{bt} w l lw 2 lt {bt+1} t "" , \\\n') + f.write("\n") + f.write("pause -1\n") diff --git a/tool/dsqss/lattice_factory/honeycomb.py b/tool/dsqss/lattice_factory/honeycomb.py index 51138d03..ea0c3fb5 100644 --- a/tool/dsqss/lattice_factory/honeycomb.py +++ b/tool/dsqss/lattice_factory/honeycomb.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,22 +14,24 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +from typing import Dict, MutableMapping + import numpy as np from scipy.special import cosdg, sindg -from ..util import get_as_list +import dsqss.util -def generate(param): +def generate(param: MutableMapping) -> Dict: dim = 2 - L = get_as_list(param, "L", extendto=dim) - bc = get_as_list(param, "bc", default=True, extendto=dim) + L = dsqss.util.get_as_list(param, "L", extendto=dim) + bc = dsqss.util.get_as_list(param, "bc", default=True, extendto=dim) basis = np.array([[1.0, 0.0], [cosdg(120), sindg(120)]]) sites = [ - {"siteid": 0, "type": 0, "coord": [0.0] * dim}, - {"siteid": 1, "type": 0, "coord": [1.0/3, 2.0/3]}, - ] + {"siteid": 0, "type": 0, "coord": [0.0] * dim}, + {"siteid": 1, "type": 0, "coord": [1.0 / 3, 2.0 / 3]}, + ] bonds = [ { "bondid": 0, diff --git a/tool/dsqss/lattice_factory/hypercubic.py b/tool/dsqss/lattice_factory/hypercubic.py index e73d6c7a..8647e1f0 100644 --- a/tool/dsqss/lattice_factory/hypercubic.py +++ b/tool/dsqss/lattice_factory/hypercubic.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -16,13 +16,13 @@ import numpy as np -from ..util import get_as_list +import dsqss.util def generate(param): dim = param["dim"] - L = get_as_list(param, "L", extendto=dim) - bc = get_as_list(param, "bc", default=True, extendto=dim) + L = dsqss.util.get_as_list(param, "L", extendto=dim) + bc = dsqss.util.get_as_list(param, "bc", default=True, extendto=dim) basis = np.eye(dim) sites = [{"siteid": 0, "type": 0, "coord": [0.0] * dim}] diff --git a/tool/dsqss/lattice_factory/kagome.py b/tool/dsqss/lattice_factory/kagome.py index 2d371f24..fc4f53a8 100644 --- a/tool/dsqss/lattice_factory/kagome.py +++ b/tool/dsqss/lattice_factory/kagome.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -15,22 +15,22 @@ # along with this program. If not, see . import numpy as np -from scipy.special import cosdg, sindg +from scipy import special -from ..util import get_as_list +import dsqss.util def generate(param): dim = 2 - L = get_as_list(param, "L", extendto=dim) - bc = get_as_list(param, "bc", default=True, extendto=dim) - basis = np.array([[1.0, 0.0], [cosdg(120), sindg(120)]]) + L = dsqss.util.get_as_list(param, "L", extendto=dim) + bc = dsqss.util.get_as_list(param, "bc", default=True, extendto=dim) + basis = np.array([[1.0, 0.0], [special.cosdg(120), special.sindg(120)]]) sites = [ - {"siteid": 0, "type": 0, "coord": [0.0, 0.0]}, - {"siteid": 1, "type": 0, "coord": [0.5, 0.0]}, - {"siteid": 2, "type": 0, "coord": [0.5, 0.5]}, - ] + {"siteid": 0, "type": 0, "coord": [0.0, 0.0]}, + {"siteid": 1, "type": 0, "coord": [0.5, 0.0]}, + {"siteid": 2, "type": 0, "coord": [0.5, 0.5]}, + ] bonds = [ { "bondid": 0, diff --git a/tool/dsqss/lattice_factory/triangular.py b/tool/dsqss/lattice_factory/triangular.py index 37d76b13..9ed73728 100644 --- a/tool/dsqss/lattice_factory/triangular.py +++ b/tool/dsqss/lattice_factory/triangular.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,17 +14,16 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import numpy as np -from scipy.special import cosdg, sindg +from scipy import special -from ..util import get_as_list +import dsqss.util def generate(param): dim = 2 - L = get_as_list(param, "L", extendto=dim) - bc = get_as_list(param, "bc", default=True, extendto=dim) - basis = [[1.0, 0.0], [cosdg(120), sindg(120)]] + L = dsqss.util.get_as_list(param, "L", extendto=dim) + bc = dsqss.util.get_as_list(param, "bc", default=True, extendto=dim) + basis = [[1.0, 0.0], [special.cosdg(120), special.sindg(120)]] sites = [{"siteid": 0, "type": 0, "coord": [0.0] * dim}] bonds = [ diff --git a/tool/dsqss/parameter.py b/tool/dsqss/parameter.py index d9733fd7..d4c49bbe 100644 --- a/tool/dsqss/parameter.py +++ b/tool/dsqss/parameter.py @@ -6,7 +6,7 @@ # 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. +# (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 @@ -16,15 +16,15 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from __future__ import print_function +from typing import MutableMapping -import codecs from copy import deepcopy -from .util import ERROR, INFO +import dsqss +import dsqss.util -def set_default_values(param): +def set_default_values(param: MutableMapping): for name, val in ( ("npre", 1000), ("ntherm", 1000), @@ -52,12 +52,12 @@ def set_default_values(param): def check_mandatories(param): for name in ("beta",): if name not in param: - ERROR("parameter {0} is not specified.".format(name)) + dsqss.util.ERROR(f"parameter {name} is not specified.") class Parameter(dict): def __init__(self, param): - super(Parameter, self).__init__() + super().__init__() if "parameter" in param: p = deepcopy(param["parameter"]) else: @@ -66,10 +66,11 @@ def __init__(self, param): check_mandatories(p) self.update(p) - def write_param(self, pfile): - with codecs.open(pfile, "w", "utf-8") as f: - for key in sorted(self.keys()): - f.write("{0} = {1}\n".format(key, self.get(key))) + def write_param(self, pfile: dsqss.util.Filename): + with open(pfile, "w", encoding="utf-8") as f: + for k in sorted(self.keys()): + v = self.get(k) + f.write(f"{k} = {v}\n") def main(): @@ -91,10 +92,11 @@ def main(): parser.add_argument( "-o", "--output", dest="pfile", default="param.in", help="Parameter file" ) + parser.add_argument("--version", action="version", version=dsqss.__version__) args = parser.parse_args() if args.input is sys.stdin: - INFO("waiting for standard input...") + dsqss.util.INFO("waiting for standard input...") p = Parameter(toml.load(args.input)) p.write_param(args.pfile) diff --git a/tool/dsqss/pmwa_pre.py b/tool/dsqss/pmwa_pre.py index 22317df6..8d38652f 100644 --- a/tool/dsqss/pmwa_pre.py +++ b/tool/dsqss/pmwa_pre.py @@ -23,6 +23,7 @@ import sys import warnings +from . import __version__ from .read_keyvalues import read_keyvalues from .util import ERROR @@ -267,6 +268,7 @@ def main(): action="store_true", help="Read from standard input (overrides --input)", ) + parser.add_argument("--version", action="version", version=__version__) args = parser.parse_args() if args.stdin: diff --git a/tool/dsqss/prob_kernel.py b/tool/dsqss/prob_kernel.py index c4d80409..48808163 100644 --- a/tool/dsqss/prob_kernel.py +++ b/tool/dsqss/prob_kernel.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,28 +14,36 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from copy import deepcopy +from typing import Sequence, Protocol import numpy as np -def heat_bath(weights, cutoff=1e-10): +class KernelCallBack(Protocol): + def __call__(self, weights: Sequence[float], cutoff: float = ...) -> np.ndarray: + ... + + +def heat_bath(weights: Sequence[float], cutoff: float = 1e-10) -> np.ndarray: """ return an array W, where W[i,j] is the probability of a transition i=>j calculated by the heat bath algorithm """ + ws = np.array(weights) ws /= ws.sum() - return np.array([ws for i in range(len(ws))]) + return np.array([ws for _ in range(len(ws))]) -def metropolis(weights, cutoff=1e-10): +def metropolis(weights: Sequence[float], cutoff: float = 1e-10) -> np.ndarray: """ return an array W, where W[i,j] is the probability of a transition i=>j - calculated by the Metropolice algorithm + (including sugestion probability T(i=>j)) + calculated by the Metropolis algorithm """ + N = len(weights) n = np.count_nonzero(weights) @@ -56,17 +64,18 @@ def metropolis(weights, cutoff=1e-10): return W -def suwa_todo(weights, cutoff=1e-10): +def suwa_todo(weights: Sequence[float], cutoff: float = 1e-10) -> np.ndarray: """ return an array W, where W[i,j] is the probability of a transition i=>j calculated by the Suwa-Todo algorithm (H. Suwa and S. Todo, Phys. Rev. Lett. 105, 120603 (2010)) """ + N = len(weights) W = np.zeros((N, N)) indices = np.argsort(weights, kind="mergesort")[::-1] - target = deepcopy(weights) + target = [w for w in weights] for i in range(N): s = weights[indices[i]] j = 1 @@ -87,7 +96,7 @@ def suwa_todo(weights, cutoff=1e-10): return W / W.sum(axis=1).reshape(-1, 1) -def reversible_suwa_todo(weights, cutoff=1e-10): +def reversible_suwa_todo(weights: Sequence[float], cutoff: float = 1e-10) -> np.ndarray: """ return an array W, where W[i,j] is the probability of a transition i=>j @@ -95,37 +104,38 @@ def reversible_suwa_todo(weights, cutoff=1e-10): under the detailed balance condition (H. Suwa and S. Todo, arXiv:1106.3562) """ - N = len(weights) - indices = np.argsort(weights, kind="mergesort")[::-1] - weights = np.array(weights) + ws = np.array(weights) + N = len(ws) + indices = np.argsort(ws, kind="mergesort")[::-1] + ws = np.array(ws) W = np.zeros([N, N]) for i in range(N): - I = indices[i] - W[I, I] = weights[I] + ii = indices[i] + W[ii, ii] = ws[ii] - S3 = sum(weights[indices[2:]]) - wdiff = weights[indices[0]] - weights[indices[1]] + S3 = sum(ws[indices[2:]]) + wdiff = ws[indices[0]] - ws[indices[1]] if wdiff >= S3: for i in range(1, N): - I = indices[i] - rst_swap(indices[0], I, weights[I], W) + ii = indices[i] + rst_swap(indices[0], ii, ws[ii], W) else: w3 = wdiff / S3 for i in range(2, N): - I = indices[i] - v = w3 * weights[I] - rst_swap(0, I, v, W) + ii = indices[i] + v = w3 * ws[ii] + rst_swap(0, ii, v, W) for j in range(N - 1, 0, -1): - J = indices[j] - v = W[J, J] / j + jj = indices[j] + v = W[jj, jj] / j for k in range(j - 1, -1, -1): - rst_swap(J, indices[k], v, W) + rst_swap(jj, indices[k], v, W) W /= W.sum(axis=1).reshape(-1, 1) W[W < cutoff] = 0.0 return W / W.sum(axis=1).reshape(-1, 1) -def rst_swap(i, j, w, W): +def rst_swap(i: int, j: int, w: np.ndarray, W: np.ndarray) -> None: W[i, i] -= w W[i, j] += w W[j, i] += w diff --git a/tool/dsqss/result.py b/tool/dsqss/result.py index a6e821aa..379d27ee 100644 --- a/tool/dsqss/result.py +++ b/tool/dsqss/result.py @@ -14,17 +14,20 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import dsqss +import dsqss.util + class Result: def __init__(self, mean, err): self.mean = mean self.err = err def to_str(self, delim="+/-"): - return "{}{}{}".format(self.mean, delim, self.err) + return f"{self.mean}{delim}{self.err}" class Results: - def __init__(self, resfile): + def __init__(self, resfile: dsqss.util.Filename): self.result = {} with open(resfile) as f: for line in f: diff --git a/tool/dsqss/std_lattice.py b/tool/dsqss/std_lattice.py index 82c09fb8..5de08906 100644 --- a/tool/dsqss/std_lattice.py +++ b/tool/dsqss/std_lattice.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,15 +14,17 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -# from .hypercubic import HyperCubicLattice -import codecs +import typing +from typing import Dict, Any -from .lattice import Lattice -from .lattice_factory import honeycomb, hypercubic, kagome, triangular -from .util import ERROR +import sys +from dsqss import __version__ +from dsqss import lattice +from dsqss import util +from dsqss.lattice_factory import honeycomb, hypercubic, kagome, triangular -def std_lattice(param, ret_dict=False): +def std_lattice(param): if "lattice" in param and not isinstance(param["lattice"], str): param = param["lattice"] if "unitcell" in param: @@ -38,18 +40,16 @@ def std_lattice(param, ret_dict=False): elif latname == "kagome": latticedict = kagome.generate(param) else: - ERROR( + util.ERROR( 'Unknown lattice: param["lattice"] = {0}'.format( param["lattice"].lower() ) ) + sys.exit(1) - lat = Lattice() + lat = lattice.Lattice() lat.load_dict(latticedict) - if ret_dict: - return lat, latticedict - else: - return lat + return lat, latticedict def main(): @@ -72,28 +72,32 @@ def main(): parser.add_argument( "-g", "--gnuplot", dest="gnuplot", default="", help="Output Gnuplot filename" ) + parser.add_argument("--version", action="version", version=__version__) args = parser.parse_args() - istoml = True + inp: Dict[str, Any] try: - inp = toml.load(args.input) + inp = typing.cast(Dict, toml.load(args.input)) + istoml = True except toml.TomlDecodeError: istoml = False + inp = {} print("Input file seems not to be TOML file.") if istoml: - lat, lattice_dict = std_lattice(inp, True) + lat, lattice_dict = std_lattice(inp) else: - lat = Lattice(args.input) + lattice_dict: Dict = {} + lat = lattice.Lattice(args.input) if len(args.out) > 0: lat.save_dat(args.out) if len(args.toml) > 0: if istoml: - with codecs.open(args.toml, "w", "utf-8") as f: + with open(args.toml, "w", encoding="utf-8") as f: toml.dump(lattice_dict, f) else: - ERROR("Cannot generate a lattice TOML file from a lattice data file.") + util.ERROR("Cannot generate a lattice TOML file from a lattice data file.") if len(args.gnuplot) > 0: lat.write_gnuplot(args.gnuplot) diff --git a/tool/dsqss/std_model.py b/tool/dsqss/std_model.py index d5f22a6a..3002e1e2 100644 --- a/tool/dsqss/std_model.py +++ b/tool/dsqss/std_model.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,25 +14,31 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from .bosehubbard import BoseHubbard_hamiltonian -from .util import ERROR -from .xxz import XXZ_hamiltonian -from .hamiltonian import Hamiltonian +import sys +from dsqss import __version__ +from dsqss import util +from dsqss import bosehubbard +from dsqss import xxz +from dsqss import hamiltonian -def std_model(param): + +def std_model(param) -> hamiltonian.Hamiltonian: if "hamiltonian" in param: param = param["hamiltonian"] if "model" in param: + ham: hamiltonian.Hamiltonian if param["model"].lower() == "spin": - ham = XXZ_hamiltonian(param) + ham = xxz.XXZ_hamiltonian(param) elif param["model"].lower() == "boson": - ham = BoseHubbard_hamiltonian(param) + ham = bosehubbard.BoseHubbard_hamiltonian(param) else: - ERROR('Unknown model: param["model"] = {0}'.format(param["model"].lower())) + util.ERROR( + 'Unknown model: param["model"] = {0}'.format(param["model"].lower()) + ) + sys.exit(1) else: - ham = Hamiltonian(param) - + ham = hamiltonian.Hamiltonian(param) return ham @@ -50,6 +56,7 @@ def main(): parser.add_argument( "-o", "--output", dest="out", default="hamiltonian.toml", help="Output filename" ) + parser.add_argument("--version", action="version", version=__version__) args = parser.parse_args() inp = toml.load(args.input) diff --git a/tool/dsqss/util.py b/tool/dsqss/util.py index dd22337a..13d95546 100644 --- a/tool/dsqss/util.py +++ b/tool/dsqss/util.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,52 +14,67 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from __future__ import print_function +from typing import TextIO, Union, List, MutableMapping, Optional +import os import sys from copy import deepcopy +Filename = Union[str, bytes, os.PathLike] -def INFO(msg, linebreak=True, file=sys.stderr): - print("INFO: {0}\n".format(msg), file=file) + +def LOG_IMPL( + msg: str, typ: str, to_be_continued: bool, linebreak: bool, file: TextIO +) -> None: + print(f"{typ}: {msg}\n", file=file) if linebreak: print("\n", file=file) + if not to_be_continued: + sys.exit(1) -def WARN(msg, linebreak=True, file=sys.stderr): - print("WARN: {0}\n".format(msg), file=file) - if linebreak: - print("\n", file=file) +def INFO(msg: str, linebreak: bool = True, file: TextIO = sys.stderr) -> None: + LOG_IMPL(msg, "INFO", linebreak=linebreak, to_be_continued=True, file=file) -def ERROR(msg, to_be_continued=False, linebreak=True, file=sys.stderr): - print("ERROR: {0}".format(msg), file=file) - if linebreak: - print("\n", file=file) - if not to_be_continued: - sys.exit(1) +def WARN(msg: str, linebreak: bool = True, file: TextIO = sys.stderr) -> None: + LOG_IMPL(msg, "WARN", linebreak=linebreak, to_be_continued=True, file=file) + + +def ERROR( + msg: str, + to_be_continued: bool = False, + linebreak: bool = True, + file: TextIO = sys.stderr, +): + LOG_IMPL( + msg, "ERROR", linebreak=linebreak, to_be_continued=to_be_continued, file=file + ) + pass -def extend_list(lst, N): +def extend_list(lst: List, N: int) -> List: last = lst[-1] for i in range(N - len(lst)): lst.append(deepcopy(last)) return lst -def get_as_list(d, name, default=None, extendto=None): +def get_as_list( + d: MutableMapping, name, default=None, extendto: Optional[int] = None +) -> List: if default is None: v = d[name] else: v = d.get(name, default) - if type(v) is not list: + if not isinstance(v, list): v = [v] if extendto is not None: extend_list(v, extendto) return v -def index2coord(index, size): +def index2coord(index: int, size: List) -> List: D = len(size) i = 0 r = [0 for d in range(D)] @@ -70,7 +85,7 @@ def index2coord(index, size): return r -def coord2index(r, size): +def coord2index(r: List, size: List) -> int: index = 0 block = 1 for x, L in zip(r, size): @@ -79,8 +94,8 @@ def coord2index(r, size): return index -def tagged(tagname, values): - if not hasattr(type(values), "__iter__") or type(values) is str: +def tagged(tagname: str, values) -> str: + if not hasattr(type(values), "__iter__") or isinstance(values, str): values = [values] s = "<{0}>".format(tagname) for v in values: diff --git a/tool/dsqss/wavevector.py b/tool/dsqss/wavevector.py index dfe20756..6c7e4e9e 100644 --- a/tool/dsqss/wavevector.py +++ b/tool/dsqss/wavevector.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,40 +14,41 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import codecs -from itertools import product +import itertools import numpy as np import scipy.special as spsp -from .util import ERROR, get_as_list, tagged +import dsqss +import dsqss.util +import dsqss.lattice class Wavevector: + dim: int + nk: int + ks: np.ndarray + def __init__(self): pass def generate(self, param, size): self.dim = len(size) - steps = get_as_list(param, "ksteps", default=0, extendto=self.dim) + steps = dsqss.util.get_as_list(param, "ksteps", default=0, extendto=self.dim) for d in range(self.dim): if steps[d] == 0: - steps[d] = size[d]//2 + steps[d] = size[d] // 2 ks = [] self.nk = 1 for d in range(self.dim): ks.append(list(range(0, size[d] // 2 + 1, steps[d]))) self.nk *= len(ks[-1]) - self.ks = np.zeros([self.dim, self.nk], dtype=np.int) - for ik, k in enumerate(product(*ks)): + self.ks = np.zeros([self.dim, self.nk], dtype=np.int64) + for ik, k in enumerate(itertools.product(*ks)): self.ks[:, ik] = list(k) - def load(self, inp): - if type(inp) is str: - with open(inp) as f: - self.load(f) - return - + def load(self, filename: dsqss.util.Filename) -> None: + inp = open(filename, encoding="utf-8") self.dim = 0 state = "waiting" for line in inp: @@ -63,7 +64,7 @@ def load(self, inp): self.dim = int(body) break if not self.dim > 0: - ERROR("invalid dimension") + dsqss.util.ERROR("invalid dimension") state = "waiting" for line in inp: @@ -77,18 +78,16 @@ def load(self, inp): state = "nk" elif state == "nk": self.nk = int(body) - self.ks = np.zeros([self.dim, self.nk], dtype=np.int) + self.ks = np.zeros([self.dim, self.nk], dtype=np.int64) state = "kpoints" elif state == "kpoints": words = body.split() kid = int(words[0]) self.ks[:, kid] = list(map(int, words[1:])) + inp.close() - def save(self, out): - if type(out) is str: - with codecs.open(out, "w", "utf-8") as f: - self.save(f) - return + def save(self, filename: dsqss.util.Filename) -> None: + out = open(filename, "w", encoding="utf-8") out.write("dim\n{0}\n".format(self.dim)) out.write("\n") @@ -100,11 +99,13 @@ def save(self, out): for k in self.ks[:, ik]: out.write(" {0}".format(k)) out.write("\n") + out.close() - def write_xml(self, filename, lat): + def write_xml(self, filename: dsqss.util.Filename, lat: dsqss.lattice.Lattice): + tagged = dsqss.util.tagged if self.dim != lat.dim: - ERROR("dimension mismatch between wavevector and lattice") - with codecs.open(filename, "w", "utf-8") as f: + dsqss.util.ERROR("dimension mismatch between wavevector and lattice") + with open(filename, "w", encoding="utf-8") as f: f.write('\n') f.write("\n") f.write(tagged("Comment", lat.name)) @@ -149,6 +150,7 @@ def main(): help='system size (specified by space separated integers like "4 4"). ' + "This option ignores the setting in [lattice] section of input TOML file", ) + parser.add_argument("--version", action="version", version=dsqss.__version__) args = parser.parse_args() @@ -156,13 +158,13 @@ def main(): inp = toml.load(args.input) if "lattice" in inp: dim = inp["lattice"]["dim"] - Ls = get_as_list(inp["lattice"], "L", extendto=dim) + Ls = dsqss.util.get_as_list(inp["lattice"], "L", extendto=dim) if "kpoints" in inp: inp = inp["kpoints"] if args.size is not None: Ls = list(map(int, args.size.split())) if Ls is None: - ERROR("system size is not specified") + dsqss.util.ERROR("system size is not specified") wv = Wavevector() wv.generate(inp, Ls) wv.save(args.out) diff --git a/tool/dsqss/xxz.py b/tool/dsqss/xxz.py index 74711dd0..f21bfdc0 100644 --- a/tool/dsqss/xxz.py +++ b/tool/dsqss/xxz.py @@ -4,7 +4,7 @@ # 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. +# (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 @@ -14,22 +14,24 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from math import sqrt +from typing import MutableMapping -from .hamiltonian import Hamiltonian, Interaction, Site, append_matelem -from .util import extend_list, get_as_list +import numpy as np +from dsqss import util +from dsqss import hamiltonian -def Splus(S, m): - return sqrt((S - m) * (S + m + 1.0)) +def Splus(S: float, m: float) -> float: + return np.sqrt((S - m) * (S + m + 1.0)) -def Sminus(S, m): - return sqrt((S + m) * (S - m + 1.0)) +def Sminus(S: float, m: float) -> float: + return np.sqrt((S + m) * (S - m + 1.0)) -class SpinSite(Site): - def __init__(self, id, M, D, h): + +class SpinSite(hamiltonian.Site): + def __init__(self, id: int, M: int, D: float, h: float) -> None: """ M: 2S+1 D: Sz^2 @@ -37,29 +39,28 @@ def __init__(self, id, M, D, h): """ S = 0.5 * M NX = M + 1 - values= [] - elements = {} - sources = {} + values = [] + elements: hamiltonian.MatElems = {} + sources: hamiltonian.MatElems = {} for st in range(NX): m = st - 0.5 * M values.append(m) - append_matelem(elements, state=st, value=-h * m + D * m * m) + hamiltonian.append_matelem(elements, state=st, value=-h * m + D * m * m) if st > 0: # annihilator - append_matelem( + hamiltonian.append_matelem( sources, istate=st, fstate=st - 1, value=0.5 * Sminus(S, m) ) if st < M: # creator - append_matelem( + hamiltonian.append_matelem( sources, istate=st, fstate=st + 1, value=0.5 * Splus(S, m) ) - super(SpinSite, self).__init__(id=id, N=NX, values=values, - elements=elements, sources=sources) + super().__init__(id=id, N=NX, values=values, elements=elements, sources=sources) -class XXZBond(Interaction): - def __init__(self, id, M, z, x): +class XXZBond(hamiltonian.Interaction): + def __init__(self, id: int, M: int, z: float, x: float) -> None: """ M: 2S+1 z: Sz_1 Sz_2 @@ -73,45 +74,45 @@ def __init__(self, id, M, z, x): Sz = [i - 0.5 * M for i in range(nx)] Sp = [Splus(S, m) for m in Sz] Sm = [Sminus(S, m) for m in Sz] - elements = {} + elements: hamiltonian.MatElems = {} for i in range(nx): for j in range(nx): # diagonal w = -z * Sz[i] * Sz[j] if w != 0.0: - append_matelem(elements, state=[i, j], value=w) + hamiltonian.append_matelem(elements, state=[i, j], value=w) # offdiagnal w = -0.5 * x * Sp[i] * Sm[j] if w != 0.0: - append_matelem( + hamiltonian.append_matelem( elements, istate=[i, j], fstate=[i + 1, j - 1], value=w ) w = -0.5 * x * Sm[i] * Sp[j] if w != 0.0: - append_matelem( + hamiltonian.append_matelem( elements, istate=[i, j], fstate=[i - 1, j + 1], value=w ) - super(XXZBond, self).__init__(id=id, nbody=nbody, Ns=Ns, elements=elements) + super().__init__(id=id, nbody=nbody, Ns=Ns, elements=elements) -class XXZ_hamiltonian(Hamiltonian): - def __init__(self, param): +class XXZ_hamiltonian(hamiltonian.Hamiltonian): + def __init__(self, param: MutableMapping) -> None: M = param["M"] - Ds = get_as_list(param, "D", 0.0) - hs = get_as_list(param, "h", 0.0) - nstypes = max(len(Ds), len(hs)) - extend_list(Ds, nstypes) - extend_list(hs, nstypes) + Ds = util.get_as_list(param, "D", 0.0) + hs = util.get_as_list(param, "h", 0.0) + self.nstypes = max(len(Ds), len(hs)) + util.extend_list(Ds, self.nstypes) + util.extend_list(hs, self.nstypes) self.sites = [SpinSite(i, M, D, h) for i, (D, h) in enumerate(zip(Ds, hs))] - Jzs = get_as_list(param, "Jz", 0.0) - Jxys = get_as_list(param, "Jxy", 0.0) - nitypes = max(len(Jzs), len(Jxys)) - extend_list(Jzs, nitypes) - extend_list(Jxys, nitypes) + Jzs = util.get_as_list(param, "Jz", 0.0) + Jxys = util.get_as_list(param, "Jxy", 0.0) + self.nitypes = max(len(Jzs), len(Jxys)) + util.extend_list(Jzs, self.nitypes) + util.extend_list(Jxys, self.nitypes) self.interactions = [ XXZBond(i, M, Jz, Jx) for i, (Jz, Jx) in enumerate(zip(Jzs, Jxys)) ] - self.name = "S={0}/2 XXZ model".format(M) + self.name = f"S={M}/2 XXZ model" diff --git a/tool/pyproject.toml b/tool/pyproject.toml new file mode 100644 index 00000000..8134fa52 --- /dev/null +++ b/tool/pyproject.toml @@ -0,0 +1,34 @@ +[tool.poetry] +name = "dsqss" +version = "2.1-dev" +description = "dsqss inputfile generator" +authors = ["DSQSS developers "] +license = "GPL-3.0-or-later" + +repository = "https://github.com/issp-center-dev/DSQSS" + +packages = [ + { include = "dsqss" } + ] + +[tool.poetry.dependencies] +python = "^3.8" +# Users should install dependencies by themselves +# numpy = "^1.17" +# toml = ">= 0.10.0" +# scipy = "^1" + +[tool.poetry.dev-dependencies] + +[tool.poetry.scripts] +dla_hamgen = "dsqss.std_model:main" +dla_latgen = "dsqss.std_lattice:main" +dla_pgen = "dsqss.parameter:main" +dla_wvgen = "dsqss.wavevector:main" +dla_alg = "dsqss.dla_alg:main" +dla_pre = "dsqss.dla_pre:main" +pmwa_pre = "dsqss.pmwa_pre:main" + +[build-system] +requires = ["poetry-core>=1.0.0"] +build-backend = "poetry.core.masonry.api" diff --git a/tool/setup.py b/tool/setup.py deleted file mode 100644 index 1dda8876..00000000 --- a/tool/setup.py +++ /dev/null @@ -1,42 +0,0 @@ -# DSQSS (Discrete Space Quantum Systems Solver) -# Copyright (C) 2018- The University of Tokyo -# -# 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 . - -from setuptools import setup - -setup( - name="dsqss", - version="2.0.4", - description="DSQSS input files generator", - url="https://github.com/issp-center-dev/dsqss", - author="DSQSS developers", - author_email="dsqss-dev@issp.u-tokyo.ac.jp", - license="GPLv3", - packages=["dsqss", 'dsqss.lattice_factory'], - python_requires=">=2.7", - install_requires=["numpy", "scipy", "toml"], - entry_points={ - "console_scripts": [ - "dla_hamgen = dsqss.std_model:main", - "dla_latgen = dsqss.std_lattice:main", - "dla_pgen = dsqss.parameter:main", - "dla_wvgen = dsqss.wavevector:main", - "dla_alg = dsqss.dla_alg:main", - "dla_pre = dsqss.dla_pre:main", - "pmwa_pre = dsqss.pmwa_pre:main", - ] - }, - zip_safe=False, -)